Skip to content

Instantly share code, notes, and snippets.

@borgwang
Last active April 12, 2024 09:28
Show Gist options
  • Select an option

  • Save borgwang/4313e9375ef233c3b812f9f80f1af2bb to your computer and use it in GitHub Desktop.

Select an option

Save borgwang/4313e9375ef233c3b812f9f80f1af2bb to your computer and use it in GitHub Desktop.
quantile_regression
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Quantile Regression Tutorial"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### 特点\n",
"\n",
"- 一般回归模型对自变量与因变量均值进行回归,假定了预测的误差的方差是固定的(实际中有时并不是)\n",
"- 分位数回归对自变量与因变量的不同条件分位数进行回归\n",
"- 可以得到一个置信区间,通过不同的分位了解预测的方差\n",
"\n",
"#### 基本原理\n",
"\n",
"假定因变量真实值为 $ y=(y_1, y_2, ..., y_n) $,目标值为 $ \\hat{y} = (\\hat{y_1}, \\hat{y_2}, ..., \\hat{y_n}) $\n",
"\n",
"分位数回归的损失函数是:\n",
"\n",
"$$ L_r(y, \\hat{y}) = (1-r)\\frac{1}{N}\\sum^{i}_{\\hat{y_i}\\geq y_i}(\\hat{y_i}-y_i) + r\\frac{1}{N}\\sum^{i}_{\\hat{y_i}< y_i}(y_i-\\hat{y_i}) $$\n",
"\n",
"其中 r 是分位数系数,这个损失函数是平均绝对误差的拓展,当 r=0.5 时退化成 LAD regression(Least absolute deviations)。\n",
"\n",
"Intuition:当 r > 0.5 时,损失函数对预测值偏小的数据惩罚更大(使 fitting 到曲线上移,值更大);当 r > 0.5 时,损失函数对预测值偏小的数据惩罚更大。即 r 控制了模型对于该分位系数下两类数据的不同惩罚程度。\n",
"\n",
"#### 梯度\n",
"\n",
"使用最简单的线性模型 y = wx + b ,上述损失函数对应 w,b 的梯度为\n",
"\n",
"$$ \\frac{dL_{w,b}}{dw} = (1-r)\\frac{1}{N}\\sum^{i}_{wx_i+b\\geq y_i}x_i + r\\frac{1}{N}\\sum^{i}_{wx_i+b< y_i}(-x_i) $$\n",
"\n",
"$$ \\frac{dL_{w,b}}{db} = (1-r)\\frac{1}{N}\\sum^{i}_{wx_i+b\\geq y_i}\\times1 + r\\frac{1}{N}\\sum^{i}_{wx_i+b< y_i}\\times(-1) $$\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PathCollection at 0x7fe26850d470>"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAG11JREFUeJzt3X9sZeV95/H3F+MUU7o1FCcaTMjQNpq0TcR440V0R9qlk6jDLjS4SUObTSNSoU5W2uzCFk0xUVVYaVdxRUmyUqtI00DDahEFwaxBkC1FzNAoVGHjwUMGOqCk4UfGzDKugtMQvMRjvvuHr+GO55x7z70+5zzPee7nJY3G9/jce54zY3/Pc77P93mOuTsiItJ8p4VugIiIlEMBXUQkEQroIiKJUEAXEUmEArqISCIU0EVEEqGALiKSCAV0EZFEKKCLiCTi9DoPdu655/rWrVvrPKSISOMdPHjwH919rNt+tQb0rVu3Mjc3V+chRUQaz8xeLLKfUi4iIolQQBcRSYQCuohIIhTQRUQSoYAuIpKIWqtcRESqNDu/wC0PP8fLS8ucNzrCnl3bmJoYD92s2iigi0gSZucXuHHfYZZXVgFYWFrmxn2HAQYmqCvlIiJJuOXh594K5uuWV1a55eHnArWofgroIpKEl5eWe9qeIgV0EUnCeaMjPW1PkQK6iCRhz65tjAwPnbRtZHiIPbu2BWpR/TQoKiJJWB/4VJWLiEgCpibGByqAb6SUi4hIIhTQRUQSUTigm9mQmc2b2YOt1xea2RNm9h0zu9vM3lFdM0VEpJteeujXAkfaXv8J8EV3fy/wKnBNmQ0TEZHeFAroZnY+cDnwldZrA3YC97Z2uQOYqqKBIiJSTNEe+peAPwTebL3+OWDJ3U+0Xh8FBndoWUQkAl0DupldARx394PtmzN29Zz37zazOTObW1xc7LOZIiLSTZEe+g7gI2b2AvBXrKVavgSMmtl6Hfv5wMtZb3b3ve4+6e6TY2NdH1otIiJ96hrQ3f1Gdz/f3bcCvwPsd/dPAgeA32rtdjVwf2WtFBGRrjZTh34D8Adm9l3Wcuq3ldMkERHpR09T/939MeCx1tffAy4uv0kiItIPzRQVEUmEArqISCIU0EVEEqGALiKSCK2HLiJSkdn5hVofuKGALiJSgdn5BW7cd5jllVUAFpaWuXHfYYDKgrpSLiIiFbjl4efeCubrlldWueXh5yo7pgK6iEgFXl5a7ml7GRTQRUQqcN7oSE/by6CALiJSgT27tjEyPHTStpHhIfbs2lbZMTUoKiJSgfWBT1W5iIgkYGpivNIAvpFSLiIiiVBAFxFJhFIuIjIQ6p61GYICuogkL8SszRCUchGR5IWYtRmCeugikrw6Z22GTO2ohy4iyatr1uZ6amdhaRnn7dTO7PxCqcfJo4AuIsmra9Zm6NSOUi4ikry6Zm2GWJCrXdeAbmZnAF8Hfqq1/73ufpOZfRX418APW7t+2t0PVdVQEZHNqGPW5nmjIyxkBO8qF+RqV6SH/gaw091fM7Nh4Btm9r9b39vj7vdW1zwRkTVNqCPfs2vbSeWRUP2CXO26BnR3d+C11svh1h+vslEiIu2aUkceYkGudrYWr7vsZDYEHAR+Efhzd7+hlXL5VdZ68I8C0+7+RqfPmZyc9Lm5uU03WkQGy46Z/ZmpjPHRER6f3hmgRfUys4PuPtltv0JVLu6+6u7bgfOBi83s/cCNwPuAfwGcA9yQ05DdZjZnZnOLi4uFT0BEZF3owcam6Kls0d2XgMeAy9z9mK95A/hL4OKc9+x190l3nxwbG9t0g0Vk8IR4+k8TdQ3oZjZmZqOtr0eADwPPmtmW1jYDpoCnq2yoiAyuEE//aaIiVS5bgDtaefTTgHvc/UEz229mY4ABh4B/X2E7RWSAhR5sbIpCg6Jl0aCoiEjvSh0UFRGR+Cmgi4gkQgFdRCQRCugiIolQQBcRSYSWzxURKVmohcQU0EVEShRyITGlXEREShTyqUXqoYtIcE1Y67yokAuJqYcuIkGFfrBy2UIuJKaALiJBhX6wctlCLiSmlIuIBJXaWuchFxJTQBeRoEI/WLkKdTyQOotSLiISlNY6L4966CISlNY6L48CuogEFypFkRoFdBEJLqU69JAU0EUkqJBT5VOjQVERCSq1OvSQFNBFJKjU6tBDUkAXkaBCTpVPTdeAbmZnmNn/MbOnzOwZM/svre0XmtkTZvYdM7vbzN5RfXNFJDWqQy9PkR76G8BOd78I2A5cZmaXAH8CfNHd3wu8ClxTXTNFJFVTE+N8/qMfYHx0BAPGR0f4/Ec/oAHRPnStcnF3B15rvRxu/XFgJ/DvWtvvAG4Gvlx+E0UkdapDL0ehHLqZDZnZIeA48AjwD8CSu59o7XIUyPzfMLPdZjZnZnOLi4tltFlERDIUqkN391Vgu5mNAv8L+KWs3XLeuxfYCzA5OZm5j4gMLk0qKk9PE4vcfcnMHgMuAUbN7PRWL/184OUK2iciJYspgGpSUbmKVLmMtXrmmNkI8GHgCHAA+K3WblcD91fVSBEpR2xPB9KkonIVyaFvAQ6Y2beBbwGPuPuDwA3AH5jZd4GfA26rrpkiUobYAqgmFZWrSJXLt4GJjO3fAy6uolEiUo3YAmiKD7cISTNFRQZIbLMyNamoXFptUaQiMQ0+rtuza9tJg5BQXgDt53z1cItyKaCLVKCq6o3NXiSqCqCbOV9NKiqPrU0Ercfk5KTPzc3VdjyRUHbM7M/MDY+PjvD49M6+PnNj0IS13vXHPjjOgWcXg/ZwqzhfeZuZHXT3yW77qYcuUoEqBh/zKlTu/OZLb83qC1XHHdtg66DSoKhIBaoYfMwLjhvvsUOUIcY22DqoFNBFKlBF9UYvwbHunnFK1Sqz8wvsmNnPhdMPsWNmf7BJV/1QQBepQBVLwmYFTcvZt+6ecSpL4MY2k7ZXGhQVaZCNVS6/9r4x7ju4cMpAaRODaQxiHdzVoKgIcdaCb0ZWid/ke85J6hxDavrgrgK6JGtQVvJrah13jBfbpi9FoBy6JCu2hajkbb3kquscpGz64K4CuiSr6bfPKSt6sa17kLLpg7tKuUiymn77nLKiF9tOgb+qINvUFBYooEuD9JpzrXIhqhBizDn362dHhllaXjll+8aLre6yeqOALo3QzwBnSiv5pTTAOzu/wI9/cuKU7cOn2SkXW91l9UYBXRqh31vvJt8+twuReqjKLQ8/x8rqqfNfzjrj9FPOJbW7rKopoEsjDPqtd0rnn9fmpddPTcGkdJdVBwV0aYRBv/VO6fx7PZeQd1lNG7dQ2aI0QtPrgzeraeffqXa8KefSxHVdugZ0M3u3mR0wsyNm9oyZXdvafrOZLZjZodaff1t9c2VQNb0+eLOadP7dAmFTzqWJE9O6Ls5lZluALe7+pJn9DHAQmAKuAl5z9z8tejAtziXSvNv4XsW6wFWvLpx+6JS15mFthcvnZy6vtS2lLc7l7seAY62vf2RmR4B0fvpEahRz+WFZF5pUBnCbOG7RUw7dzLYCE8ATrU2fNbNvm9ntZnZ2yW0TSU6st/Fl5otTeXpRU3L97QoHdDM7C7gPuM7d/wn4MvALwHbWevC35rxvt5nNmdnc4uJiCU0Waa5Ye695F5rr73mq56DexECYpSm5/naFyhbNbJi1YH6nu+8DcPdX2r7/F8CDWe91973AXljLoW+2wSJNFuttfN4FZdW955RQSrXjTZuY1jWgm5kBtwFH3P0Lbdu3tPLrAL8JPF1NE0XS0W3mY6gB07wLDfQ3I7VpgTAVRXroO4BPAYfN7FBr2+eAT5jZdtYeOv4C8JlKWiiSkE6915ADplkXmnahU0Kxiq1iqUiVyzfIfhbt18pvjkj68nqvIddrWf/86+95itWMUubQKaEYxVixpJmiIpEIPWA6NTHOrVddlMSAZh1irFhSQBeJRAzlfk2s7Agl9AU4ixbnEolELEvFakCzmBgrltRDF4lE7L3jOh/W3AQx1turhy4SkVh7xzEOAIYWY729ArpIYLGVvmVJ6YlJZYrtAqyALhJQU3q+MQ4AVqUJF9g8yqGLBBRj6VuWGCpw6tDEh1q0U0CXZDRx0K4pPd8YBwCr0JQLbB6lXCQJoVMX/d6mx1j6liXGAcAqNOUCm0cBXZIQctBuMxeTWGrPi4htALAKTbnA5lHKRZIQsmd18wPP9H2bHnvt+aBpempJPXRJQqie1ez8AkvLK5nfK3oxaVLPt8kVIEU0PbWkgC5JCJW66NQLb8ptelGhxynq0qQL7EZKuUgSQqUuOvXCm3KbXlTTK0AGgXro0ghFbvVD9KzyUj1nnznc2F5enqZXgAwC9dAlejFM9sircc8bRLvpN36ltrbVZVAmFzWZArpEL/StfqcLSt2pnpCTp5peATIIlHKR4LqlU8q81e+nSqNbjXtdqZ7Qg5JNrwAZBAroElSRIFVWSWK/ATGW3HEMKx42uQJkECjlIkEVSaeUdavfb+omltxx2ReWJq59I511Dehm9m4zO2BmR8zsGTO7trX9HDN7xMy+0/r77OqbK6kpEqTKylPnHWthabljQIsld1zmhSWGgWYpX5GUywngend/0sx+BjhoZo8AnwYedfcZM5sGpoEbqmuqpKhoOqWMW/28Y0Hn9EssueMyJ0/FkL6R8nUN6O5+DDjW+vpHZnYEGAeuBC5t7XYH8BgK6NKjMoNUtwHPrGO16xTQQuSOs87n8x/9QCkXlljGBaRcPQ2KmtlWYAJ4AnhXK9jj7sfM7J2lt06S12vvNy9oFxnwbD9WXk89loCWdz6f/+gHeHx656Y/v+mrCko2c/diO5qdBfwt8N/cfZ+ZLbn7aNv3X3X3U/LoZrYb2A1wwQUXfPDFF18sp+UycDYGOVjrza/3WrMC1PjoSGYA3DGzv6f961Z1+zr9WyrlEh8zO+juk932K1TlYmbDwH3Ane6+r7X5FTPb0vr+FuB41nvdfa+7T7r75NjYWLHWi2TolPftNOCZNdAXy0BnnqpTIlq2N01dUy5mZsBtwBF3/0Lbtx4ArgZmWn/fX0kLJUn9TPDpFOQ6DXhmDXbGMtCZZXZ+gdPMWM24ey4zJaKa8vQUyaHvAD4FHDazQ61tn2MtkN9jZtcALwEfr6aJkpp+J/h0yvt2GvDMG+yMMaCt/9tkBfOY7iCqVsa666mv3Z6lSJXLNwDL+faHym2ODIJ+S+Y6VcSsv++6uw9lvjeWwc5usv5tAIbMBiYlUsYSB6GXSQhFM0Wldv3mh7vlfacmxhmPZFZnv/L+Dd50TzoQtStjMbbQC7qForVcpHb9lMxtvH3+4m9vzwxwTXrocpamlhOWmd4oY0B4UOvs1UOX2vVaYdLLNPWmV2/EXn2TpexlBMpY4iCW9Xfqph66VCav19ZrhUm32+dOx6hSv73STu+LufomT9nLCJRxl9X0O7V+FZ5YVIbJyUmfm5ur7XgSTpkTVy6cfoi8n9KR4aEgk2P6Pb/UJvTMzi/kDkQb8PzM5X1/rqpc3lZ0YpECupSm/Rcor466n5mOebMmh0o8Rq/6nckZ+wzVXmRdnNo18ZxiVepMUZFuNuZRswIt9DcolZdXLvMYvep30C2lwbq8EkuA4dOM139yQmut10wBXUrR6Ze7XT+DUnkDnSFLFPsddEtpsK7jRcjg1ddXtNZ6zRTQpRR50+7bbWZQampinMend/L8zOU8Pr2TqYnxoBUh/R67iVUsefIuQkNmrKyefPc0CDXgMVBAl1IMWd5kYiorHwxZotjvsZteVtkuxlTYoNOgaMPFMpK/dfqh3O+90Gelg8Qv6+ev16WMpbuig6KqQ2+wmNarGM+Z4ZiX55Y05NX8D2INeAyUcmmwmNarSCk3LJuTUlqpadRDb7BOJXB1p2LqmOEYS3pJuotxaeJBoIDeYHkLOf3syHCQVEyVv8QxpZdEYqWUS4PlpTnMiCYVU5aY0ksisVJAb7C8XOXS6yuZ+ze5bKzTM0M1G1FkjVIuDZeV5sgrG2vibMR1nZ4Z2j4bEZSCkcGlHnqCUqw4yTqnjZSCkUGnHnqCQq+pXUU1ysZzypsO1+S0kshmKaAnKlTZWFXVKBsvEj9+4wRLy6eOFTQ5rSSyWV0DupndDlwBHHf397e23Qz8PrDY2u1z7v61qhopccqb9p1VjXL9PU8B/QX1rIvE8JAxfJqx8ubbffWmp5VENqtIDv2rwGUZ27/o7ttbfxTMB0zecyTzBi5X3fteQjXrIrGy6px1xumajSjSpmsP3d2/bmZbq2+KdBLbLMm8nnjeU4TWv9/Pcybz8uJLr68w/8e/3tNniaRsM1UunzWzb5vZ7WZ2dt5OZrbbzObMbG5xcTFvN+mg7KeqlyEvyK66d6xGKbJu+kYpPRRCpEr9BvQvA78AbAeOAbfm7ejue9190t0nx8bG+jzcYItxlmReMF1PfeStj27Q84UoxTJMkSr0FdDd/RV3X3X3N4G/AC4ut1lxmp1fYMfM/tpnJvbyHMq62tgpyE5NjHPrVReRFdIder4QafU+kWL6Kls0sy3ufqz18jeBp8trUpxCLg6VN0tyYy+5zjZ2q3WfmhjnursPZb63n1pxrd4n0l2RssW7gEuBc83sKHATcKmZbWetw/UC8JkK2xiFTmmPqgPNnl3bCj0woO42dguyeQ+9UO5bpBpFqlw+kbH5tgraErVe0h5lKzrzM2QbsxS9EIlIObSWS0F5vcrRM4crP3bRksXYqkGU+xapl6b+F7Rn1zb23PsUK6sn11i/9v9OMDu/EMWDHWLsESv3HZfY5jNIudRDL2hqYpzh006t21h50ystH+ylZFE9YukkxvkMUi710AuanV/g9ZU3M79XZY6604Mdsu4M1COWPCEH9qUe6qEX1KkXXmWOutNnq3clvYht0FzKp4BeUKcf+ipz1J0e7BBqtmioCVayObENmkv5FNALyq1yGRmu9HZ1PS+ep+7elfKwzaUlFNKngF5Q3i/DFRdtqby3OjUxzngkvasY15WRYjRonj7znKVOqzA5Oelzc3O1Ha9sfzR7mLue+D6r7gyZccnPn82TL/3wpABnrE2fHS+5JGxj+SKsXVDq/oW8cPqhzMe/GfD8zOW1tUNkkJjZQXef7LbfwFa59FqPOzu/wH0HF95a63vVnb/7hx+cEtzWX2+sF99s/W/o54SuK7qujIjUbyADet5knbkXf8CBZxczA2ZWqqHbvU17KqKMRbNiKEmMcfKSiKwZyBx6Xh74zm++lDvY1+/g48tLy0nlnZWHFYnXQPbQ84Lzxh53+6SLvFRDN+eNjnScHLRjZn/jpl/HcKcgIqcayB56L/ne9WDcqR48z3oqotPxVPZXnOrfRTpLIqD3+oueFZyzH5j2dvBfTzXkPVpt4+e0pyK6XQyamn6pk+rfRbprfMola4DzursPcfMDz3DzR34lMzWQVTHya+8b476DCx0H+9bft3FQsFupYvvx8tI2mn7dmdYhEemu8QE96xcdYGl5pWMlSVYeePI953QtC+y3fHD9eDtm9qvsrw9ah0Sku8YH9E6/0L324IoO9m1mUFBlf/1R/btId40N6OsTdbrVgpfdg6tygpAePpBPF0KR7hoZ0P9o9jB3fvOlrsEcyu3B9fL0oE6yevhlfXaqYpkpKxKzrmu5mNntwBXAcXd/f2vbOcDdwFbgBeAqd3+128E2u5bL7PwCNz/wDEvLK4X2Hxke4mMfHM+d/dmrvPz3+OgIj0/v7Osz6/hsEWm2Mtdy+SrwZ8D/aNs2DTzq7jNmNt16fUM/De1mPQ2xsLT8VjVJEUNmfOyD4ydVrmy211vlwJwG/URks7rWobv714EfbNh8JXBH6+s7gKmS2wWcXHsMxYP5yPAQt151EQeeXSx1yn2VDwiI9eEDmswj0hz9Tix6l7sfA2j9/c7ymvS2vJLETton9JTd663yAQExPnxAk3lEmqXyQVEz2w3sBrjgggt6em/RwGvAJy+5gP86dfKTfcoudatyYC7GQT9N5hFpln4D+itmtsXdj5nZFuB43o7uvhfYC2uDor0cpMiCWGefOcxNv5E9I7SKUrcqF6aKbdEr5fVFmqXflMsDwNWtr68G7i+nOSfrtObK+OgIX/rt7cz/8a/nBkEt9bo5seb1RSRb1x66md0FXAqca2ZHgZuAGeAeM7sGeAn4eBWNKyMNEVuvt0k0mUekWfRMUelIs1dFwtMzRaUUusMRaQ4FdBEZSCnefSqgi8jASXXtpCSeWCQi0ouUHtzeTgFdRAZOqnMsFNBFZOCkOsdCAV1EBk6MayeVQYOiIjJwYlw7qQwK6CIykFKcY6GUi4hIIhTQRUQSoYAuIpIIBXQRkUQooIuIJKLW5XPNbBF4sc+3nwv8Y4nNaQqd9+AYxHMGnXcR73H3sW471RrQN8PM5oqsB5wanffgGMRzBp13mZ+plIuISCIU0EVEEtGkgL43dAMC0XkPjkE8Z9B5l6YxOXQREemsST10ERHpIPqAbmaXmdlzZvZdM5sO3Z46mNm7zeyAmR0xs2fM7NrQbaqTmQ2Z2byZPRi6LXUxs1Ezu9fMnm39v/9q6DbVwcz+c+tn/Gkzu8vMzgjdpiqY2e1mdtzMnm7bdo6ZPWJm32n9ffZmjxN1QDezIeDPgX8D/DLwCTP75bCtqsUJ4Hp3/yXgEuA/DMh5r7sWOBK6ETX778Bfu/v7gIsYgPM3s3HgPwGT7v5+YAj4nbCtqsxXgcs2bJsGHnX39wKPtl5vStQBHbgY+K67f8/dfwL8FXBl4DZVzt2PufuTra9/xNovd1rrfOYws/OBy4GvhG5LXczsnwH/CrgNwN1/4u5LYVtVm9OBETM7HTgTeDlweyrh7l8HfrBh85XAHa2v7wCmNnuc2AP6OPD9ttdHGZDAts7MtgITwBNhW1KbLwF/CLwZuiE1+nlgEfjLVqrpK2b206EbVTV3XwD+FHgJOAb80N3/JmyravUudz8Ga5044J2b/cDYA7plbBuYshwzOwu4D7jO3f8pdHuqZmZXAMfd/WDottTsdOCfA1929wngx5Rw+x27Vs74SuBC4Dzgp83sd8O2qtliD+hHgXe3vT6fRG/JNjKzYdaC+Z3uvi90e2qyA/iImb3AWnptp5n9z7BNqsVR4Ki7r9+F3ctagE/dh4Hn3X3R3VeAfcC/DNymOr1iZlsAWn8f3+wHxh7QvwW818wuNLN3sDZg8kDgNlXOzIy1fOoRd/9C6PbUxd1vdPfz3X0ra//X+909+R6bu/9f4Ptmtv6E4g8Bfx+wSXV5CbjEzM5s/cx/iAEYDG7zAHB16+urgfs3+4FRP1PU3U+Y2WeBh1kbAb/d3Z8J3Kw67AA+BRw2s0OtbZ9z968FbJNU6z8Cd7Y6Lt8Dfi9weyrn7k+Y2b3Ak6xVds2T6KxRM7sLuBQ418yOAjcBM8A9ZnYNaxe3j2/6OJopKiKShthTLiIiUpACuohIIhTQRUQSoYAuIpIIBXQRkUQooIuIJEIBXUQkEQroIiKJ+P88VW7ZxYsSVAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# generate data with varying variances\n",
"X = np.random.uniform(0, 10, (100, 1))\n",
"Y = 2 * X + 10 + np.random.normal(0, 0.8 * X)\n",
"\n",
"plt.scatter(X, Y)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fe2906a8208>]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+QXGWd7/H3N8Mow49lgEQNE9jJ3XUju6AJO1C4qVU3WDe4RBlxRdT1osWKlr/AtWISCoWt0kqAVcFSuTcrLFCyGAriQAUx6xIQN5ask0wgshHxSsAMuTAKIwJjnEy+94/unkxmzuk+3XN+9+dVlcrM6dN9nk5mvufp7/N9nsfcHRERKb45WTdARETioYAuIlISCugiIiWhgC4iUhIK6CIiJaGALiJSEgroIiIloYAuIlISCugiIiVxWJoXmzt3rvf29qZ5SRGRwtu2bduv3X1eo/NSDei9vb0MDg6meUkRkcIzsyejnKeUi4hISSigi4iUhAK6iEhJKKCLiJSEArqISEmkWuUiIpKkgaFhrtn8GE+PjnFCdxcrly+if0lP1s1KjQK6iJTCwNAwazbuZGx8AoDh0THWbNwJ0DZBXSkXESmFazY/NhnMa8bGJ7hm82MZtSh9CugiUgpPj441dbyMFNBFpBRO6O5q6ngZKaCLSCmsXL6Irs6OQ451dXawcvmijFqUPg2Kikgp1AY+VeUiIlIC/Ut62iqAT6eUi4hISSigi4iUROSAbmYdZjZkZpuq3y80s4fM7HEz22Bmr0iumSIi0kgzPfRLgF1Tvr8K+Iq7vxZ4HrgozoaJiEhzIgV0M1sAnAN8s/q9AcuAO6qn3Az0J9FAERGJJmoP/Vrgs8CB6vfHA6Puvr/6/R6gfYeWRURyoGFAN7MVwLPuvm3q4YBTPeT5F5vZoJkNjoyMtNhMERFpJEoPfSnwDjPbDXybSqrlWqDbzGp17AuAp4Oe7O7r3b3P3fvmzWu4abWIiLSoYUB39zXuvsDde4ELgC3u/n7gfuDvqqddCNyVWCtFRKSh2dShrwL+0cx+QSWnfkM8TRIRkVY0NfXf3R8AHqh+/UvgjPibJCIirdBMURGRklBAFxEpCQV0EZGSUEAXESkJrYcuIpKQgaHhVDfcUEAXEUnAwNAwazbuZGx8AoDh0THWbNwJkFhQV8pFRCQB12x+bDKY14yNT3DN5scSu6YCuohIAp4eHWvqeBwU0EVEEnBCd1dTx+OggC4ikoCVyxfR1dlxyLGuzg5WLl+U2DU1KCoikoDawKeqXERESqB/SU+iAXw6pVxEREpCAV1EpCSUchGRtpD2rM0sKKCLSOllMWszC0q5iEjpZTFrMwvqoYtI6aU5azPL1I566CJSemnN2qyldoZHx3AOpnYGhoZjvU4YBXQRKb20Zm1mndpRykVESi+tWZtZLMg1VcOAbmaHAw8Cr6yef4e7X2FmNwFvBn5bPfWD7r4jqYaKiMxGGrM2T+juYjggeCe5INdUUXro+4Bl7v6imXUC/2lm91YfW+nudyTXPBGRiiLUka9cvuiQ8khIfkGuqRoGdHd34MXqt53VP55ko0REpipKHXkWC3JNZZV43eAksw5gG/CnwNfdfVU15fJGKj34+4DV7r6v3uv09fX54ODgrBstIu1l6botgamMnu4utq5elkGL0mVm29y9r9F5kapc3H3C3RcDC4AzzOwUYA3wOuB04DhgVUhDLjazQTMbHBkZifwGRERqsh5sLIqmyhbdfRR4ADjb3fd6xT7gX4EzQp6z3t373L1v3rx5s26wiLSfLHb/KaKGAd3M5plZd/XrLuCtwM/MbH71mAH9wE+TbKiItK8sdv8poihVLvOBm6t59DnA7e6+ycy2mNk8wIAdwEcTbKeItLGsBxuLItKgaFw0KCoi0rxYB0VFRCT/FNBFREpCAV1EpCQU0EVESkIBXUSkJLR8rohIzLJaSEwBXUQkRlkuJKaUi4hIjLLctUg9dBHJXBHWOo8qy4XE1EMXkUxlvbFy3LJcSEwBXUQylfXGynHLciExpVxEJFNlW+s8y4XEFNBFJFNZb6ychDQ2pA6ilIuIZEprncdHPXQRyZTWOo+PArqIZC6rFEXZKKCLSObKVIeeJQV0EclUllPly0aDoiKSqbLVoWdJAV1EMlW2OvQsKaCLSKaynCpfNg0Dupkdbmb/ZWYPm9mjZvZP1eMLzewhM3vczDaY2SuSb66IlI3q0OMTpYe+D1jm7m8AFgNnm9mZwFXAV9z9tcDzwEXJNVNEyqp/SQ9rzzuVnu4uDOjp7mLteadqQLQFDatc3N2BF6vfdlb/OLAMeF/1+M3AlcD18TdRRMpOdejxiJRDN7MOM9sBPAt8H/i/wKi776+esgcI/N8ws4vNbNDMBkdGRuJos4iIBIgU0N19wt0XAwuAM4CTg04Lee56d+9z97558+a13lIRKaWBoWGWrtvCwtX3sHTdlsKugz7DxAR8/vNgBmeemcolm5pY5O6jZvYAcCbQbWaHVXvpC4CnE2ifiMQsT7MySzmp6BvfgI9//NBjc9IpKIxS5TLPzLqrX3cBbwV2AfcDf1c97ULgrqQaKSLxyNvuQKWZVPTcc5WeuNmhwfz1r6889qMfpdKMKLeN+cD9ZvYI8BPg++6+CVgF/KOZ/QI4HrghuWaKSBzyFkALP6noYx+rBPHjjz/0+Kc/De7w8MNw7LGpNSdKlcsjwJKA47+kkk8XkYLIWwAt5OYWP/85LAqpkR8dhWOOSbc9U2imqEgbyduszMJMKnKH00+v9ManB/NvfrPyuHumwRy02qJIYvI0+FizcvmiQwYhIb4A2sr7zf3mFjfcAP/wDzOPd3bCyy/DYfkKoflqjUhJJFW9MdubRFIBdDbvN3eTil5+GY48MvixBx6AN7851eY0wyoTQdPR19fng4ODqV1PJCtL120JzA33dHexdfWyll5zetCESu/6XX/Zw/0/G8m0h5vE+01dTw88HVJ9nWKcDGJm29y9r9F5yqGLJCCJwcewCpVbf/xU5mWIeRtsjezxxw+WG04P5t/5zsHceEEooIskIInBx7DgOD3cZFGGmLfB1oZqQfzP/mzmY7Ug3t+ffrtmSQFdJAFJVG80ExzT7hkXolrl6qsPBvLpHn98MpAXeSkCDYqKJCCJwcegChUjeBGltHvGua5WCQrgNdPSKUVfikCDoiIFMr3K5W9eN487tw3PGCht+/XEu7vht78NfuyFF+DoowMfyuvgbtRBUfXQpdTyWAs+G0Elfn1/fFyp3mPLRkfDp9mfeio88kjDlyjs4G6VArqUVtE/PkeVuzruiGK72TaRUmmkkEsRTKFBUSmtvC1EJQc1s+pj4CDlli3hA5wXXdRyuWEhBnfrUA9dSqvoH5/LrN7NdmovffqnrK1rzoI1IS8aw3hgrgd3I1BAl9Iq+sfnMot6s71m82Os3biO/v/+QfALfe97sHx5rG0ragoLFNClQJrNuSa5EFUWyjTAe0xXJ6Nj4zOOH3KzNWNryPMXrtrEE+vOSaZxBaaALoXQygBn0T8+T1WmAd6BoWFe+sP+Gcc751jdlMppn7yV546oLE/bo09ZgRTQpRCi5lynK/LH56laff95dM3mxxifOJjvnvvS8wx+7QOh5598+b2l+ZSVNAV0KYR2H+As0/uvtXn3VSvCTzpwYLKCZW2JUk1JU0CXQmj3Ac7SvP9rr+WJqz4d/nhApUqWn7KKNm6hOnQphKLXB89W0d7/9NrxyZrxT88M5r2rNnHy5fcysH1PBi0N10ytfF407KGb2YnALcBrgAPAene/zsyuBD4MjFRPvczdv5tUQ6W9lWmAsxVFev+1QLjrC28LPefnH/wYH1r0Lp4eHaMnp++liOMWDRfnMrP5wHx3325mRwPbgH7gfOBFd//nqBfT4lwixfsY37QYp+JnaeHqewJXsjRIvWQytsW53H0vsLf69e/MbBdQop8+kfTkufxwVjeaOkH8bz78f3jiuJ5KIIynqako4rhFUzl0M+sFlgAPVQ99wsweMbMbzSxkmTMRqcnr+jIt5YuHh8PXU6GSG+9dtYknjqvcFPIcCIMUbdwCmgjoZnYUcCdwqbu/AFwP/AmwmEoP/kshz7vYzAbNbHBkZCToFJG2kdfyw7AbzWduf3hmUK8F8QULZr7QgQMMbN/DyZffe8jhvAfCIP1Lelh73qn0dHdhVCYz5X2d+UgbXJhZJ7AJ2OzuXw54vBfY5O6n1Hsd5dCl3eV1A4WwfDFUgvHdT93Fa2/53+EvELDzT6nHCVIWWw7dzAy4Adg1NZib2fxqfh3gncBPW22sSLtotL5MVoEwLF9cd/JPnc5gWWboFk2UiUVLgQ8AO81sR/XYZcB7zWwxlS0NdwMfSaSFIiVSr/wwywHTqTeaukH8r/8aHnww0bYUSd4+iWhPUZGcyDwdU6dSZena+zJNCeXR9BswJLefq/YUFSmYTAZM6wTxd79vHT858ZRKkCrYgGYa8jjxSAFdJCdSq3veswdOPDH04YHteybTCHmdxZkHeaxYUkAXyYnEN+SoN4NzYgLmVKqY+8l+klMR5HHikRbnEsmJROqeP/7xupN/JjdTntM4FARu1tzG8jjxSD10kRyJrdwv5vVU8rxkQVbyuGCaArpIxmIrfasXxLu74fnnW25jHgcA8yBv9fYK6CIZiqXnm8LqhnkcAExK3mrLm6EcukiGWl6sq5YXDwrm99xzMDcek7CBvqItuNVIETe1mEoBXUqjiIN2TfV8G6xuOBnE//ZvY2xhRR4HAJOQ19Uwo1LKRUoh60G7Vj+mRyp9q5dS2b8fOjrCH49JHgcAk1D01JICupRCloN2s7mZhNWef2vbTWBnhT8xg51/8jYAmIQ81pY3QwFdSiHLntWVdz/a8s1kes/3iRZXN5R4JD65K2EK6FIKWfWsBoaGGR0bD3ws6s2k/7QF9Nc7IUeBvMgVIFEUPbWkgC6lkFXPqt5gWcObScE2U856nCItRU4tqcpFSiGr7cLq9cIDbyb1yg3/7d9iLzeMU9ErQNqBeuhSCFE+6mfRswpL9Rx7ROfBtjzzDLzmNeEvktMAPl3RK0DagXroknt5mOwRVuMeVp99xdv/4mBPPCiY79uX6954kHaZXFRk6qFL7mW9jkiU3HHt08MXf3QL7/vh7fCFkBebZQDPclCy6BUg7UABXTLXKEjF+VG/lYDY6IbSv6SH/tMWhL9ATL3wrAcli14B0g4U0CVTUYJUXCWJrQbEsBvH1jVnwZo6F4w5nZL1JxUodgVIO1AOXTIVpXIirnVEWq3SmH7j2H3VCnaHTQCq5cUTyI3HPShZxLVvpL6GPXQzOxG4BXgNcABY7+7XmdlxwAagF9gNnO/urS+4LG0pSpCK66N+2LWGR8dYum5L6GuuXL6ofkpl/Xr48Iebaksr4pw8lXX6RpIRJeWyH/iMu283s6OBbWb2feCDwH3uvs7MVgOrgVXJNVXKKGqQiuOjfti1ICSg/eY3MHdu+CzOlCtU4hyUzEP6RuLXMOXi7nvdfXv1698Bu4Ae4Fzg5uppN0P92csiQeJclrVRCiHoWlNNpl9q5YZz5wacNJZaueH09wPENnlKNeXl1NSgqJn1AkuAh4BXu/teqAR9M3tV7K2T0ms2nRJWpRIlhTD1WtN76p/9wU187Md3hDc05d542PtZe96pbF29bNavX/RVBSWYecQfVDM7CvgB8EV332hmo+7ePeXx59392IDnXQxcDHDSSSf95ZNPPhlPy6XtTA9yUOnNrz3v1MAgDZVebFAAXLpuC8OjY+GDm5DppJ9a+6YLez/NqvdvqZRL/pjZNnfva3RepB66mXUCdwK3uvvG6uFnzGx+tXc+H3g26Lnuvh5YD9DX11ecaXGSO/XyvvUGPAeGhg8NUmZsrXehHMzeTDoloprycopS5WLADcAud//ylIfuBi4E1lX/viuRFkoptTLBp16QqzfgOZl6qVOpsnTtfbkJaANDw8wxYyLgxhJnSkQ15eUTpYe+FPgAsNPMdlSPXUYlkN9uZhcBTwHvTqaJUjatlszVy/sGVYAAB1MqQVPxr70WLrkEoH6PPUW1f5ugYN5O0+zjWOKg7Gu3B2kY0N39P4GwhZvr7JElEqzVkrl6ZXu15126YQd/9PsXeeS6C8IbkIOUSpigfxuADrO2yW/HUSPfrnX2mikqqWs1P9xozfP+0xaw+6oVgcH8rH/aVIjVDcP+DQ64lzoQTRXHuuvtuna71nKR1LVSMjf94/NX3rO4EuA+9zk4LWxpQ+hdtalSvfH2U2Npe9KKWk4YZ3ojjgHhdq2zVw9dUtfsZKKg9dD7T1tQmfzzhZnBfGD7HpauvY+FqzaltnNRXOKcaJWWuNerj2Pd9XZdu109dElMWK+t2ZK52sfnujXjHFqpktY66a30Sus9r4jlhHEvIxDHEgftunZ75IlFcejr6/PBwcHUrifZiXXiSp3NlE++/N5MJse0+v7KNqFnYGiYSzfsCHzMgCfWndPy66rK5aCoE4sU0CU2U3+BwuqoI890rBPEr37T/+IbbzyfjtleYxZancmZ9AzQNAXdnKYq4nvKq1hnioo0Mv2XOyjQQoNBqZdegqOOCn24d9Wmya+7OjtCA0kaA1+tDrqVabAurMQSoHOO8fIf9rNw9T2F7x0XiQZFJRb1frmnChyUqq1uGBTMX3gB3BnYvmdGuWJPhgNfrQ66lWmwru5NyOD5l8cz29S7XSmgSyzCpt1Pdcig1Pr1BwN5kFrN+NFHA5XBwq2rl/HEunPYunoZ/Ut6Mq0IafXaRaxiCRN2E+owY3zi0E9o7VADngdKuUgswvLZUBkcm/zYHeNmyllWhLR67SJWsYQJqyTJMhXW7jQoWnB5GcnvXX1P6GONyg3zPntTwgX9/DW7lLE0pkHRNpCn9Sp6AmY45nWtcYlPWM1/O9aA54Fy6AWWp/Uqarnh3VetmPwzw9VXF2I9FZmdRmvuSHLUQy+weiVwqaZi9u2j/7QFiW+mnJf0kjSmtdazoYBeYGELOR3T1ZlOKqbO5B9GR+GYY2K7VJ7SSyJ5pZRLgYWVwJmRXCpm48Zo5YYxBnPIV3pJJK8U0AssLFc5+vJ44PmzKhurBfF3vWvmY7UgnmBuvN6eoQtX38PSdVs0cUXanlIuBReUqwwrG2t6NuKCBTBcJ0imOLhZb8/QqbMRQSkYaV/qoZfQrGcj1nrjQcE8hd54kKD3NJ1SMNLu1EMvoZZmI9Yb4Lz+evjoRyNfP4lqlOnvKex2otmI0s4U0EsqUtnY+Di84hXhj7fQC0+qGmX6TeKlffsZHZs5VlDERa5E4tIwoJvZjcAK4Fl3P6V67Ergw8BI9bTL3P27STVSYlavN/7889DdHellwqZ9B1WjfOb2h4HWgnrQTaKzw+icY4wfOHjT0WxEaXdRcug3AWcHHP+Kuy+u/lEwz7vNm6OVGzYRzIP2kQwbuJxwb3kJ1aCbxPiEc9Thh2k2osgUDXvo7v6gmfUm3xSpp+W8dL3e+CwGNsN64vVWXWx1n8mwvPjoy+MMff5/NvVaImU2myqXT5jZI2Z2o5kdG3aSmV1sZoNmNjgyMhJ2mtTR9K7qb3pTtN74LIQF2Qn3utUoUdZNn65Mm0KIJKnVgH498CfAYmAv8KWwE919vbv3uXvfvHnzWrxce4s8S7IWxH/4w5kvEnO5YVgwraU+OkJuJgZNp13KtCmESJJaCuju/oy7T7j7AeBfgDPibVY+DQwNs3TdltRnJtbdh7IWxAMC6Jff/kkGtu9JpGa8XpDtX9LDl85/A0Eh3aHpWnGt3icSTUtli2Y23933Vr99J/DT+JqUT1kuDjV9luScAxP88ppzQ88/ZDPlhNrYqNa9f0kPl27YEfjcVmrFtXqfSGNRyhZvA94CzDWzPcAVwFvMbDGVDtdu4CMJtjEX6qU9kg40ta2+dn3hbeEn/eY3LF2/Y0aOOsk2NgqyQZtegHLfIkmJUuXy3oDDNyTQllyrm/ZI0uAg/aefHmmt8czaGCJsz0nlvkWSobVcIgrrVXYf0ZnMBWt58dNPn/HQZF58Wm48b9Ugyn2LpEtT/yNauXwRK+94mPGJQ4Poi7/fz8DQcDxB6iMfgfXrAx8aPfwoFl/ybSA8L57HHrFy3/miXZ/KTQE9ov4lPVy28ZEZAX38gM8+R11n8s/StfdFzou3tCiXtA3t+lR+CugRDQwN8/L4gcDHWspR15vB+a1vwfvfX3nt1fcEnjJc3Tc0KKjrl1OCZDmwL+lQQI+oXu105Bz1gQPQUWdN74B68XobO6h3Jc3I26C5xE+DohHV+6FvmKOuDXAGBfPnnqs7g7Pexg5ZbeiQ1QQrmZ28DZpL/BTQIwqtcunqDO4h79wZbT2VY0OXwQEOVoqESbt31fS6MpIbWkKh/BTQIwr7ZVjxhvmH9FYng/jrXz/zRVpcT6V/SQ89OeldRV5XRnJHZaTlpxx6RP1Lehh88jlue+hXTLjTYcZpJx3DnduGWfXdr/PB7ZuCn3j88fDrX8/6+nkpSVQettg0aF5ubRvQm63HHRga5s5tw5NrfU+4c+vFfxV6fu+qTXR1dlR6QC1cb7q8lCSGDdIqDyuSvbYM6GH1uINPPsf9PxsJDJi1VMMtGz7Hm3YPBb7uZ8/+FLe/4eCGC1NTEXHU/+ahd5WXTwoiMpN5Akurhunr6/PBwcHUrhdm6botgb1Mg0N2k5/sYS8+AeaEDzdMXd0w6DXDerU93V1sXb2siZbng2YbiqTLzLa5e1/D89oxoC9cfQ9R3vXuq1aEPnbqpRv43SuPbPgaPd1dPF2tCAl7XAFRROqJGtDbMuVSb7JO99gL7Pjq+wIfG+06msWfui3ydWqpiGs2PxZ6PU2/jk6fDETqK0VAb/YXPSgPXK83XiszfGBomI7bHw7dBBkOpm2m97ynX28qTb9uTOuQiDRW+IAe9It+6YYdXHn3o1z5jr8I/GWvHfuP6zfwtX/5TODrfvWN7+H6sz44WaUy9XnTg3NYEJ9+vXo9dZX91ad1SEQaK3xAD/pFBxgdGw/vwZnRD4GbRixde99kT39tg+DczEf/WoVK2ICsyv7qU/27SGOFD+j1fqEP6cGtWgVXXx184vbtsGQJAFsjXHM25YMq+2uN6t9FGitsQK/lzetWq7izdc1ZsCb88Vavm8QEIQ36hdONUKSxQgb0ywd2cuuPnwoN5juuu4Du378Y/OCLL8KRjcsNg8Q1MBfUw9egX315mSkrkmcNA7qZ3QisAJ5191Oqx44DNgC9wG7gfHd/PrlmVgwMDXPl3Y8yOjY+47Gj973EzmvfE/i8J151Ess+9I1KEPj5KP1LWgvoSQ7MadCvsTzMlBXJsyg99JuArwG3TDm2GrjP3deZ2erq96vib97BFMfw6NiMmZwA9974CU4e2R343Mu/8wh3bhuOrdeb5MCcBv1EZLYaBnR3f9DMeqcdPhd4S/Xrm4EHSCCgT09D1IL5MWO/4+GvvjfwOWuWf4KB089h7Xmncn/Mvd4kB+byOuinvL5IcbS6Hvqr3X0vQPXvV8XXpIOmpyGW/eK/2H3VisBg3rtqE72rNvHgW945ucZz3L3eJDcIyOPmA9rMQqRYEh8UNbOLgYsBTjrppKaeWwu8YT3y11/ybV44/CgM+PszT+IL/Yfu7BN3rzfJgbk8Dvopry9SLK0G9GfMbL677zWz+cCzYSe6+3pgPVQW52rmIrWAfPj+fQD89pVH8r4Lvsijr/nTyXOOPaKTK94ePCM0iVK3JAfm8jbop7y+SLG0GtDvBi4E1lX/viu2Fk1RC8jPHD13cona2g6dUVYpzGOvt0jymtcXkWBRyhZvozIAOtfM9gBXUAnkt5vZRcBTwLuTaFwcATlvvd4i0WQekWJpy/XQJTpVuYhkT+uhSyz0CUekOBTQRaQtlfHTpwK6iLSdsq6d1OrEIhGRwqo3x6LIFNBFpO2UdY6FArqItJ2wuRRFn2OhgC4ibSePayfFQYOiItJ2yjqLXAFdRNpSGedYKOUiIlISCugiIiWhgC4iUhIK6CIiJaGALiJSEqkun2tmI8CTLT59LvDrGJtTFHrf7aMd3zPofUfxx+4+r9FJqQb02TCzwSjrAZeN3nf7aMf3DHrfcb6mUi4iIiWhgC4iUhJFCujrs25ARvS+20c7vmfQ+45NYXLoIiJSX5F66CIiUkfuA7qZnW1mj5nZL8xsddbtSYOZnWhm95vZLjN71MwuybpNaTKzDjMbMrNNWbclLWbWbWZ3mNnPqv/vb8y6TWkws09Xf8Z/ama3mdnhWbcpCWZ2o5k9a2Y/nXLsODP7vpk9Xv372NleJ9cB3cw6gK8DbwP+HHivmf15tq1KxX7gM+5+MnAm8PE2ed81lwC7sm5Eyq4DvufurwPeQBu8fzPrAT4F9Ln7KUAHcEG2rUrMTcDZ046tBu5z99cC91W/n5VcB3TgDOAX7v5Ld/8D8G3g3IzblDh33+vu26tf/47KL3e51vkMYWYLgHOAb2bdlrSY2R8BbwJuAHD3P7j7aLatSs1hQJeZHQYcATydcXsS4e4PAs9NO3wucHP165uB/tleJ+8BvQf41ZTv99Amga3GzHqBJcBD2bYkNdcCnwUOZN2QFP0PYAT412qq6ZtmdmTWjUqauw8D/ww8BewFfuvu/55tq1L1anffC5VOHPCq2b5g3gO6BRxrm7IcMzsKuBO41N1fyLo9STOzFcCz7r4t67ak7DDgNOB6d18CvEQMH7/zrpozPhdYCJwAHGlmf59tq4ot7wF9D3DilO8XUNKPZNOZWSeVYH6ru2/Muj0pWQq8w8x2U0mvLTOzb2XbpFTsAfa4e+1T2B1UAnzZvRV4wt1H3H0c2Aj8VcZtStMzZjYfoPr3s7N9wbwH9J8ArzWzhWb2CioDJndn3KbEmZlRyafucvcvZ92etLj7Gndf4O69VP6vt7h76Xts7v7/gF+ZWW2H4rOA/86wSWl5CjjTzI6o/syfRRsMBk9xN3Bh9esLgbtm+4K53lPU3feb2SeAzVRGwG9090czblYalgIfAHaa2Y7qscvc/bsZtkmS9Ung1mrH5ZfAhzJp2XBJAAAAZUlEQVRuT+Lc/SEzuwPYTqWya4iSzho1s9uAtwBzzWwPcAWwDrjdzC6icnN796yvo5miIiLlkPeUi4iIRKSALiJSEgroIiIloYAuIlISCugiIiWhgC4iUhIK6CIiJaGALiJSEv8fVa/kGbcv+t8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# normal linear regression\n",
"from sklearn.linear_model import LinearRegression\n",
"lr = LinearRegression()\n",
"lr.fit(X, Y)\n",
"plt.scatter(X, Y)\n",
"plt.plot(X, lr.predict(X), color='red')"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"# 一个无敌简易版的分位数回归\n",
"class QuantileRegression:\n",
" def __init__(self, r=0.5):\n",
" self.r = r\n",
" self.batch = 8\n",
" self.w, self.b = 0.0, 0.0\n",
" self.lr = 0.1\n",
" \n",
" def fit(self, X, y):\n",
" for n in range(5000):\n",
" random_idx = np.random.randint(0, len(X), self.batch)\n",
" train_X, train_y = X[random_idx], Y[random_idx]\n",
" pred_y = self.w * train_X + self.b\n",
" grad_w = np.mean((pred_y > train_y) * (1 - self.r) * train_X) + np.mean((pred_y < train_y) * self.r * - train_X)\n",
" grad_b = np.mean((pred_y > train_y) * (1 - self.r)) + np.mean((pred_y < train_y) * self.r * -1)\n",
" self.w -= self.lr * grad_w\n",
" self.b -= self.lr * grad_b \n",
" print('training done! w: %.4f b: %.4f' % (self.w, self.b))\n",
" \n",
" def predict(self, X):\n",
" return self.w * X + self.b\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"training done! w: 0.9305 b: 9.5375\n",
"training done! w: 2.1781 b: 10.0625\n",
"training done! w: 2.8371 b: 10.1125\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fe23efff940>"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNXd+PHPyQIJa9gUSAjBB0VF9khFlIJacUEFd2mt1io+v2qrVq2gFVyqULVi+9TSovXBWp+iAiLigq1IsYogGCIiAkpQEsCwBSEJSSY5vz9mJkwydzJL7jrzfb9evkzO3Jl7JiHfOfd7vudcpbVGCCGE96U53QEhhBDmkIAuhBBJQgK6EEIkCQnoQgiRJCSgCyFEkpCALoQQSUICuhBCJAkJ6EIIkSQkoAshRJLIsPNk3bt31wUFBXaeUgghPG/dunV7tdY9oh1na0AvKChg7dq1dp5SCCE8Tyn1dSzHScpFCCGShAR0IYRIEhLQhRAiSdiaQzdSV1dHaWkpR44ccborKSErK4u8vDwyMzOd7ooQwmSOB/TS0lI6duxIQUEBSimnu5PUtNbs27eP0tJS+vXr53R3hBAmczygHzlyRIK5TZRSdOvWjT179jjdFSEssbiojMeXbWZnRTW9c7K5e/wAJg7LdbpbtnE8oAMSzG0kP2uRrBYXlTFt0Qaq6+oBKKuoZtqiDQApE9RlUlQIkRQeX7a5MZgHVdfV8/iyzQ71yH4S0G3w1FNPUVVV1fj9BRdcQEVFBQAdOnSI+vy3336bAQMG0L9/f2bNmmV4zMqVKxk+fDgZGRksWLDAnI4L4SE7K6rjak9GEtBt0Dygv/nmm+Tk5MT03Pr6em655RbeeustPv/8c/7xj3/w+eefhx2Xn5/PvHnzmDx5smn9FsJLeudkx9WejCSgA4888ggDBgzgnHPO4ZprruGJJ55g7NixjdsU7N27l+AeNNu3b+fMM89k+PDhDB8+nA8//BCAFStWMHbsWC6//HJOPPFEfvjDH6K15g9/+AM7d+5k3LhxjBs3DvBvgbB3796wfjz++OOceuqpDB48mBkzZgCwZs0a+vfvz3HHHUebNm24+uqree2118KeW1BQwODBg0lLk1+pSE13jx9AdmZ6k7bszHTuHj/AoR7ZzxWTokFLP93Jrgpz69F75WQxYXDviI+vW7eO+fPnU1RUhM/nY/jw4YwYMSLi8ccccwz//Oc/ycrKYuvWrVxzzTWNgb+oqIiNGzfSu3dvRo8ezQcffMAvfvELnnzySd577z26d+8e8XXfeecdtm7dypo1a9Bac/HFF7Ny5UrKy8vp06dP43F5eXmsXr06gZ+EEMktOPEpVS4p7P3332fSpEm0a9cOgIsvvrjF4+vq6rj11ltZv3496enpbNmypfGxkSNHkpeXB8DQoUPZvn07Z5xxRkz9eOedd3jnnXcYNmwYAIcPH2br1q106tQp7FipVBHC2MRhuSkVwJtzVUBvaSRtJaMAmZGRQUNDA0CTVayzZ8/m2GOPpbi4mIaGBrKyshofa9u2bePX6enp+Hy+mPugtWbatGncfPPNTdpXrVrFjh07Gr8vLS2ld29nfk5CCHdL+YTrmDFjePXVV6murubQoUO8/vrrgD8nvW7dOoAmVSMHDx6kV69epKWl8cILL1BfX2/4uqE6duzIoUOHWjxm/PjxPPfccxw+fBiAsrIyysvLOfXUU9m6dSslJSXU1tYyf/78qFcRQojUFHNAV0qlK6WKlFJLA9/3U0qtVkptVUq9pJRqY103rTN8+HCuuuoqhg4dymWXXcaZZ54JwF133cWcOXM4/fTTm0xg/uxnP+P555/ntNNOY8uWLbRv3z7qOaZMmcL555/fOClq5Nxzz2Xy5MmMGjWKQYMGcfnll3Po0CEyMjL44x//yPjx4znppJO48sorGThwIADTp09nyZIlAHz88cfk5eXxyiuvcPPNNzceI4RIHUprHduBSv0SKAQ6aa0nKKVeBhZprecrpf4MFGut57T0GoWFhbr5DS42bdrESSedlFjvLfDAAw/QoUMH7rrrLqe7Yhm3/cyFEC1TSq3TWhdGOy6mEbpSKg+4EHg28L0CzgKCuYjngYmJdVUIIYQZYp0UfQr4FdAx8H03oEJrHZz1KwWSYmr5gQcecLoLQgiRkKgjdKXUBKBca70utNngUMPcjVJqilJqrVJqrezyJ4QQ1okl5TIauFgptR2Yjz/V8hSQo5QKjvDzgJ1GT9Zaz9VaF2qtC3v0iHrTaiGEEAmKGtC11tO01nla6wLgamC51vqHwHvA5YHDrgPC16MLIYSwTWvq0O8BfqmU+hJ/Tv2v5nRJCCFEIuIK6FrrFVrrCYGvt2mtR2qt+2utr9Ba11jTRe+zY/vcefPm0aNHD4YOHcrQoUN59tlnzem8EMIzUn6lqB3s2D4X4KqrrmL9+vWsX7+eG2+80ZS+CyG8QwI6ybF9rhBCuGpzLj5bBN+VmfuanXLhlEsjPpxM2+cuXLiQlStXcsIJJzB79uwmzxNCJL+UH6GHbp/bqVOnmLbPvemmmxg0aBBXXHFFk/RHcPvctLS0xu1zYxW6fe7w4cP54osv2Lp1K0ZbMxjtDnnRRRexfft2Pv30U8455xyuu+66mM8thEgO7hqhtzCStlIybJ/brVu3xq9vuukm7rnnnpjPLYSwxuKiMltvuJHyI/Rk2T53165djV8vWbJENt8SwmGLi8qYtmgDZRXVaKCsopppizawuMjktHIId43QHRC6fW7fvn2bbJ975ZVX8sILL3DWWWc1Hv+zn/2Myy67jFdeeYVx48bFtX1ur169eO+99wyPOffcc9m0aROjRo0C/OWMf//73znmmGMat8+tr6/nhhtuaLJ9bmFhIRdffDF/+MMfWLJkCRkZGXTt2pV58+a18icjhGiNx5dtpiZ9O227fQwNbagpv5DqunoeX7bZslF6zNvnmkG2z3UHt/3MhUg2H+36iOtfeqZxg6uGml7UHfAP1hRQMuvCuF4v1u1zU36ELoQQZvn6u6+Zt3EeAO3aZFBZ66N2z7no+qMLCHvnZFt2fgnozcj2uUKIeH1X+x2z181u0nbVwPHM+1c7dMg8W3ZmOnePH2BZPySgCyFEgnwNPh5Z/UiTtr4d+3L9KdcDMKCTvVUuEtCFECIBD656MKxtxqgZTb6fOCzX0gDenAR0IYSIw4ubXuTLii+btN37vXvJTMt0qEdHSUAXQogYfFD2Af/65l9N2m4ffjud23Z2qEfhUn5hkR3s2D7366+/5uyzz2bw4MGMHTuW0tJSczovRJJYXFTG6FnL6Tf1DUbPWh7zAp+i8iIeXPVgk2B+7cnXMmPUDFcFc5CAbgs7ts+96667+PGPf8ynn37K9OnTmTZtmmn9F8LrElm1WV5VzoOrHmTJV0sa287JP4cZo2ZwXOfjbOh1/CSgkxzb537++eecffbZAIwbN0622BUixOPLNlNd13SbjuCqzeaqfdU8uOpB5hTPadI+Y9QMRueOtrSfreWqHPrb29/m28pvTX3NY9sfy3kF50V8PFm2zx0yZAgLFy7ktttu49VXX+XQoUPs27evyaZdQqSqnRXVMbXHUrkSjd0bcoVyVUB3Quj2uUBM2+feeuutrF+/nvT0dLZs2dL4WHD7XKBx+9wzzjgjpn6Ebp8LcPjwYbZu3UqnTp3CjjXaHfKJJ57g1ltvZd68eYwZM4bc3FwyMlL+1ysE4F+dWWYQ1IOrNo0C+dSRU2mb3jasvSXB1E7waiCY2gFsCequ+otvaSRtpWTYPrd3794sWrQI8H8YLFy4kM6d3TVhI4RT7h4/oEmgBf+qzRMGvs2Dq95vcuxPBv6E/E75CZ2npdSOHQE95XPoybJ97t69exs/gGbOnMkNN9wQ2w9AiBQwcVguMy8dRG5ONgronruK00Yup6D70d1ST+99OjNGzUg4mEPsqR2rRB2hK6WygJVA28DxC7TWM5RS84DvAwcDh16vtV5vVUetkizb565YsYJp06ahlGLMmDE8/fTTrf3RCJFUJg7LpX+fAyzcuhDIDPznF2+ePJJoqR2rRd0+V/nzEe211oeVUpnAf4DbgP8GlmqtF7T4AiFk+1x3cNvPXIhYtGay8WDNQZ765KmwdrMCeWgfjVI7My8d1KqUi2nb52p/xD8c+Db4sWbfJupCiJSX6GSj1pqHPnoorN3sQB4U7Iurq1yUUunAOqA/8LTWerVS6v8BjyilpgPvAlO11jXWddUesn2uEO6TyGSjUeXKPafeQ1ZGlsHR5rF7Q65QMQV0rXU9MFQplQO8qpQ6BZgG7AbaAHOBe4Cwj0Kl1BRgCkB+vvFkg9basNJEmM/OO1QJYZZ4JhuNAvnVA65mQFfr9iF3i7jKFrXWFUqpFcB5WusnAs01Sqn/BQyTzlrrufgDPoWFhWHRJCsrq3EBjAR1a2mt2bdvX5NSSyG8IJbJRqNA3j27O7cMvcXSvrlJLFUuPYC6QDDPBs4BfquU6qW13hWYNJ0IfJZIB/Ly8igtLWXPnj2JPF3EKSsrq3HxkxBeEamO/O7xA1i6bSnrvl0X9hyr8uRuFssIvRfwfCCPnga8rLVeqpRaHgj2CliPv+olbpmZmfTr1y+RpwohUoTRZOPNZ3Wn+MizcKTpsakYyIOili2ayahsUQgh4mF35YobmFa2KIQQbuFU5YpXSEAXQrieUSC/7PjLOKX7KQ70xr0koAshXMsokHds05FfjvilA71xPwnoQgjXeebTZ9hZuTOsPZnz5GaQgC6EcI3tB7fz/OfPh7VLII+NBHQhhOMadAMPf/RwWLtXA7lTdy2SgC6EcJRZdwtyCyfvWiQBXQjhCKNAfkG/Czi156kO9MY8Tt61SAK6EMJWRoF8+95Ktmw8j3mvldM7Z7mtW86azcm7FklAF0LYwiiQAwzJupH5azZQXecPeHbfWNlsTt61SAK6EMJSH+36iGXbl4W1Byc8R89a7uiNlc3W0kZiVpOALoSwRF1DHY+ufjSsvXnlitM3Vjabk3ctkoAuhDCdUXpl2shptElvE9bu9I2VreDUXYskoAshTGMUyMf1GceYvDERn+NkiiLZSEAXQrRapAnPWBYGOX1jZcvVVoJugLYdLT+VBHQhRMJaE8hDOXljZcuUfQKfhGxjcNHvLT+lBHQhRNxe/+p1Pin/JKw90aX6Ti2VN52vFt66O7x9yGRbTi8BXQgRs2pfNY99/FhYe2v2XHFyqbxpdn8GHz8T3j76duhq3y02JaALIWJilF6593v3kpmW2arXdXKpfKutmAWHdoW3j38U2rS3vTsS0IUQLTIK5GNyxzAuf5wpr++5OvSq/fCuwdxBryFQeIP9/QkhAV0IYcisCc9oPFOHvvkt2PJ2ePsZv4Qufe3vj4GoAV0plQWsBNoGjl+gtZ6hlOoHzAe6Ap8A12qta63srBDCenYF8iBX16FHmuQEuHA2pKXZ258oYhmh1wBnaa0PK6Uygf8opd4CfgnM1lrPV0r9GfgpMMfCvgohLDRn/RzKq8vD2q2+yYQr69A/fRm+/iC8feAkOG6s3b2JWdSArrXWwOHAt5mB/zRwFhCsxXkeeAAJ6EJ4zsGagzz1yVNh7XbeLcg1deiv32bcftb90L67vX1JQEw5dKVUOrAO6A88DXwFVGitfYFDSgHD34ZSagowBSA/P7+1/RVCmMgovfLr7/2a9LR0B3rjkP3b4IMIi35sWAxkppgCuta6HhiqlMoBXgVOMjoswnPnAnMBCgsLDY8RQtjL6sqVeDi2qCjSaLzf9+GUS60/vwXiqnLRWlcopVYApwE5SqmMwCg9D9hpQf+EECZ6cNWDbN9bSfGOg1TV+mjXJoMhfTrzvxeFLxayg+2Lilqa5Dz/Mcjw5n1Mg2KpcukB1AWCeTZwDvBb4D3gcvyVLtcBr1nZUSFE4oIj8u17K1lTsh9fg/9ief83F/PRrnQW55U5ksO2bVHRqqdh7xbjxzyWVmlJLCP0XsDzgTx6GvCy1nqpUupzYL5S6jdAEfBXC/sphEhA89RK8Y6D+Bo0NbuPphScXJVp+aKiSGmVwp9Cr8HmnMNFYqly+RQYZtC+DRhpRaeEEK1TXlXOnOLworMD31xsONnl1KpMSxYVffs5rPmL8WNJNBo3IitFhUgyRhOe00+bjlKKd/693FWrMk1dVBRpNA5JH8iDJKALYRG7qzeMAvk5+ecwOnd04/dWrspM5P22elFRvQ/evNP4sbNnQLuu8bwFz5OALoQFrKreMAqaxUeeNTzWaGGQVasyW/N+E1pUtOK3cChCYV2KjMaNKP9CUHsUFhbqtWvX2nY+IZwyepZxaiM3J5sPpp6V0Gs2D5ptey4iI03Rr0d7dh440liGOH3UdNsnOK14v4YipVXyT4chV5l3HpdRSq3TWhdGO05G6EJYwIrqjWCJX9ueixrbfA2aL789jAZqdl/KEWDaLvtvDmFptUr5Jlj9Z+PHJjwFSrX+HElCAroQFrCiemNX5Q7a9lwZ1n5kd9NVjU6UIVpSrSKTnHGTgC6EBcyefHxw1YPk9NpJZcgG1TW7JwHGo1O7yxBNe78NDfDGHcaP2XQ7Ny/f31QCuhAWMGvyMbRyZUifzqwp2c+RA0Oorz4O8Idzo1kwu8sQW/1+XTIa9/r9TWVSVAgXinSTiSFZNzYJmuNO7MHCdWVhI+OZlw7yRACKGMizOsMPHrK3L9g4uRsnmRQVAu9dPsdyt6Dm/S/s29VT75FdxbD2OePHLnwSHNy613P3N21GArpIWl66fG7Nbd9cc3OIaJqNxkv2VlJcepCqGh9Pd/i5/4PI4X3YPXN/0wgkoIukZdtOfq3wxf4veGnzS2Htdt4tyFJaw9Lbw5pL9lby2LZ+vFU31N/QwoetnVdZrr6/aQwkoIuk5fbLZ6NRedIE8iiTnD+atZyyuqa/B6MPW7uvslx5f9M4SEAXScutl89GgfzyEy5nYLeBDvTGZDFWq8T6YevEVZZnUlgGJKALz4j30tttl8+tyZODiyd4930FH/7B+LHxj0Kb9mHNnbMzqaiuC2tv/mHr9qsst5GALjwhkUtvt1w+tzaQg0sneBOsHV9cVEZlrS+sPTNNhX3YuvUqy60koAtPSPTS28nLZzMCeZCrJngjBfJux8Ppt0Z9+uPLNlNXH77+pUNWRth7cdtVlttJQBee4KVL76LyIpZ8tSSsvTUTno6/fxNXckbqc0VVeArGLVdZXiEBXXiCFy69tdY89FH46kYzKlcce/8WLMmP9704eZXl2nmLCCSgC09w+6W3UXrl2pOu5bic40x5fVvf/3e74N+zjB/7/j3QqXfUl2gpELr9dxnkynmLKKIGdKVUH+BvQE+gAZirtf69UuoB4CZgT+DQe7XWb1rVUZHa3HrpbWaevCW2vH+TRuPRAqFbf5fNuWreIkZRN+dSSvUCemmtP1FKdQTWAROBK4HDWusnYj2ZbM4lkkVrArnrLuNNTqu4dYOrePWb+obhTpYKKJl1oa19MW1zLq31LmBX4OtDSqlNgDs/noSwmBm15K64jH/zV1Bf06QpuLfKHZXXH/2gSeClHZ/ANYkX5m2aiyuHrpQqAIYBq4HRwK1KqR8Da4E7tdYHzO6gEG7w7x3/ZkXpirD2eFMrjl/GRxiNl+yt5IJtl5vyQePFQGjEK7n+UDEHdKVUB2AhcLvW+jul1BzgYfz76z8M/A64weB5U4ApAPn5+Wb0WQjbmF254sjo9chB+Od048dGXA+9h/GjWcupNthb5c6Xi4H4groXA6ERr+T6Q8UU0JVSmfiD+Yta60UAWutvQx5/Blhq9Fyt9VxgLvhz6K3tsBB2MUqv3DToJnp3iF7lEYmto9c4cuORPlDqtY57pO7FQBiJ1/Z1iaXKRQF/BTZprZ8Mae8VyK8DTAI+s6aLQtjLysqVaKNXUyZME5jkjPRBA4mlhLwWCJNFLCP00cC1wAal1PpA273ANUqpofhTLtuBmy3poRA2saMEsaXRa6smTD/8I+zbavzYhKdAGd9MOsjogyaU1yY07eK2iqVYqlz+g/GtxaXmXCQFu2rJgyKNXhOaMDWp5DD4+ne+XEy9QSmz1yY07eCaiqUQslJUpKw3tr3B2m/D10U4dZOJmCdM6+vgzbuMX+T4c+HExGqkg0EoGSY07eB4xZIBCegi5dQ31POb1b8Ja3f6bkFRJ0wt2FeluWSa0LSaG+vtJaCLlGKUXrll6C10z+7uQG+aijRh+vdj/w9ef834SSYF8lAyoRkbN9bbS0AXKcHuPHkiQkfHI757l7OytzAkrzP9uje7488Fv4N0+/903TYB6DQ31ttLQBdJzQuBPNTE0seY2Hhr0Wb17haMxmPlxglAp7kxPSUBXSQlLwXyxZ+Uwht3UFXjo13bjKaj8i794Izbne0g7pwAdAO3packoIuk8uKmF/my4suwdjcGcl6/jZK9ldSU7Ke+wV8qWFXjY03JfoqHPuiqQOHGCUCreDm1JAFdJIXa+lpmrpkZ1u7WQB5UXHqwMZgD3Ou7EYBcl4183TgBaAWvp5YkoAvPC6ZXtu+tpHjHQapqfXStmcSvzh3mcM9ClLwPny0Ia66q8fGobzKHadek3W0jXzdOAFrB66klCejCs0Lz5Nv3VrKmZD8+XwY15ZeyE2wdWUW8TI9SO/70xuUc9sDI140TgFbwempJArrwHKMJz+IdB6ncOalJm10jK6PL9JpXf07J+q7hJYfQpFrFSyNft00AWsHrqSUJ6MIzWqpcmffaG4aP2TGyemDJRqrr6nk049nGtvoGf368MaBHKDlMlZGvV3jpA9aIBHTherFUrjg1slpcVMav6uYY/iVV1fhiqh330sjXyxUgsfD6B6wEdOFa1b5qHvv4sbB2o8oV20dWe7+EVf8D63eGPfQ/vknsohu5OdlMtubsjvB6BUisvPQB25wEdOFKRumVe069h6yMLMPjbRtZNZvkrKrxNX4dLDkM8spleqy8XgGSCiSgC1cxCuQndDmB7OoxnP3Ehy0Ga0tHVhGqVdq1zeD2yuvD2ru0y0y6IOf1CpBUIAFduEJLE56OXeqHBPGSvZUUlx5sXJ7PhKf8584rI9sg1TPjooFGr+hpXq8ASQUS0IWjYtlzxfZL/Waj8ZJAjXt9g/anVXyQ3ewDxa5JNCcnJb1eAZIKJKALR8wpnkN5VTnQdIVntyM/DAsQZl7qRwyIlftg+UOGz/nNrpG8W9s0aIZ+oNg1ieb0pKTXK0BSgQR0YatDtYd4ct2Tjd8HV3gGFwWVER6kzLrUT3QB0PIPnKtxD+WGSUkvV4CkAgnowjZG6ZXNG39AZUVtk7bmQcqsS/3QgBhcBBS2AAjCasfdkjs2e1Iy2WvKU1HUgK6U6gP8DegJNABztda/V0p1BV4CCoDtwJVa6wPWdVV4lVEgH3HsCCYcNyGmFZ5mXepfd/hZOmSEB7+qGh+jN14S8TXdkjs284PF6fSNsEYsI3QfcKfW+hOlVEdgnVLqn8D1wLta61lKqanAVOAe67oqvCaWCc9Yg1SrLvUDk5zHtK2jqqbpQ4214y0ENLfkjs38YHFD+kaYL2pA11rvAnYFvj6klNoE5AKXAGMDhz0PrEACuiC+uwWZGaRCUwgFndN5rueCJqmUIXmdWVOyn+W+wSxrODXs+S0FNCdyx0YpkZmXDjLlg0VqypNTXDl0pVQBMAxYDRwbCPZorXcppY4xvXfCU+asn0N5dXlYe0s3mYh39Bsp7xtMIdyv/+L/V10Ja0oUQGNQ7/eTZykuKuOzZZvB5QEtUkpk5qWD+GDqWa1+fbfMCwhzKa119KMApVQH4N/AI1rrRUqpCq11TsjjB7TWXQyeNwWYApCfnz/i66+/NqfnwjUqjlTw+6LwTajMvltQ8yAH/tH8zEsHwdLbmyzDD2rXNoOJ978S1j561nLDgJabk21KwGwtq/vX0s9SUi7uo5Rap7UujHZcTCN0pVQmsBB4UWu9KND8rVKqV2B03gsIH5oBWuu5wFyAwsLC2D49hGcYpVfuP+1+0lSa6edqnvf9Yfq/GKi3w9KMsGD+a98NNJDmnwEqKgsLUm6Z6IzE6pSIW+YFhLliqXJRwF+BTVrrJ0MeWgJcB8wK/P81S3ooXMkokJ/b91xG9R4V0/MTKZkLBrPQfceBxuX4VTW+sA2ywPjORW4OaIuLykhTinqDq2czUyJSU558YhmhjwauBTYopdYH2u7FH8hfVkr9FPgGuMKaLgo3iWfCM5KESua0Znb7eYZplQNZebS7cCoPL9oA1Ic9Hmmy040BLfizMQrmbrqCsJoZNfKpWGcfS5XLfwAV4eGzze2OcCszAnlQXCVzIfuqBKtU6hv8we5e343+vO+FR/O+t7+0HiNumeyMxuhnA5CuVMrkt82okU/VOntZKSpaNHvdbL6r/S6svTUTnjHlhw22qw1Wq/zo28nsrKgmt9moa+KwXB5fttnT1RuRfjYNWid1IAplRo18qtbZS0AXhvZU7eFPxX8KazejciVSydz5nbZF3Hf89Z4/Z9a73wQun2H2VUNdvaozUV4tJzQzvWHGhLBb6uy11uzYX03HrAy6tG9j+fkkoIswRumV6adNxz8/3nrNg+6jGc+SnqYY2asr0GyTrIt+H9fls5snO2PhxQ8ks9MbZnyoOfHBuOtgNX//6Gv2V9YZPj7z0kGWnTtIArpoZBTIJ/WfxOAegxN6vUijtsY/8kDteLu2GQzJ63x0VWdWZ/jB0a1sW7p8Dj4e8RwWSnRU2tLzvPiBZHZ6w4wPNSs/GOvqG3hzwy4+2rY/puOP7dSWCYN7t/q8sYh5YZEZCgsL9dq1a207n4iNmROeQZEWrrx53ALjrWohbJfDoH5T3yDSv9LszHRHFsckujAn2Rb0LC4qizgRrYCSWRcm/LpuqHL54Mu9LP10V1zPGTegB2efdCzpaeZc0YLJC4tEcjI7kIf+ATWvo34041nQUFyaER7QIwTyoEiXz+lKOTbxleioNJkm64IfTpG0Jr1hxlVWPK9ReqCKp9/7Kq7Xz83J4uqR+XTv0DaR7llCAnoKemjVQ2iDMa+ZI/J6remrdnPQCGCoAAASm0lEQVRz+tImxzXWkY+7DzrEtv1PpMtno/I+sGfiK9FJN7dM1pkhUoklQGaaoqrWR7+pb7gqbdTQoLlv8WdxP+/7J/Rg/MBjTZtHsooE9BRSeqiUv37217B2MypXjG4eYeTpDj9n8kXx7UUSKa/sZIliopNuXq1iMdLih5CCA1X+yUGnasDnr/mG4tKDcT/vlz84gR4d3TPqjocE9BRhlF4xc/OssorqFgN54yKgBCelIl0+O1URkuikmxerWCJpKRVWV9/0CtDKtNK+wzU88c6WuJ93cu9OXHtaX9P74yQJ6EnOKJBPPnEyx3c53pwTBOrGZ2buoPn8enBfFQVhi4DM4GRFSKLn9mIVSyROpMJaytm3ZPqEk8luk97q87udVLl4XKSZfCsqV5potgDo/1Z/0/h18w2ytidY6SDcz+jfX6RUWDxb//57yx7e/mx33P0Z0bcLl4/Ii/t5bidVLinAaEHHfe/P4LXSrhQ0qyQxJZBX7Yd3jT8oFre7nDXf5YS153owNyxi19pUmK++gftf25jQuR+ddIrrJyntJgHdw0InItv29G9T72uA4h0HGwO6KYE8wnJ8oLHkcHJeGRuSJDcsWidSWml1yX5Wl8S2GCfU5JH5DMrrbHY3k5IEdA/bWVGNytxLm24rm7RX1foYknUjjy/bzLzXWlE2FkMgD7IjN5yK26F60bY9h1ldsp8xJ/RobIs1kNuxPD6ZSQ7dw4b/6WYqa5vuD16z+1JysjOp8TUkthpx9V+g/HPjx6IsALJSsq2wTBaJTlLecc7xHNMpy+TeJC/JoSex4ITnkD7+/cF9DZrafWPRdV3JzkxHKeJfjRjHaNwJybTC0oseX/ZFxE2nopFRt30koHtI88qVgu7tyUxry8YN49hZd3R/8DtivcmDrwbe+pXxyYZdC3lRBwS2iVTyVlZR7brViF4mk5TeJgHdA6KWIDarCoy6gtLlo3EjkRaxAGhS5440Zko0XXJyr45cO6rA3M4IU0hAd7FE91yJtODj78f+H7we4V7eLg3kQUbvqTlJwRgr3lHB/I93JPRcSZd4iwR0F9p2cBsvfP5CWHusJYihFSenfLeS87M3Nt1vPGjCU2DBJbIV1SjNq2giTeV7cZMrMyU66r7pzH4c16ODyb0RdpOA7iJaax766KGw9kRqySeWPsbEgcHvmm2ub+Fo3Kqb8zb/kKis8VFRHT5J58VNrhJx36sbaEiwQE1G3ckrakBXSj0HTADKtdanBNoeAG4C9gQOu1dr/aZVnUwFRnnyW4feSrfsbrG/SEMDvHGH8WMnnAcDzk+wd8YiLfs2qka58+ViILGgbvQhkZmuyExT1IVEtWRcyFRX38D0BCcpH5l4Cmkm3mRBuF8sI/R5wB+BvzVrn621fsL0HqUYo0A+oMsArj7x6thfxIFJzkgj8Ug57nqtEx6pG31I1NVrurTLpF2bjKRZaJRouqRLu0x+dd6JJvdGeFHUgK61XqmUKrC+K6kl3s2zDPPSpY9FPoHFk5yRRuLpze5U1PzxRCYtI+XFK6rqKJp+blyv5Qbrvj7AgnWlCT1X0iWiJa3Jod+qlPoxsBa4U2t9wOggpdQUYApAfn5+K06XHBLZBTF0NDxIbeOaw8upeVVR0q9r04nO8x+HjDZmd9lQpCBbr3WLW6hGKj1siZdvCpHoqPtHp+UzsLfsXyLik2hAnwM8jL8E+GHgd8ANRgdqrecCc8G/9D/B83ne5v2bmb95flh7LBOejy/bzP36L01+W/UNmuLSg/6A7kDJYaQgG1zcdOfLxYYjdYX/AyqeUboXbgqRaOAGGXUL8yQU0LXW3wa/Vko9Ayxt4fCkkUg5XqsrV16/jVsOfxPWvK7hBBZVjmHiRU1XFdm1gVVLQTZ4vjteWh9WXqgh7rSLm24K0ZqVlL+ZeIqpd4IXormEArpSqpfWelfg20lA/Hdd9ZhEyvGM0it3jLiDTm06tXyy12+HkFDYrm1G482VQ28e0XyvcatKBo1EC7ITh+Vye6xbEMR4PrsDeKKjbqXg0Uky6hb2i6Vs8R/AWKC7UqoUmAGMVUoNxR91tgM3W9hHV4hncyijQD60x1Au6X9JyyeJUK0yJK8zF2y7PGrKwe4NrKIF2VyP5L5lJaVIFrFUuVxj0Bx+6/gkF2lUGdqe0G3fvtsJ//6t8WPn/gbadqQfMDOGVEosfbSTG3PfiY66rxnZh8F54XdkEsJNZKVojCJNAua0y0wskMdROx5rXtxt1SBO5r5lklKkIgnoMbp7/ADuXlBMXf3R3HZaVhlHOq9m+95uTe7hmVAg7zkYTv1pWHM8eXE3joitzn23ZpLy4UsGkpGeZnKP3E3u+pTcJKDHaOKwXO5d9GkgoDfQtudiABr00Xt4RgzkH/4R9m01fixKyWE8eXE3VYNYQUbdrWPnpLlwhgT0GC0uKqOqrqHxZsyhDuw4jxk/mxT+JBOW5Ld0Ywejem4nqkHM9sXu73j+w68Teq4E7sjkrk/JTwJ6jB5a9RBteza9f6fv0MnUV57YtHywugL+FWGkPu4+6HBMXOdt6cYOyTC6SnTUfUVhHsPzu5jcm+TmtklzYT4J6FEEJzyrQm7GrOvbU7tnfOP3d48fYNkGWS3d2MGp0VUiedhH39zEoSO+Fo+JREbd5nDbpLkwnwT0CBZsWcDGfUcn29q1yaCy1kfN7kubHPdE1nNMLM0Lf4F23eHs+1vdj2CgNHORTmtEy8M2NGjuW5zYOrNUnKS0kxsnzYW5JKA3s/PwTp7Z8ExY+/RR0wOBq57RaRu4MG016WmKvt3as3j9TqpqfLRrmwETnrJkVWbU+4TaJFIe9t5XN7C6ZH9Mr9EmXfHgJadY0T3RgmSfNBcS0Bs16AYe/ujhsPbQypXc/0zly/JKtNYopejeoQ3b9lRS36D9S/J9oF5az+0vrW/cpMqsPxYnR1fb9hzmmfdLgMi7JVbVGu+uKOkSd0mGSXMRWcoG9NA8cJf8JQzp07lJLfl937uPjLQMqPfBm3dSsreSbXv8wRz8m249cGA8O3TTSc5glXrzVERr63/tGl1Fm6Rs1ybdMHgf07GtBG8hHJaSAT2YB27o9gptekJlLawp2c+ewzXsKx3N7n3t6f7uHVyQW92433hx6UHqA7c7C90gqyXBSUvAlPpfM0dXz76/ja/2VMb9vCF5OazfURF2pXDvBSeZ0i8hROJSMqA/tOohGro1rbioqcxn484RPJLxrP+ncgTWlPi3Ou3XvT3/qh7A0oZRcZ9rZ0W1o/W/VkxSympDIdwppQL6q1tf5dO9nzYpQQTovvtMbsl4DTKKmrTXN2h+tPsaPvjJ2RRtXA4J3m2npcVBo2ctNy0gJlrTfWynttx+zgkxHy95WCHcKSUC+o7vdvDcxucavw+WIM7YG6jKyHgt7DnBtIo6eARouR48kuCkZaQKFUgs/VJ6oIqn3/sq5n6E8nKeW64MhGhZUgT0SH/ovgYfj6x+pOnBGhZm17Nm9wGah+bZvsvYQ9PVh8GywGDgiHRrtSDlP0VYlUtLHwYtpV8SHXX/9Ix+9D+mQ0LPdSPZh0SI6Dwf0I3+0G9/aT3TP3iAEX27HK1c2V/C/fWdSFMKQiY6q2p8PN3h54w7sQeH15VBC2WBkYJzpCDe/HmxjtTj5eVRd6xkHxIhovN8QG/+hx7cPKvG569c6bZvHXfl9Kdnemd/5A3oN/wH9PvJVQBMDrQV9u0a9ZI+0fLBYN759JnvsjOQxgnVrk161PeayispZR8SIaLzfEAP/kGH7oKYRS2D07bRtyqbwoNd6dkt6+gTLpwNacZBMdbJvliPMxpx9+3Wnm8P1TSWQAKkpymGhNwN55TcTvzwe32jvn4qkX1IhIjOswE9mDdPy95ORudPABiRtoV0GgA4v9y/4KeKQEVLKzbIMjpv6Oj89P7dmP3PCPudN1PQmO6poKq2vkmaJvja89d8w5/e+0om/ULIPiRCROfJgP7rxRt4cc0XtDnmbTKyYGTaF42PBQM5wJ/rL6K+cwGTLzrLlPMuLirjzleKG0fXZRXV3PlKMSMLujZZZWrkv79/HH27RT5GJv1aJvuQCBGd0i1UbAAopZ4DJgDlWutTAm1dgZeAAmA7cKXW+kC0kxUWFuq1a9cm3NnFRWXMWLKBI51foT1HGJi2vfGx0EAeLDnMzkznshG5vPfFnriDwNZvD/HcB9ubtL22vsxw2Xu7NulcMvToayYySTl61nLDlEJuTjYfTDXnA0kI4U1KqXVa68Jox8UyQp8H/BH4W0jbVOBdrfUspdTUwPf3JNLRaIJpiLKKajLaf0GPzqsZnFbW+Ph55T1QKJY0nM5HDSc3tqcrxWUjclm4rqzFUa/WmntfjW0lZaQNqKpr61tdaSKTfkKI1ooa0LXWK5VSBc2aLwHGBr5+HliBBQH9aBqiluN7PUNPdXR71h/s6UGGVtznuwFN00nO7Mx0Zl46yJStXgGuOrUPQ/vksHLLHssm5tw66SeLeYTwjkRz6MdqrXcBaK13KaXiu69ajIIB+cSMLfRQ+0nXirF7u/GC70Kma+OgkpuTzS/O7s/qkv1xb/UKLadLrJyYc+Okn+T1hfAWyydFlVJTgCkA+fn5cT03mG7Y7BtAh12X8qn+L1aQaXhs/2M6cGpBVwDW7zgIRN7qtV2bdKaedyKd2xm/ViRWTsy5cdJPFvMI4S2JBvRvlVK9AqPzXkB5pAO11nOBueCfFI3nJME0hCaNtfpEw2PaZKQxIr+LYZXJkLwc1n69n7r6o6fNzkzn0UmD4g7mQVZuTOW2Ta8kry+EtyQa0JcA1wGzAv8P393KBC1tiNWuTTpD8nIaA3nX9pnc+YMBpKWpJsdJDjhxbs3rCyGMRQ3oSql/4J8A7a6UKgVm4A/kLyulfgp8A1xhReeCgfe3b3/BroNH6JydyU9GF/Df3/8vsjKjL5UPvoYE8MS4Ma8vhIgsah26mVpbhy7sJ1c4QjjPzDp0kcLkCkcI75CALoRIScl49SkBXQiRcpJ1jUVqbq4thEhpLa2x8DIJ6EKIlJOsaywkoAshUk6ktRReX2MhAV0IkXLuHj+A7GZrWZJhjYVMigohUo4b904ygwR0IURKSsY1FpJyEUKIJCEBXQghkoQEdCGESBIS0IUQIklIQBdCiCRh6/a5Sqk9wNcJPr07sNfE7niFvO/UkYrvGeR9x6Kv1rpHtINsDeitoZRaG8t+wMlG3nfqSMX3DPK+zXxNSbkIIUSSkIAuhBBJwksBfa7THXCIvO/UkYrvGeR9m8YzOXQhhBAt89IIXQghRAtcH9CVUucppTYrpb5USk11uj92UEr1UUq9p5TapJTaqJS6zek+2Ukpla6UKlJKLXW6L3ZRSuUopRYopb4I/N5HOd0nOyil7gj8G/9MKfUPpVSW032yglLqOaVUuVLqs5C2rkqpfyqltgb+36W153F1QFdKpQNPA+cDJwPXKKVOdrZXtvABd2qtTwJOA25JkfcddBuwyelO2Oz3wNta6xOBIaTA+1dK5QK/AAq11qcA6cDVzvbKMvOA85q1TQXe1VofD7wb+L5VXB3QgZHAl1rrbVrrWmA+cInDfbKc1nqX1vqTwNeH8P9xJ9c+nxEopfKAC4Fnne6LXZRSnYAxwF8BtNa1WusKZ3tlmwwgWymVAbQDdjrcH0torVcC+5s1XwI8H/j6eWBia8/j9oCeC+wI+b6UFAlsQUqpAmAYsNrZntjmKeBXQIPTHbHRccAe4H8DqaZnlVLtne6U1bTWZcATwDfALuCg1vodZ3tlq2O11rvAP4gDjmntC7o9oCuDtpQpy1FKdQAWArdrrb9zuj9WU0pNAMq11uuc7ovNMoDhwByt9TCgEhMuv90ukDO+BOgH9AbaK6V+5GyvvM3tAb0U6BPyfR5JeknWnFIqE38wf1Frvcjp/thkNHCxUmo7/vTaWUqpvzvbJVuUAqVa6+BV2AL8AT7ZnQOUaK33aK3rgEXA6Q73yU7fKqV6AQT+X97aF3R7QP8YOF4p1U8p1Qb/hMkSh/tkOaWUwp9P3aS1ftLp/thFaz1Na52ntS7A/7terrVO+hGb1no3sEMpFbxD8dnA5w52yS7fAKcppdoF/s2fTQpMBodYAlwX+Po64LXWvqCr7ymqtfYppW4FluGfAX9Oa73R4W7ZYTRwLbBBKbU+0Hav1vpNB/skrPVz4MXAwGUb8BOH+2M5rfVqpdQC4BP8lV1FJOmqUaXUP4CxQHelVCkwA5gFvKyU+in+D7crWn0eWSkqhBDJwe0pFyGEEDGSgC6EEElCAroQQiQJCehCCJEkJKALIUSSkIAuhBBJQgK6EEIkCQnoQgiRJP4/LmvyTrbgCzAAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.scatter(X, Y)\n",
"for r in [0.1, 0.5, 0.9]:\n",
" qr = QuantileRegression(r)\n",
" qr.fit(X, Y)\n",
" plt.plot(X, qr.predict(X), label='quantile'+str(r), alpha=0.6)\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8VOWd+PHPNzdIAhIQVEhIE6viBY2sUfGH+rPYn2jViq5211qX2lq67XbXW6nAruJdWq3a3bauLNZLa738FFHReqnUtloBuUVUpFKCkIBcDAFJINdn/5iZOJM5M3NmMmfOZb7v18uXmWfOzHkGyHee832+z3PEGINSSin/K3C7A0oppbJDA7pSSgWEBnSllAoIDehKKRUQGtCVUiogNKArpVRAaEBXSqmA0ICulFIBoQFdKaUCoiiXJxs5cqSpqanJ5SmVUsr3VqxYsdMYMyrVcTkN6DU1NSxfvjyXp1RKKd8TkY/tHKcpF6WUCggN6EopFRAa0JVSKiA0oCulVEBoQFdKqYDIaZWLUko5aeGqZu56ZR1bWvcxpqKUGVPGMXVCpdvdyhkN6EqpQFi4qplZC9awr6sHgObWfcxasAYgb4K6plyUUoFw1yvr+oJ5xL6uHu56ZZ1LPco9DehKqUDY0rovrfYg0oCulAqEMRWlabUHkQZ0pVQgzJgyjtLiwpi20uJCZkwZ51KPck8nRZVSgRCZ+NQqF6WUCoCpEyrzKoD3pykXpZQKCA3oSikVELYDuogUisgqEVkUflwrIktF5CMReVJESpzrplJKqVTSGaFfBayNevxj4F5jzOHALuDb2eyYUkqp9NgK6CJSBZwLzA8/FmAy8HT4kEeAqU50UCmllD12R+j3AT8CesOPDwRajTHd4cdNQP5OLSullAekDOgich6w3RizIrrZ4lCT4PXTRWS5iCzfsWNHht1USimVip0R+iTgqyKyEXiCUKrlPqBCRCJ17FXAFqsXG2PmGWPqjTH1o0alvGm1UkqpDKUM6MaYWcaYKmNMDfCPwGJjzGXAH4CLw4dNA55zrJdKKaVSGkgd+vXAtSKynlBO/cHsdEkppVQm0lr6b4x5A3gj/PMG4KTsd0kppVQmdKWoUkoFhAZ0pZQKCA3oSikVEBrQlVIqIHQ/dKWUcsjCVc05veGGBnSllHLAwlXNzHrxJXqHLKHkoBKat5/LrAVrABwL6ppyUUqpLFu6dSm3vH0LvUOWANDbdSAg7Ovq4a5X1jl2Xh2hK6VUlmzas4mH3n8IgPbO0N6FnTvOwvQM6TtmS+s+x86vAV0ppQbos87PuGfFPTFtB5g6tm/7YtyxYypKHeuHBnSllMpQT28Pty29LaZt7NCxfGv8t6gb3MysBWvY19XT91xpcSEzpoxzrD8a0JVSKgM3v31zXNuNE28kdP+fzyc+tcpFKaU86rG1j7G+dX1M2+yTZ1NcUBx37NQJlY4G8P40oCullA1vb3mbVz9+NabtqglXUTG4wqUexdOArpRSSTTsaGDh+oUxbZcfdTmHVhzqUo8S04CulMoL6a7a3NG+g182/DKmbfLYyZxWdZrTXc2YBnSlVOAtXBVbcdLcui/hqs2Ong7mLpsb9x5zTpnjfEcHSAO6Uirw7nplXUz5INC3ajM6oFtVrvghkEdoQFdKBV6i1ZmRdqtAfv2J1zO4aHDa58r1hlzRNKArpQJvTEUpzRZBfXj189z89vKYtmlHT6NmWE1G50knteME3ZxLKRV4M6aMo7S4sO9x8fC/UD7mWerGDutrO+mQk5hzypyMgzkkT+3kgo7QlVKBFxkdz319Ma1Fb1BWUkTd2BHUjCwHspcnT5XacVrKgC4ig4E/AYPCxz9tjJkjIg8D/xfYHT70m8aY1U51VCmlMrW3cy8N++dzziSA6r72bE94JkrtOLkhVzQ7I/QOYLIxZq+IFANvisjvws/NMMY87Vz3lFIqJNPJxlxWrsyYMi7nG3JFSxnQjTEG2Bt+WBz+zzjZKaWUipbJZKNVIJ9RP4Oy4jLH+unGhlzRJBSvUxwkUgisAA4DfmGMuT6ccjmF0Aj+dWCmMaYj2fvU19eb5cuXJztEKaXiTJq72DKVUVlRylszJ8e0WQXyiw+/mGNGHuNY/5wmIiuMMfWpjrM1KWqM6QGOF5EK4FkRGQ/MAj4BSoB5wPXALRYdmQ5MB6iuru7/tFJKpWRnstEqkJcXlfPDE3/oWL+8Jq0qF2NMq4i8AZxtjLk73NwhIg8Bln9qxph5hAI+9fX1mqpRSqUt2WTj4k2L+XPzn+Oe89MKz2yxU+UyCugKB/NS4MvAj0VktDFmq4R2c58KvOdwX5VSecpysnHQfo445o/8ubk85th8DOQRdkboo4FHwnn0AuApY8wiEVkcDvYCrAb+2cF+KqXyWOxkYzvDq1+gbuywvjpyyO9AHmGnyuVdYIJF+2SLw5VSyhFTJ1TSsH9++NGYvvYf1v+Q8uJy6xflGV0pqpTyPKsJz3Nrz6X+kJSFH3lFA7pSyrOsAnkBBdxwyg0u9Mb7NKArpTznqXVPsbZlbVy75smT04CulPKMT9o+4YF3H4hr10BujwZ0pZTrjDHcsiRuXaIG8jRpQFdKucoqT+73yhW37lqkAV0p5QqrQH5a5WlMrvZ3RbSbdy3SgK6UyimrQA7BSa/YvSG1EzSgK6VyIlkgX7iqmUlzF7uy5Wy2uXnXIg3oSilHfdjyIU+uezKuPTIid/vGytnm5l2LNKArpRxht3LFzRSFE9y8a5EGdKVU1lmlV6454RoOKDkgrt3tGytnm5t3LdKArpTKGqtAftzI47jw8AsTvsbtGys7YeqESleuLjSgK6UGbCCVK27fWDlINKArpTKWjRJEt2+sHCQa0JVSaXvnk3d4qfGluPZMa8ndSlEEjQZ0pZRtvaaXW5fcGtc+0EVBbi2VDxoN6EopW6zSK9eecC1DS4YO6H2DVofuJg3oSqmkrAL5YRWHcdlRl2Xl/YNWh+4mDehKKUu52nMlaHXobtKArpSKkevNs4JYh+6WlAFdRAYDfwIGhY9/2hgzR0RqgSeAEcBK4HJjTKeTnVVKOefljS+zdOvSuHand0HUOvTssTNC7wAmG2P2ikgx8KaI/A64FrjXGPOEiPw38G3gfgf7qpRyQFdvF3csvSOuPVfb2WodevakDOjGGAPsDT8sDv9ngMnA18PtjwA3oQFdKV+xSq9cd8J1DCkZktN+BLYOfe92WPrfMObv4KjzHD+drRy6iBQCK4DDgF8AfwNajTHd4UOaAMu/DRGZDkwHqK6uHmh/lVJZYBXIxw4dy7fGf8uF3gRMTze8dF1s2/rXvBPQjTE9wPEiUgE8CxxldViC184D5gHU19dbHqOUyg0v3i0oMIuKGv8M7z0d3z7hcqiqz0kX0qpyMca0isgbwESgQkSKwqP0KmCLA/1TSmVBdCDfuLONhs27ae/s5sD9l7k6+ej7RUWd7fDKLOvnzpwDZSNy2h07VS6jgK5wMC8Fvgz8GPgDcDGhSpdpwHNOdlQplb6n1j3F2pa1fY837mxjWWMLbVtC29k2424A9e2iovcWQOMf49sPPQOOSbxVsNPsjNBHA4+E8+gFwFPGmEUi8gHwhIjcBqwCHnSwn0qpNHT0dDB32dy49r++fzZt/Wq+3QygvlpU1PYpLI6/AxMAU+6EkrLc9seCnSqXd4EJFu0bgJOc6JRSKnNWefIZ9TMoKy7j4edetHyNWwHUF4uKXrjKuv3Yr0HNpNz2JQVdKapUQFgF8poDaph2zLS+x14LoJ5dVPTRa/DhIuvnvvJTKPRm6PRmr5QKgFxVb6RTueJkAM3k83pqUVFvL7x4jfVz4/8eak/PbX8yoAFdKQc4Vb0RHTSHVz9P3dhh1IwsjzkmWQmiUwF0IJ/X9UVFSx+A7R9YP3fefSCS2/4MgAZ0pRzgRPVGJGh2D32DkkN20NYJyxpb2LG3gy279rNr01cZU1FK3eDmpOdwIoD6rlql4zN49T+snzv2Eqg5Nbf9yRIN6Eo5wInqjZ+8+i69Bz5NQVRbd6/h/Yaz+lb1uVXH7ZtqlUQTnADn/yx3/XCIBnSlHJDtyceb376ZlkGbYto6tp0PpjjuWDdGxl6bbI2xcz28/V/Wz530XTj46Nz2x0Ea0JVyQLYmH6MnPMtKimjr7Ka3cxRdLaclfV2uR8aerFbJcDTu560INKAr5YCBTj5aVa7UjR3GkmWT6YgKmoL1Jkq5Hhl7plpl9eOweYn1c2fMgqGHJH2537cikNDuuLlRX19vli9fnrPzKeU3qUoQ+48ev3TkKJ5Z0Rw3Mr7zomN9EYCyJku58UlzF1umjiorSnlr5uRMepYVIrLCGJNyhy8doatA88vl893v3E1bd1tce/8SRKsKlfovjPDFZ8y6RdeA6bV+7uwfQ/HgtN/SN5O7CWhAV4Hlh8vnPZ17uHfFvXHt6Wxn63odd4Yy+rLt7oTfzUj8/AArVTw9uWuDBnQVWF6vjbZKr8w+aTbFhfGVK0GTzpftwlXNsOhq2ju6KRtURF3VMGqjF1NlcfGPJyd306ABXQWWVy+frQL5MQcew8VHXOxCb9xh68t2dzONC+bQ0dhCT29orq+9o5tljS3sOWAcdX+fZKSeIc9M7mZIA7oKLK9dPnvxbkFuSfplGzXB2dC0uy+YA8zuvhKAyo9Kecuhvvk1hQUa0JWPpJtz9crlc7YCuV8meO0YVlpM676uvscXFvyZEwvWUTaoCBjT197e0c2i3on8pXd8zOvdvsryKg3oyhcymeB0+/L59iW30913H/XPZTIi98MEr10LVzXT1hn6c7mjaH5fe4EIdVXDPj/w/J/xi/etywj9MkmZaxrQlS9kOsHpxuXzrv27+M9V/xnXPpDUitcneNOy6Gpulu646FNUKNReMBtGHNrX5pWrLL/QgK58wasTnP1ZpVf+/eR/p6hgYL9qfvn8CRkDi64GQmmU/mZ3X4l0w8VRwRzcv8ryGw3oyhe8NsHZn1Ugrz+4nnMPPTcr7+/1z5+QxQrOskFFtHd0c2v3N9jH54t/En0WNycp/TZvoQFd+YJXL71zVbni1c9vad8uGp+YQUPTbuva8fPu49YFa9iHtz+LH+ctUgZ0ERkLPAocAvQC84wxPxORm4DvADvCh842xrzkVEdVfvPapXeuSxC99vkthUfjjTvbWNavdvwrGy7mzuND+8tMDR/u6c+CP+ctUm7OJSKjgdHGmJUiMhRYAUwFvgbsNcbcbfdkujmX8rtsBHK/XcYntf51WPt8TNPC1Vv68uSRunFwf4OrdNXOfNFyJ0sBGudmJ5VmV9Y25zLGbAW2hn/+TETWAj7916dUZra3b+f+hvvj2jOpJffqZXxaXzRJdje8pu2bloHQNxO4YX6ct0grhy4iNcAEYCkwCfiBiPwTsBy4zhizK9sdVMptVqPyGybeQIEUWBydnFcv42190STbovbg8XDSdwAYE5DacV/NW4TZDugiMgR4BrjaGLNHRO4HbiW0v/6twE+Bb1m8bjowHaC6ujobfVYqJ6wC+SmjT+GsmrMyfk+vlh8m+qK57qkG6lbPid0MK5rF7oZ+DIRWfDFv0Y+tgC4ixYSC+WPGmAUAxphtUc//D7DI6rXGmHnAPAjl0AfaYaWc5uSEp1cv4/t/oUSv4FzWGNrJsC+on/gdOCR2KX40PwbCRPy2r4udKhcBHgTWGmPuiWofHc6vA1wIvOdMF5XKjVxUrqQavbo1YTqmopQtrW3cXvSruOd6eg0NTbupvWK+xSut+S0QBoWdEfok4HJgjYisDrfNBi4VkeMJpVw2At91pIdKOSyXJYjJRq+uTZi+cBW/ObiNZe0t9PS7AdDN3f9EByVIN33lhupzXqtYslPl8iahSp3+tOZc+Vrz3mbmr4kfdTq9nW2i0WtOJ0xbNsBbn+e/I+mUJRtaMMbElBuC+ykhL/JixZKuFFV5yWpUfuPEG5Es3fkmEzmZME1SqVJ7xXwawkEKj6/i9AIvVixpQFd5xSqQf2nslzi96nQXehPLsQnTpQ/A9g8SPx9VqRKkCU2nebFiSQO6ygt+uFtQ1sv9ktWNJ7mZsk5o2uPFiiUN6CrQ/BDII7IyOk4WxIfXwqlXZ9w/r00Aus2L9fYa0FUg+SmQR8t4dJzhaNwuL04Aus2L6SkN6CpQNuzewK8/+HVcu5cDecYj32RB/OTvwUFHZq2PXpwA9AKvpac0oKvA8GLlSippj3y79sPL1yd+wyyMxq14cQLQKX5OLWlAV75nFcin1Exh4uiJLvQmPbZHvslG41/5KRQ6+6vsxQlAJ/g9taQBXflW/0C+cWcbDZt3s2vTV3m1op0ZU5o9/0uYdOTbvAJWPpr4xQ6Nxq14cQLQCX5PLWlAV75jNSLfuLONJcsmuzayyvQy3Wrke0fRfMoGFcHKMfEvyGEQj+bFCUAn+D21pAFd+UayypVJcxezryv2ly5XI6uBXKZHRr4zzXyKw6szCwuEuqphsQe6FMijeW0C0Al+Ty1pQFeet6F1A79em7xyxc2R1U3Pv5/xZfrUpp9Qd2gbDU1CewexN1T2QBDPN35PLWlAV55mNSq3KkF0a2S1cFUzrfu6LJ9L+GXSb4KzdmT553uNjzoKJv5zNruYVX6uALHD76klDejKk6wC+SVHXMLRBx5tebxbI6u7XlmX8Lm4LxOHF/84ze8VIHb5ObWkAV15SqYrPN0aWSVL6cyYMi55ED/1Ghhek/1OOcTvFSD5QAO68oRUgdzOpb4bIyurVM8Q2rlp8BNMbaqyfpEPRuNW/F4Bkg80oCtX2RmRe+FSP9EXSnSqJ3IfzsIC4YQvjIh9g/PuAw+vWLXD7xUg+UADunLFhy0f8uS6J+ParVIrbl/qJ/1CGbqWukOfpqFpd3yVCmR9NO7mpKTfK0DygQZ0lXNWKzz/+v7ZbGndx6t/XBwXpLJ5qZ9JQLT6QrnBPACLiuD4MbFVKuBYSsXtKxW/V4DkAw3oKmes0iuj+X88sWxP36IgqyCVrUv9TANi5IsjklKJaO/ojj3Q4dy421cq4O8KkHygAV05LvUKz+RBKluX+pkGxHvLH44P3oTSK7mc4Mz2pGTQa8q9oLfX8MHWPYwoL8nJXEPKgC4iY4FHgUOAXmCeMeZnIjICeBKoATYCXzPG7HKuq8pv7Ex42glS2brUT3Su5tZ9TJrbL9UTVW5YVzWMZY0t9PQaAN7uPZrfF57Gnecdm9b5Byqbk5Jup2+CavnGFp5Z2Wz53J0XOf/vxc4IvRu4zhizUkSGAitE5DXgm8Drxpi5IjITmAkk2ahZ5Yt0asntBqlsXOonOhdEAtq71K2eE5sPh77H39j29b4vlDtdGM1mc1LSC+kbv9u2Zz8/X7ye7vAXfSJfHFXOOceOzkmfUgZ0Y8xWYGv4589EZC1QCVwAnBE+7BHgDTSg57V1Let4Yt0Tce3JFgVlM0ilSiFYnQuicuMGGpqKYgP6l/4DhoyiFngr7R4NjNXnufOiY7OSJtGa8vR0dPfw5DubWbv1s5THHlBaxHdP/yIjykty0LNYaeXQRaQGmAAsBQ4OB3uMMVtF5KCs9075gjGGW5bcEtdu57Zv6aZTEgVtOymE6HN1tm7lmqKn496/L1fu8uKfRJ/nzouO5a2Zkwf8/lpTntyNz71HV0/ykXfEZSdXM75yWOoDc0CMsddpERkC/BG43RizQERajTEVUc/vMsYMt3jddGA6QHV19Qkff/xxdnquPMEqvXLlsVdSOST7l+39gxyERvORUatVgKqsKI0NgOHc+MLVW+ImOmd3Xxl/vEsmzV1s7/NkKNmfZb6lXN78aCcvrtlq69iJh47g/OPGUFCQ20ViIrLCGFOf6jhbI3QRKQaeAR4zxiwIN28TkdHh0floYLvVa40x84B5APX19fa+PZTnZbrnykAky/smm/Bc8fz9nCAfxrRHT3TO7r4S8NYiGadTIvlaU759z37u/f1Hto+/+suHc/ABgx3sUXbZqXIR4EFgrTHmnqinngemAXPD/3/OkR4qT8lWIM+kZC5ZkEt05x+A9cuFEbUjYnLjtVfMpyHcB/FYQFu4qpkCEXosrp6zmRIJek25MYbZz75n+/jjqoZx6UnVDvbIeXZG6JOAy4E1IrI63DabUCB/SkS+DWwCLnGmi8oLsjkiz7RkLlneNzLheYN5IO75nl5DQ9PuuJtGeDGgRf5srIK5l64gnJbJF37/vPfGnW00NLXS3tlDWUkhdVUV1ER9qd96wTEUFRY49hncYKfK5U0gUcLozOx2R3nNezvf45mPnolrH0hqJdOSuWQVMZE7/7z9t/jXze6+EumGqeefm3Gfc8XqzwagUCRv8tt2vvA/2LKHXy9JPB+3cWcbyzZ+vnagvbOH1Ztb+YcTxwb6z1BXiipLA6lcSSXT/HD/vO+95Q+HNsJq+h0QqhcPbZLVzcu9J/Kn3rq+1/qleiPRn0GvMYEORNESfeHPfnYNSxtbbL3Hppb2vmAe/R5Br7PXgK7iWKVXvl/3fUaVjcrK+2dSMhe5BN/a2sY95Y9SN24YtSPHxB943n3cumAN+3r9uSOgX8sJs7mNQKLFX+2d8VcuALdPHR9XdVI780XLY4NeZ68BXfWxCuRDi4dybf21WT1PuouJFq5qpuPZf+Vfeg0UQXsHLAuP1GpHlsNZt8GgoQBMDb/Gr9UbftyidiDbCDy9ookVH8fuGFJWUmgZvMtKCrl84hc4eswBKfvk1y/GgdKArhwrQUw0arNdMrfrY3jzHli9Je7yuafX8I1tX2fG8eO46953Ep7DSZmOSpO9zo/lhHbnRPZ19nDLog9Svl9dVUVM/htCX2p3XHisrWAO/vxizAbbC4uyob6+3ixfvjxn51PJOVlLPqCFK/3uw/nbpZv6fo7UjEe/pxuLYzL9fEFb0LNwVTNXP7k64fPplAHeceF4JHxXp2ykcIK0m6TdhUUa0PPQym0reWHDC3Ht2RyRJ6qjTrjScflDsNU6MCxcvYWr274Z116Y7jmyKNOVnE6vAM0lqy+naGUlhVxwvHUAPffY0Zx6+EgnuxcoWV0pqoLBycqV/r/cVoEWLCal+o3GY0RqxquaKbUY1SYKJLmY+Mq0UicIm2Lt3NvBT1/9K8+tbk74d1AAdPcaHl+2ibKSULrEr6NjP9GAnies0iv/NuHfGD44bvudjCSqn+5vTEVp8iA+dDScMTOmKVFeOdH+LbmY+Mp00s2Pk3WRCc7+ElWdABQWCp3dvX3H6V7ruaEBPeCsAvlBpQfxveO/l9XzJCo1i/bjkgc56eARQHn8kyl2N0w00enWxFemk25en6xLFLytJKpGKRSJ26kwH2rAvUADekDlevOsRPnsyH4qZYOKqKuK3U+FE78Dh4zP+JxuVoRkem4vVbGs3/4ZD765Ma3XRN915+TaEZZfTm6mwvKdTor6XP+Z/COOeTlmv4oIJ3dBBKiJWshRRDe3FD3c9/jrJ/erdHB5r/F8lc7oG+zdMs2qksT2VsbKNp0UzQPRE5GFpRvZOXglrY2hsq9IUHc6kEdUVpTyL3v/K669bFD4n9hX7obC4pz0RaUfvC89aSzHVVWkPrAfr6XC8p0GdB8LTUR2MeiQhX1t3b2Ghs27eej8n+SmEzvWwZJf8puD21jWLjGLQQoLBM67DzRv6qglGz7ludVb0nqNkzcs9lJaKd9oQPexTwc/xqBDYts6tn+Fzt7Bzi+q6FepEsmNNzTt5pq2bzpyziAtFBkIJ1In2ebFrYnzgQZ0H4pMeJaVFNHWGbqNWm/HaLp2nQLAsNLijPfWSGrZ/8C2xDcMqL1iPrV8vp9KNg1kvxA/Szd4f+e0Wg4dNcSh3iiv04DuI/0rV+rGhm6j1rblwr620uJCRMhov/GE7Cz+cVime6j7idVGVam4MfpW3qUB3QcSlSA+dP5PLNMQ1yTYWyOtsrFkQfwLk+C4r9l/ryxIds/Q2pkv+i4Fk+7t0UCDt0pNA7qH2dlzxSpXOaAVlB4YjVtJtMISwOD9FEy6qRO/3ZxYeYMGdA/q6e3htqW3xbXbLUFMezVisiB++o9gmPsB0uoz9eeVFMwv/rCepl3pLaLR0bfKBg3oHmOVXplRP4Oy4jLb72GrbKxrH7w8M8E7MKDRuBPVKP0/U6LlcLlejdjTa/iPhZo6Ud6gAd0jrAL5xNETmVIzJaP3S1g2lmw0fu69UDCwu6A7VY3S/0uiraOb1n1dccc5vclVuqmT2V85kqGDdUGVyo2UAV1EfgWcB2w3xowPt90EfAfYET5stjHmJac6GWQ52XNlawMs/1Xi5zMcjSda9m1VjXLdUw1AZkHd6kuiuFAoLhC6+t3VJpurEdMN3qCjb+UuOyP0h4GfA4/2a7/XGHN31nuUJ3ISyB2c4Ew0Ek+U4+4xJuORutWXRFePYXhZMWUlRVlJ7di9PVo0Dd7Ka1IGdGPMn0Skxvmu5Id3PnmHlxrjL2ZSBXLbeekl/w071iZ+oyxVqiQaiSfadTHyfCaTlony4q3tXay68ay03isi3dH3LRccQ3HhwNJRSjltIDn0H4jIPwHLgeuMMZYrIkRkOjAdoLra/v0Fg6art4s7lt4R125nRG4rL53jcsNEQbbHmKRbqNrZN72/gd4UIt3gfUBpEbPOOSqt1yjlBZkG9PuBWwmVAN8K/BT4ltWBxph5wDwIbZ+b4fl8zSq9MuukWZQUlth6faLRMIuuhqYx1i8afzHUnpZ2X+1KFGQrw1cP1z3VYDlSF0JfUOmM0tMpw4zcHi0dmjpRQZFRQDfGbIv8LCL/AyzKWo88LN1yPKtAPnnsZE6rSi/Q9h8NR24a0d5h0ceqH4X6+NYexlQsdmz1ZLIgGznfNU+ujisvNJB22iVZGWa6o+/oO8srFTQZBXQRGW2M2Rp+eCGQXiGuD6VTjpftCc8xqfYa//JNUDo8pxtYpap1nzqhkquzsQVB1PmWNrb0PV7a2BLz2MphBw3h26fWpn0upfzKTtni48AZwEgRaQLmAGeIyPGEBlwbge862EdPsLM5VNYrVzr2wqv/bnuv8VxvYJVqi9TKAeS+N7e088s3/pZWfzR1ovKdnSrLQ/fXAAAMM0lEQVSXSy2aH3SgL56WaFS5pXUff2n+C69tei3uuYwDedK9xqcxpqLMMpWSrI9uSCf37Yc9vpXyOl0papP1JGA3Q6sW8dqmqpjWjAL5J2vgnfmWTzXubOMb277OlrbkufuBVoNkW6K0jJ10SbRzxh/C6UeMcqqbSgWGBnSbZkwZx4ynG+jqCaU9Bh2yAICuHti4s42akeXMPnk2xQVpLvNOUW64cFUzs5atYV9XKFAny4unvSlXDhxx8NCYYGwnkOvo2zl616dg04Bu09QJlcxe8C4Fo56Jae818MFHh/HQ+dPtv9nSebD9fevnxvwdnDCt72E6eXEv3MtRUyfela93fconGtBtuuKFH9Fz4Kdx7R2fXMR2u2+SweKfZDd2sKrnzuW9HNMN3pedXM34ymEO9Ualkg93fcp3GtBTePT9R2nc00jD5t0x7R2fXNT3c9IcdbIgfvL34KAjk54/2Y0dcjm6WtbYwrOrmtN6jY6+vcVrk+Yq+zSgJ7ChdQO/Xvvrvsft4ZsxRwfyiLgctTGhVZyJpLEUP9mNHZwcXSUbfW/c2UZDUyvtnT2UlRRSV1XB49MnZr0PKru8Nmmusk8Dej/dvd3cvvT2uPYD919m+ctQUVpsbz+Vs+dCcfq/OJH3zuYinf7SSZ1s3NnGyk276OjuBaC9s4fVm1vTXs6vcs+Lk+YquzSgR7FaGHTDxBsokALqBjdb/jJcfMwQFt56Ce0d3ZQNKqKualhf3TiQlY2xpk6oHNh9QqO8+v4n/GHdjtQHRolOnUyau7gvmEdoHtYfvDBprpylAR3rQH7lsVdSOeTzf+hTJ1Sy/OMWHl+6mR5juLP4QQ4aWsLOdztpD6/gbO/o5u2/fcqX1l3Yt0nV1Cz1MZPRlRN3ltc8rL/lctJc5V7eBvSFq5q55e1baO/spqykiLqxw6gZWc7Jo0/m7JqzLY9/e8Uqbi18Fgilybftid0da3b3lX0/9y8JG2j9r53RVbpVJzPPOZJhpenVzWseVinvEpPgZgROqK+vN8uXL8/Z+RKZ/ft5PPfBMrqj9kYpKhAurL6KP3y4Iz5gvnAVC1dvob2jO/69ooK4lchI3Wp0fedFx2Y8WnpudTNLNthfbVlcKNxywfiMzhWtfy0zDPyzKKWSE5EVxpj6VMfl1Qi9cXcjj37wKC//dUtMMO/45CI6gMe2bOrb7rVmzzI6nv05jatHUDuyPCaY/7ZnMu+ZQ22dc0vrvgHX/3rpzvKah1XKu/IioPe/W1CiEkTD53uNA/T0hjbEqh1ZTtmgIq5u+2ba5x5TUZp0cdCkufF7lqebOrnpq0czqKgw7b5lSvOwSnlT4AO61YTniP1fZ0vr/r7H/1i4mONkg+Xrb2q7mKnn/wNUNVOa5CbIViKTlokqVCAU1K/7/w08+c5maqKrYxI4rmoYl56Un7fy031IlEouEAHd6he9YX/8zoXfr/s+o8pGhUsQ3+UGMy/he0Zy45Xhyb5I4Eh0a7UIITTSr+wXcBItDoJQSqWhqdUyoOtqyxDdh0Sp1Hwf0Pv/ou8c/BjX/xEGFRVwwheGUzOynFMrT+XM6jNDL3jhKqYCdYe20dBU1Fc/XlkxmGnbL6W96/Ng3b8sMFFwThbEI7sLHj+2om91pZX2zh5unzqeggK9PZoV3YdEqdR8H9Ajv+hFw1ZQWPpxX3tHdy/LGlu4oOoqzjx4mOVNI2pHlsOoI2Hi9wC4w8YlfaJJwW179rPx0/aEe33XjCynZmQ5z61utgzqlRWlGsyT0Pp3pVLzfUDfuncbgw75fVx7xycXMadofmhPlaYx8S+0WMFpd7LvrGMOjru/ZSqR1MnJtSN0+XUGtP5dqdR8G9AXrNzEbUtup3hkbG34qG2n8v3C5yFcrRJTO37MhXDoGWmfK7rqxGpjqv6572R3lk9W9qeTfonpPiRKpebLgP7Vx6/lr9t2Ez01eePOTxEECp+PObZsUFFa+6n8Zf1OXnh3q+VzG3e2sWxjS9/Nmts7e1i5aRf/cOLYtFd99j9eJ/2S0/p3pVJLuVJURH4FnAdsN8aMD7eNAJ4EaoCNwNeMMbtSnWygK0V/8dYbPLDqsb7NoQ6gnR982k6xKYg79v6e89lZNIa/P6HSevUnoXr0WxettX3+l9ZsZfe+rrj2yopS3po5OcNPFTJp7mLLlEI23lsp5W/ZXCn6MPBz4NGotpnA68aYuSIyM/z4+kw6msrCVc385JW1fDr48b5qknGymfNbBlHeUwTEBvNIuWGhCJeeUMkzK5pjRr3p1Hz3Lxl8Ytkmy+OyMTGnk35KqYFKGdCNMX8SkZp+zRcAZ4R/fgR4AwcCeiQNsb93N0MGd3O4NHFOSxnlPbHB+KbuaXTy+SZTpcWFXHZyNb9dtimu1C1RzfePpoxjeHlJ0v44OTHn1Uk/zesr5R+Z5tAPNsZsBTDGbBWRg7LYpz6RksTh9DJj5x7ggL7nHu6Zwl/N2LjXRCYq2zp7ktZ8Z7Jgx8mJOS9O+mleXyl/cXxSVESmA9MBqqvTW7IeSTe0UUqD+SKbzEG83Xs0oaU8sQ47aAgn1oyIaSsrKUxY850JJyfmvDjpp4t5lPKXTAP6NhEZHR6dj4bEN743xswD5kFoUjSdk0TSEJ0U82TPlyyPKSkq4ITq4X0plBvOO4qyktDHcqLm28mNqby26ZXm9ZXyl0wD+vPANGBu+P/PZa1HUZLdIHnIoCKuO+sIrphUm/D1Xhz1+olX8/pKKWspA7qIPE5oAnSkiDQBcwgF8qdE5NvAJuASJzoXHZCbW/fF7ZVi9z00gGfGi3l9pVRieXnHImWfVrko5T69Y5HKCr3CUco/NKArpfJSEK8+NaArpfJOUNdYxG+CopRSAZdsjYWfaUBXSuWdoK6x0ICulMo7idZS+H2NhQZ0pVTemTFlHKXFhTFtQVhjoZOiSqm8E9RV5BrQlVJ5KYhrLDTlopRSAaEBXSmlAkIDulJKBYQGdKWUCggN6EopFRA53T5XRHYAH2f48pHAzix2xy/0c+ePfPzMoJ/bji8YY0alOiinAX0gRGS5nf2Ag0Y/d/7Ix88M+rmz+Z6aclFKqYDQgK6UUgHhp4A+z+0OuEQ/d/7Ix88M+rmzxjc5dKWUUsn5aYSulFIqCc8HdBE5W0TWich6EZnpdn9yQUTGisgfRGStiLwvIle53adcEpFCEVklIovc7kuuiEiFiDwtIh+G/95PcbtPuSAi14T/jb8nIo+LyGC3++QEEfmViGwXkfei2kaIyGsi8lH4/8MHeh5PB3QRKQR+AZwDHA1cKiJHu9urnOgGrjPGHAVMBP4lTz53xFXAWrc7kWM/A142xhwJ1JEHn19EKoF/A+qNMeOBQuAf3e2VYx4Gzu7XNhN43RhzOPB6+PGAeDqgAycB640xG4wxncATwAUu98lxxpitxpiV4Z8/I/TLHax9PhMQkSrgXGC+233JFRE5ADgdeBDAGNNpjGl1t1c5UwSUikgRUAZscbk/jjDG/Alo6dd8AfBI+OdHgKkDPY/XA3olsDnqcRN5EtgiRKQGmAAsdbcnOXMf8COg1+2O5NChwA7goXCqab6IlLvdKacZY5qBu4FNwFZgtzHmVXd7lVMHG2O2QmgQBxw00Df0ekAXi7a8KcsRkSHAM8DVxpg9bvfHaSJyHrDdGLPC7b7kWBHwd8D9xpgJQBtZuPz2unDO+AKgFhgDlIvIN9ztlb95PaA3AWOjHlcR0Euy/kSkmFAwf8wYs8Dt/uTIJOCrIrKRUHptsoj8xt0u5UQT0GSMiVyFPU0owAfdl4FGY8wOY0wXsAD4Py73KZe2ichogPD/tw/0Db0e0N8BDheRWhEpITRh8rzLfXKciAihfOpaY8w9bvcnV4wxs4wxVcaYGkJ/14uNMYEfsRljPgE2i0jkDsVnAh+42KVc2QRMFJGy8L/5M8mDyeAozwPTwj9PA54b6Bt6+p6ixphuEfkB8AqhGfBfGWPed7lbuTAJuBxYIyKrw22zjTEvudgn5ax/BR4LD1w2AFe43B/HGWOWisjTwEpClV2rCOiqURF5HDgDGCkiTcAcYC7wlIh8m9CX2yUDPo+uFFVKqWDwespFKaWUTRrQlVIqIDSgK6VUQGhAV0qpgNCArpRSAaEBXSmlAkIDulJKBYQGdKWUCoj/BeyohD1jWg0DAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# 和 StatsModel 的包对比一下\n",
"import pandas as pd\n",
"import statsmodels.formula.api as smf\n",
"\n",
"data = pd.DataFrame(np.hstack([X, Y]), columns=['x', 'y'])\n",
"mod = smf.quantreg('y ~ x', data)\n",
"\n",
"plt.scatter(X, Y)\n",
"for r in [0.1, 0.5, 0.9]:\n",
" res = mod.fit(q=r)\n",
" y_ = res.params['x'] * X + res.params['Intercept']\n",
" plt.plot(X, y_, label='quantile'+str(r), alpha=0.6)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Ref\n",
"- [这份 tutorial 的前几页](https://support.sas.com/resources/papers/proceedings17/SAS0525-2017.pdf)\n",
"- [StatsModels Example](https://www.statsmodels.org/dev/examples/notebooks/generated/quantile_regression.html)"
]
},
{
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment