Created
July 12, 2017 15:20
-
-
Save rutj3/66b3edc665c0897ea4e777e54ec6a985 to your computer and use it in GitHub Desktop.
Numba_regression
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 23, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"'0.33.0'" | |
] | |
}, | |
"execution_count": 23, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"import numpy as np\n", | |
"import numba\n", | |
"import matplotlib.pyplot as plt\n", | |
"from scipy.stats import linregress\n", | |
"%matplotlib inline\n", | |
"numba.__version__" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 30, | |
"metadata": { | |
"collapsed": true | |
}, | |
"outputs": [], | |
"source": [ | |
"@numba.vectorize\n", | |
"def calc_slope(x1,y1,x2,y2):\n", | |
"\n", | |
" xd = x2-x1\n", | |
" \n", | |
" if xd == 0:\n", | |
" slope = 0\n", | |
" else:\n", | |
" slope = abs(y2-y1) / abs(xd)\n", | |
" \n", | |
" return slope\n", | |
"\n", | |
"@numba.njit\n", | |
"def siegel_repeated_medians(x,y):\n", | |
" # Siegel repeated medians regression\n", | |
" \n", | |
" n_total = x.size\n", | |
" \n", | |
" slopes = np.empty((n_total), dtype=y.dtype)\n", | |
" ints = np.empty((n_total), dtype=y.dtype)\n", | |
" slopes_sub = np.empty((n_total-1), dtype=y.dtype)\n", | |
" \n", | |
" for i in range(n_total): \n", | |
" for j in range(n_total): \n", | |
" if i == j:\n", | |
" continue\n", | |
" \n", | |
" slopes_sub[j] = calc_slope(x[i],y[i],x[j],y[j])\n", | |
" \n", | |
" slopes[i] = np.median(slopes_sub)\n", | |
" ints[i] = y[i] - slopes[i]*x[i]\n", | |
" \n", | |
" return np.median(slopes), np.median(ints)\n", | |
"\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 25, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"np.random.seed(0)\n", | |
"n_samples = 200\n", | |
"\n", | |
"x = np.random.randn(n_samples)\n", | |
"y = 5 * x + 2 + 0.1 * np.random.randn(n_samples)\n", | |
"y[-20:] += -20 * x[-20:]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 26, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"xnew = np.linspace(x.min(), x.max(), 10)\n", | |
"\n", | |
"# scipy linregress\n", | |
"slope, intercept, _,_,_ = linregress(x,y)\n", | |
"y_scipy = xnew * slope + intercept" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 31, | |
"metadata": { | |
"collapsed": true | |
}, | |
"outputs": [], | |
"source": [ | |
"slope, intercept = siegel_repeated_medians(x,y)\n", | |
"y_sieg = xnew * slope + intercept" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 32, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"<matplotlib.legend.Legend at 0x21b250dd8d0>" | |
] | |
}, | |
"execution_count": 32, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAGoCAYAAADVd+V5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtclGX6+PHPMweG44AHkIFMUDRRUEtbjTIlM60sETHL\nStvadju52Zr2w2qtzNxM13a1ctstszKskHC/WlurSaFhZpsHPIUHPHLyxFFgmHl+fwwzMDKYKTDA\nXO/Xi5c5zz3P3EMyF/fhum5FVVWEEEIIT6VxdweEEEIId5JAKIQQwqNJIBRCCOHRJBAKIYTwaBII\nhRBCeDQJhEIIITyaBEIhhBAeTQKhEEIIjyaBUAghhEeTQCiEEMKjSSAUQgjh0XTu7sDlUhRFAcKA\nUnf3RQghhNsFACfUX1FIu80HQmxB8Ji7OyGEEKLVuAI4frGN20MgLAU4evQoRqPR3X0RQgjhJiUl\nJXTt2hV+5QxhewiEABiNRgmEQgghfjXZLCOEEMKjSSAUQgjh0SQQCiGE8GjtZo1QCCHcyWKxYDab\n3d2Ndk2v16PVapv8vhIIhRDiMqiqSn5+PmfPnnV3VzxCUFAQoaGh2FLIm4YEQiGEuAz2IBgSEoKv\nr2+TfkCLOqqqUlFRQWFhIQAmk6nJ7i2BUAghLpHFYnEEwU6dOrm7O+2ej48PAIWFhYSEhDTZNKls\nlhFCiEtkXxP09fV1c088h/173ZTrsRIIhRDiMsl0aMtpju+1TI02M4vFQmZmJnl5eZhMJoYOHdos\nu56EEEJcGhkRNqO0tDSioqKIj49n0qRJxMfHExUVRVpamru7JoQQopYEwmaSlpZGUlISsbGxZGVl\nUVpaSlZWFrGxsSQlJUkwFEK43dGjR3nwwQcJCwvDy8uLbt268eSTT3Lq1ClHm+HDhzNt2rRG7/HN\nN99w00030bFjR3x9fenZsydTpkyhurq6Jd5Ck5BA2AwsFgvTp09nzJgxpKenM2TIEPz9/RkyZAjp\n6emMGTOGp59+GovF4u6uCiE81MGDBxk0aBA5OTmkpKSwf/9+li5dyvr167nuuus4ffr0L95j9+7d\njB49mkGDBvHtt9+yc+dOFi9ejJeXV5v6fJM1wmaQmZlJbm4uKSkpaDTOv2toNBqSk5OJi4sjMzOT\n4cOHu6eTQohmoaoq58zuCQI+eu1FbyZ5/PHH8fLy4quvvnKkJVx55ZVcffXV9OjRg2effZa33nrr\ngvf46quvCA0NZf78+Y7HevTowejRoy/9TbiBBMJmkJeXB0BMTIzL6/bH7e2EEO3HObOFPn/+0i2v\nvfulUfh6/fLH+unTp/nyyy+ZO3euIwjahYaGcu+99/Lxxx/z5ptvXvA+oaGh5OXl8e2333LjjTde\nVt/dSaZGm4G94kF2drbL6/bHm7IyghBCXKycnBxUVSU6Otrl9ejoaM6cOUNRUdEF7zNhwgTuuece\nhg0bhslkYty4cSxZsoSSkpLm6HazkRFhMxg6dCgRERG88sorpKenO02PWq1W5s2bR2RkJEOHDnVj\nL4UQzcFHr2X3S6Pc9tq/hqqql/V6Wq2WZcuW8fLLL/P111/z/fff88orr/Dqq6+yZcuWNvPLvowI\nm4FWq2XhwoWsWbOGhIQEp12jCQkJrFmzhgULFkg+oRDtkKIo+Hrp3PJ1seuDUVFRKIrCnj17XF7f\ns2cPHTp0IDg4+KLuFx4ezv3338+SJUvYtWsXlZWVLF269KK/Z+4mgbCZJCYmkpqays6dO4mLi8No\nNBIXF0d2djapqakkJia6u4tCCA/VqVMnRo4cyZtvvsm5c+ecruXn57NixQomTpx4SVVcOnTogMlk\nory8vKm62+xkarQZJSYmMnbsWKksI4RodZYsWUJcXByjRo3i5ZdfJjIykl27djFjxgzCw8OZO3eu\no21RURHbtm1zer7JZCI9PZ1t27Yxbtw4evToQWVlJe+//z67du1i8eLFLf2WLpkEwmam1WolRUII\n0er07NmTrVu3Mnv2bO666y5Onz5NaGgoCQkJzJ49m44dOzrafvTRR3z00UdOz58zZw633347Gzdu\n5JFHHuHEiRP4+/vTt29f0tPTGTZsWEu/pUumXO5iqbspimIEiouLizEaje7ujhDCg1RWVnLo0CEi\nIyPx9vZ2d3c8woW+5yUlJQQGBgIEqqp60VtXZY1QCCGER5NAKIQQwqNJIBRCCOHRJBAKIYTwaBII\nhRBCeDQJhEIIITyaBEIhhBAeTQKhEEIIjyaBUAghhEeTQCiEEMIlRVFIT09v0dfMyMhAURTOnj3b\nYq8pgVAIITxUUVERjz76KFdeeSUGg4HQ0FBGjRrFpk2bAMjLy+PWW291cy+bnxTdFkIIDzV+/Hiq\nq6tZvnw53bt3p6CggPXr13Pq1CkAQkND3dzDliEjQiGEaEqqCtXl7vn6FYconD17lszMTF599VXi\n4+Pp1q0bv/nNb0hOTubOO+8EGk6NHj16lLvuuougoCA6duzI2LFjyc3NdVyvqanhj3/8I0FBQXTu\n3JlZs2YxZcoUEhISHG2sVivz5s0jMjISHx8f+vfvT2pq6uV/3y+DjAiFEKIpmSvglTD3vPasE+Dl\nd1FN/f398ff3Jz09nSFDhmAwGC7Y3mw2M2rUKK677joyMzPR6XS8/PLLjB49mh07duDl5cWrr77K\nihUrWLZsGdHR0fztb38jPT2d+Ph4x33mzZvHhx9+yNKlS+nZsyfffvst9913H8HBwW47ukkCoRBC\neCCdTsd7773Hww8/zNKlS7nmmmsYNmwYd999N/369WvQ/uOPP8ZqtfKvf/3LcXL9smXLCAoKIiMj\ng1tuuYXFixeTnJzMuHHjANvhv59//rnjHlVVVbzyyiusW7eO6667DoDu3buzceNG/vGPf0ggFEKI\ndkHvaxuZueu1f4Xx48dz++23k5mZyebNm/niiy+YP38+//rXv3jggQec2m7fvp39+/cTEBDg9Hhl\nZSUHDhyguLiYgoICfvOb3ziuabVaBg4ciNVqBWD//v1UVFQwcuRIp3tUV1dz9dVX/6q+N6VmDYSK\nojwKPApE1D60C3hJVdUvaq8rwIvAw0AQsAl4VFXVnObslxBCNBtFuejpydbA29ubkSNHMnLkSJ5/\n/nl+97vfMXv27AaBsKysjIEDB7JixYoG9wgODr6o1yorKwNg7dq1hIeHO137panZ5tTcm2WOAf8P\nGAgMAr4GViuK0rf2+kzgj8AjwGCgHPhSURQ56lkIIdygT58+lJeXN3j8mmuuIScnh5CQEKKiopy+\nAgMDCQwMpEuXLvzwww+O51gsFv73v/853dtgMHDkyJEG9+jatWuLvD9XmnVEqKrq/5330LO1o8Qh\niqLsBqYBL6uquhpAUZTJQAGQAKx0dU9FUQxA/V8dAly1E0II0bhTp04xYcIEHnzwQfr160dAQABb\nt25l/vz5jB07tkH7e++9l9dee42xY8fy0ksvccUVV3D48GHS0tKYOXMmV1xxBVOnTmXevHlERUXR\nu3dvFi9ezJkzZxxrigEBATz99NM89dRTWK1WbrjhBoqLi9m0aRNGo5EpU6a09LcBaME1QkVRtMAE\nwA/IAiKBUGCdvY2qqsWKonwPXEcjgRBIBmY3b2+FEKJ98/f3Z/DgwSxatIgDBw5gNpvp2rUrDz/8\nMLNmzWrQ3tfXl2+//ZZnnnmGxMRESktLCQ8PZ8SIERiNRgCeeeYZ8vPzmTx5MlqtlocffphRo0ah\n1Wod95kzZw7BwcHMmzePgwcPEhQUxDXXXOPyNVuKov6KvJNLegFFicUW+LyBMmCSqqqfK4oSh21N\nMExV1bx67T8BVFVVJzZyP1cjwmPFxcWO/xlCCNESKisrOXToEJGRkXh7y4rO+axWK9HR0dx1113M\nmTOnSe55oe95SUkJgYGBAIGqqpZc7D1bYkS4DxgABAJJwHJFUS55j6yqqlVAlf3v9iG3EEII9zp8\n+DBfffUVw4YNo6qqiiVLlnDo0CEmTZrk7q5dULNXllFVtVpV1f2qqv6oqmoysB14EsivbdLlvKd0\nqXdNCCFEG6HRaHjvvfe49tpruf7669m5cyfr1q0jOjra3V27IHfkEWqwTW0ewhbwRgDbABRFMWLb\nPfqWG/olhBDiMnTt2tVRsLstae48wnnAF8ARbGt5k4DhwChVVVVFUV4HnlMUJQdbYJwDnABa9twP\nIYQQHqu5R4QhwPuACSgGdmALgv+tvT4f2y7St7El1G8ERquqWtnM/RJCiCbT3JsORZ3m+F43dx7h\nQ79wXQX+XPslhBBtil6vB6CiogIfHx8398YzVFRUAHXf+6YgtUaFEOISabVagoKCKCwsBGy5drKT\nvXmoqkpFRQWFhYUEBQU55SZeLgmEQghxGeyH19qDoWheQUFBTX5gsARCIYS4DIqiYDKZCAkJwWw2\nu7s77Zper2/SkaCdBEIhhGgCWq22WT6kRfNr9oR6IYQQojWTQCiEEMKjSSAUQgjh0SQQCiGE8GgS\nCIUQQng02TUqGrBYLGRmZpKXl4fJZGLo0KGyG04I0W7JiFA4SUtLIyoqivj4eCZNmkR8fDxRUVGk\npaW5u2tCCNEsJBAKh7S0NJKSkoiNjSUrK4vS0lKysrKIjY0lKSlJgqEQol1S2nrV9NozDIuLi4sx\nGo3u7k6bZbFYiIqKIjY2lvT0dDSaut+RrFYrCQkJZGdnk5OTI9OkQohWqaSkhMDAQIBAVVVLLvZ5\nMiIUAGRmZpKbm8usWbOcgiDYTp1OTk7m0KFDZGZmuqmHQgjRPCQQCgDy8vIAiImJcXnd/ri9nRBC\ntBcSCAUAJpMJgOzsbJfX7Y/b2wkhRHsha4QCkDVCIUTbJ2uE4rJotVoWLlzImjVrSEhIcNo1mpCQ\nwJo1a1iwYIEEQSFEuyMJ9cIhMTGR1NRUpk+fTlxcnOPxyMhIUlNTSUxMdGPvhBCiecjUqGhAKssI\nIdqiS50alRGhaECr1TJ8+HB3d0MIIVqErBEKIYTwaBIIhRBCeDQJhEIIITyaBEIhhBAeTQKhEEII\njya7RkWTk/QLIURbIiNC0aTkYF8hRFsjgVA0GTnYVwjRFkllGdEkpGi3EMLdpOi2cCs52FcI0VZJ\nIBRNQg72FUK0VRIIRZO42IN9CwoKSElJISMjA4vF0mL9E0KIxsgaoWgSv7RGOHjwYLZt20ZNTY3j\n8YiICBYuXCjHOwkhmoSsEQq3utDBvoMHD2br1q0MGDBAdpMKIVodGRGKJpWWlsb06dPJzc11PKbT\n6RgwYADff/+97CYVQjSbSx0RSiAUTa5+ZZmCggKeeuopsrKyGDJkSIO2WVlZxMXFsWHDBjkDUQhx\nWeRgXtFq1D/YNyUlBZDdpEKI1kvWCEWzutjdpPZ2QgjR0mRqVDQrqTgjhGgpsmtUtEoX2k2akJDA\nmjVrWLBggQRBIYTbNGsgVBQlWVGUHxRFKVUUpVBRlHRFUa46r42iKMpLiqLkKYpyTlGUdYqi9GzO\nfomWlZiYSGpqKjt37iQuLg6j0UhcXBzZ2dmkpqZKHqEQwq2adWpUUZT/ACuBH7BtzHkFiAH6qKpa\nXtvmGSAZmAIcAuYAsbVtKi/iNWRqtI2QcwqFEM2pTaRPKIoSDBQCw1RV/VZRFAU4ASxUVXVBbZtA\noAB4QFXVlRdxTwmE7YgESyHEpWora4SBtX+erv0zEggF1tkbqKpaDHwPXOfqBoqiGBRFMdq/gIBm\n7K9oQXKorxDCHVosECqKogFeBzapqmrfSx9a+2fBec0L6l07XzJQXO/rWBN3VbiBHOorhHCXFpsa\nVRTlLeBW4AZVVY/VPhYHbALCVFXNq9f2E0BVVXWii/sYAEO9hwKAYzI12nZJioUQoim06qlRRVGW\nAGOAeHsQrJVf+2eX857Spd41J6qqVqmqWmL/AkqbvMOiRcmhvkIId2ru9AmlNgiOA25SVfXQeU0O\nYQt4I+o9xwgMBrKas2+i9bjQob4Wi4UzZ84AsH79ejnDUAjR5Jp7RPgGcB8wCShVFCW09ssHbHOf\n2NYNn1MU5U5FUWKB97HtJE1v5r6JVqKxMmz2zTO33347AC+//LJsnhFCNLnmDoSPYtspmgHk1fuq\nv/Y3H1gMvI0t39AfGH0xOYSifRg6dCgRERG88sorWK1WoG7zTExMDDfccAMRERFs3LhRNs8IIZqc\n1BoVrYI98I0ZM4aZM2dy77330rVrVzp06MDatWsdFWhk84wQojFtIqG+OUggbD9cHeobGRnJggUL\nnMqwyRmGQghXWvWuUSEuRmJiIvv37+e5554D4PPPPycnJ6dBLVI5w1AI0ZQkEIpWRavVMmKEbRNx\nhw4dXE59yhmGQoimJFOjotWRBHshxKWQqVHRbsgZhkKIlqRzdweEcMV+huH06dOJi4tzPB4ZGSln\nGAohmpRMjYpWTY5lEkJcrEudGpURoWjVtFqtpEgIIZqVrBEKIYTwaBIIhRBCeDQJhEIIITyaBEIh\nhBAeTQKhEEIIjyaBUAghhEeTQCiEEMKjSSAUQgjh0SQQCiGE8GgSCIUQQng0CYRCCCE8mgRCIYQQ\nHk2KbgshhHCb1nDCjIwIhRBCuEVaWhpRUVHEx8czadIk4uPjiYqKIi0trUX7IYFQCCFEi0tLSyMp\nKYnY2FiysrIoLS0lKyuL2NhYkpKSWjQYysG8QtA6pmeEaE/sP1PHjx+nqKiI4OBgwsPDGTp0KABR\nUVHExsaS/tG7aHZ/BmdyYdRcrFYrCQkJZGdnk5OT86t+DuVgXiEuUVpaGtOnTyc3N9fxWEREBAsX\nLiQxMdF9HROijXL1M2UXERHBQw9Mpr/hGO/eGo3mr73BagZFA3FT0QSEkpycTFxcHJmZmS1yMLdM\njQqP1pqmZ4RoD+w/U507d0ZRFG699Vb++c9/cuutoxkUpuEvw6w8Uv5X0u/2pWPBJlsQDI2FW14G\nnTcAMTExAOTl5bVIn2VqVHgsi8VSNz2Tno5GU/d74eVMzwjhqew/UzExMWRnZ9t+tpYvQZOdirp9\nJcrJfY62eaVWrDETCL/1TxAa43SfrKws4uLi2LBhw68aEV7q1KiMCIXHyszMJDc3l1mzZjkFQQCN\nRkNycjKHDh0iMzPTTT0Uom2x/0zdMeomhhqP88HIYjR/6wfrX0Q5uQ+rxouPdppZHTCZrovKuOtf\n+7CG9HG6h9VqZd68eURGRjrWE5ubBELhsezTLvZpmPO19PSMEG2a1YJ1/9csT/Dmt6de4f1xPgSe\n2gao0O0GuHMJ5Y9t5960c5zqcDUW1TbyS0hIcFqWSEhIYM2aNSxYsKDFZmJks4zwWCaTCYDs7GyG\nDBnS4Hp2drZTOyGEC0X7YHsK7PiEm0qOQ38vUKvJOWXBe/Bv6Xr709ChGwDZWVkAnDt3DoAXXniB\nZcuWERcX57hdZGQkqampLbpRTdYIhceSNUIhLlH5ScheZQuAJ35yPKx6B7JieyXb1GhWfX+Y2Nh+\njp8t+8/Uzp07iYmJYdeuXeTk5AA0WeqSpE8I8StptVoWLlxIUlISCQkJJCcnOxb5582bx5o1a0hN\nTZUgKNoti8VCRkYG69evZ+vWrfj6+nLjjTfyxBNP4OXl5dy4pgp+/hK2r4ScL8FaY3tco4OokTDg\nHpSeo/D9v8/5a1ISAwcOZM2aNYwZM4aEhATS09P54osvGDRoEGvXrnX62WqJFIkLkRGh8Hiucp4i\nIyNZsGCB5BGKdqN+0YiQkBC++eYbFi5cSEVFRYO2Go2G6dOnM//VV+HYVtvIL3sVVJ6ta2QaAP3v\ngZjx4B/s9PwL5RE258/WpY4IJRAKwYUry0jVGdHWXSgw2fn4+BAeHs7+/fu5MlDhvn56pg3vQrCm\nuK5RgAn6TYT+d0NI9AVf80KVZZrr50cCoQRC0Qyk6oxo6+wJ7t27d+fAgQMu21x11VWcOLSPpD56\nHru+I4M6n3NcU/W+KNF32IJf5DDQtN5fAiWPUIgmJlVnRFtlX/tbsWIFjz/+OJGRkY4gqCiKo51G\ngZHdtTzf5wiFMwN5d6wPgzqfw6qqrD9Yw5T0c7zpOw0S34YeN7XqIHg5ZEQohAuyo1S0VampqUyd\nOpX8/HynxxVFwf553zdYw+T+eu6N1RNurPu3ve+kheXbzXy4w8zRElvbJ554gsWLF7fcG7gMsmtU\niCZkr5CRkpLSaNWZliwKLMTFmDlzJq+99prLa5194J5YLyb30zMwrO6Xt1MVVlbuqmH5tmp+OGFt\n8LwePXo0W39bCwmEQrggVWdEW2HflPLZZ5/x97//nZCQED777DP69evHiuXvsO7Np5ncT8+tPXXo\nNLZp0WoL/OeAlWU/VfJ5Tg3VFud7GgwGzGYziqLw2GOPueFdtSwJhEK4IFVnRFvgajOXVquh+kAm\n/qc+4felafxhgq/j2pbjFlbsspKtXsXXm7c1et+qqioAZsyY0TCfsB2SNUIhXJA1QtHa2TdzjRkz\nhtGjR/Pas09wfz89j1wXRJh3paNdha4Dr2/I54MdZvaerJv69Pf3p6yszOW9HXmE8+c3+/toSq1y\n16iiKDcqivJ/iqKcUBRFVRQl4bzriqIoLymKkqcoyjlFUdYpitKzOfskxMWwV51Zs2ZNqygKLER9\nFouF6dOnc9edo0n/cyJJpe9w6MkAXor3Jsy7knMWDan7DVjuS8d31kG+87nZKQgCDYKgVqvl2muv\nZeHChZw7d67NBcHL0dzpE37AduDxRq7PBP4IPAIMBsqBLxVF8W7mfgnxixITE0lNTWXnzp3ExcVh\nNBqJi4sjOzu7xYsCC89gT3tISUkhIyMDi8XiolENu9Jf55VrCvjomh/QrHmSkHP7saoq3x7T8Nfc\n3mwbmcqEFUVkHlOgdnNXYwwGAy+88AJVVVVs2bKFP/3pTx4xHVpfs64Rqqr6BfAFOOeu1P5dAaYB\nL6uqurr2sclAAZAArGzOvglxMRITExk7dqxUlhHNprq6msWLF/PJJ5+wa9cuysvLHdecijfk77TV\n+dzxCf3KC+kXqwdrNXS+Cmu/u7n+0b9R4xvCjz/+wOjChQAcPHgQvV7P+PHjAUhJSeH48eNs3LgR\nf39/7r//fkaMGOHx/55bbI1QURQVGKeqanrt37sDB4CrVVXdVq/dN8A2VVWfbOQ+BsBQ76EA4Jis\nEQoh2hKLxcKkSZP49NNPOf9zOCAggFmzZrF7SwbB+Rt4/s6eBFUedVyv1htZurGI+Gn/JPaW+0BR\nHGuGAwcO5NChQ5w6dcrpnjNmzGj3052tco3wF4TW/llw3uMF9a65kgwU1/s61vRdE0KI5pOWlkZA\nQACffPKJUxD08/Ojc5A/t115jn47X2J5v80svMWboMqjqFoviL4T7k5BO+NnFu018ewbn2Ktfb59\nKr+oqMgpCIaGhvLJJ5+0+yB4OdpiibV5QGC9ryvc2x0hhLh4aWlpjB8/3nE4be/evVGAKfG9WHST\nhf1/UFiZ5MttPfUoqJQE9uaRNefYNORdmPgB9L4NrZePy81cJpOJ2NhYFEVh2rRpbNiwgWPHjjFh\nwgT3vulWzp15hPb6P12A+lnJXYBGE1xUVa0Cqux/P3/tUQghWiuLxcKf/vQndDodNTU19OigsHTi\nlXQ9c4zuHfIB2yaVIyWwfFsVH2w38/oHs/nHj7cz7KTzLk/7CHD69OluP+G9rXNnIDyELRiOoDbw\n1eYEDgbecmO/hBDisjR2dFdmZibFBYf5bT89U/r7cv2VOmAzdNBQo/XldOgNTJj7GZmHLdgnTH/4\n4QfAdfEG2czVNJp1s4yiKP5AVO1ffwL+BGwATquqekRRlGeA/wdMwRYY5wD9gD6qqla6uKWr15CE\neiFEi2ss2Lmq9hIV2Y3lf57MFacyCTnzI94620yWxapSEjyIFz7by/GAq3nn/Y8ICgpyep24uDjy\n8vKkeMNFaK1FtwdhC3x2f639cznwADAfW67h20AQsBEYfbFBUAgh3OHTTz9l6tSpFBTU7fWLiIhg\nwoQJLFiwgDFjxpDy0Uf0C1Eo/vYtfPavJSj3b7aGOoWdBRY+2FnDRztrGDDUlwce+zt33XUXBWPG\nOL2Ot7c3WVlZpKamShBsRlJiTQghGnH+qC8uLo7Jkyfz8ccfO7UzmUyEhYXx448/Mvr6Aaz9yxQ0\nOz6Gwt2ONmfNetJytKTm6MjYdwaTycTBgwcBiI2N5bbbbmPhwoXU1NQ4nhMSEsJbb73lEet9qqpe\n9p4POaFeAqEQogm5muKsf6YfQJcuXfjjH//I/77fhOHQV0zu58XN3XVo7fvxtQbofRv0v4fNRb5c\nd8ONvPjii8yePRuA8PBw8vLysFqdy595eXnx7LPP8uyzz7bLkWBFdQ1780vZk1fC7hMl7Mkr4WRZ\nNd/MGH5ZwbC1To0KIUSbY09x6NOnD0lJSRQUFJCZmekUBBWgt/dJQr+fw4cD/PC+uu6Uh2JjNIHD\nHoE+CeBjW/PrayoFoGfPnqxatYpHHnmE48ePO72uwWBg5syZzJ49u10EQFVVKSytYveJEnbn2b72\n5JVw6GQ5rsZgRaVVhBhbvsKmBEIhhMerPwUaEhLCb3/7WxRFYffu3ezevdupba9OGn4/2Mh9/XR0\nMVTXPmqmqMaPNzae5sMd1cx540nuGXiP0/PqH901fPhwxo4dS0ZGBhkZGQAMHz6c4cOHt9kAaLZY\nOVhUzu68YvbklTqC3+nyapftgwMMRJuM9DEZ6RNmpI8pgE7+Bpdtm5sEQiHcoLEdh6LluZoCBdsx\nRQsWLOC5557DUnaSiTF6pvTXM+QKHWAFqim3aPnsZ4WlWcXkVBkoLLSlOHfp0sXpXlarlXnz5hEZ\nGcnQoUMB22kPI0aMYMSIES3wLptW8Tkze2pHd7tPlLAnv4Sf88uotjQ84V6rUeje2Y8+YUZH4Is2\nGQkOcE/Qc0UCoRAtzNUHr1NxZdGs6v8SkpOTwwsvvGDb5ZmSQq9evQgODsZqtVJVUUZg/iaWxpdz\nx1UBeGlta1c1VpUv99cQeOMf0MeM5f6X42vvXOR4jUWLFuHj40NMTAzZ2dnMmzePNWvWtLndn6qq\ncuzMOXY65WJlAAAgAElEQVTVruPtrg18x8+ec9ne36Aj2hTgGOVFm4z06hKAt751v2cJhEK0oPqH\nqaakpDg+KF955RWSkpKkIkgzsVgs/Pe//2XmzJns3r3b6Xgjb29vJk+ezJAhQ3h90SKu7gJTBnhz\n/wAfgtTV0EcPwE95FlL363nn+2IKylWuy/2el166vcFrTZ8+nVWrVrW5ai+VZgs/F9TfwGL779Kq\nGpftw4N8nEZ5fUxGrujgg0bT9qp9ya5RIVqInHrfsuxn+73xxhusXr26wc5MO51ORxcfC/+Z/1s6\nHfsvJl2x41peqZUVO838+7AfmT+fZubMmcyfPx9FUejSpQv5+flO95o4cSIrV65s9VPfJ8uqHLs1\n7aO8gyfLsVgbxgMvrYaeXfwdU5p9woxEhxoJ9NW7oecXJukTEghFK5eRkUF8fDxZWVkMGTKkwfWs\nrCzi4uLYsGEDw4cPb/kOtgP2ALR69WqWLVtGcXGxy3Z+fn5QXc64aD2T++kZ0V2Lpnbb/jmzSrYl\ngudTd7PuoIWAwCAiIiLYtm0bWq0Wi8WCr68vFRUVjvvpdDo6duzIiRMnWlXAs1hVDp0sdwQ7e+Ar\nKq1y2b6jn1dtwAtwjPZ6BPuj17aN8xkkfUKIVi4vz1ZbPiYmxuV1++P2duLXSU1NZerUqQ1GafVp\nFBgeoWVyPwtJfY341RvUfJNbQ8f4xxj++Ovo/Qt5/PHZfPnnP3P27Fm2bbOdA2CfUq0fBE0mE3l5\nebz11ltuDYJlVTXszas3yssrZV9+CZXmhiNhRYHITn5Eh9VNa/YJMxISYPDIgwwkEArRQuxFk7Oz\ns12OCOtvrxcXZh/5HT9+nKKiItatW8fatWsbbd+7s4b7++l5aJAfXXzq1rxyTll4f4eZD3eYyT2r\n8lxXbx567Clee+01Xn75ZTp16sQbb7zB22+/TUZGhsvpVW9vb1atWtVi63+qqnKiuJI9J+ry8nbn\nlXD4VIXL9j56Lb1rN7DYpzZ7hwbg6yUf/3YyNSpEC5E1wqZxMSM/gE4+CnfH6Jjc34vfhNd9P8+c\nU/l4l5n3t5vJOmZxes6MGTMYN24c999/PwcOHHC6ptPpGD9+PHfccQdFRUUEBwcTHh7erOt/1TVW\ncgpLHXl59qBXfM7ssn2o0bt2SjOAPqZAok0BdOvkh7YNbmC5FLJGKIFQtAH1d40mJye73F7fmncW\nukP90d8777zDhg0bnK7XL3sWFODL8PBqpvTXc3tPHfp6KQ+f59Sw+Vwkf12zmypLg5dxEhkZybx5\n88jLy+PAgQP06NGDxx57DC8vr2Z5jwBnyqvrTWva1vQOFJVhtjT8jNZpFKJCztvAYjLS0a/5+tcW\nSCCUQCjaCFd5hJGRkSxYsECCYD0Wi4U5c+awaNEiSkoafqb5+PgwadIk3nnnHX4TrmVyfz33xHjR\n0aeuzU8FKu/9VEXKTjNFFRf+rAsKCmLJkiXNPsqzWlUOn65wqrO5O6+EvGLXh+4YvXVOwa6PyUjP\nLv4YdDJrcD4JhBIIRRvS2rfXu1tqaipTpkxx2pRiFx4ezvHjx7npmiiG+B1hcj89V3Wu+96dKIMP\ntttOdz9a5e8yiLrSHOt856ot7M2vt5Z3ooS9+aVUVLsekl7Z0fe8UV4A4UE+HrmB5VLIrlEh2hCt\nVispEuex5/0999xzbN682WUbfy+4uXMhk0f4clNkIWAr0FxhhrQ91Szfbuaq0Q/zxrqltc/45c/C\npjjqyFFc+rxRXu7Jclyk5mHQabgqtOEGlgDv1peb5wkkEAoh3C4tLY1HHnmEoqIip8e9vb2prqrk\nnsEmRpvOkBitx1dfNzo6qovk+dTdbKvqyva9tWf7GWxzo+fn+tUXEhLCAw88wC233PKrC13bi0vb\ng519pHeqkeLSnf29HMHOnqoQ2dkPXRvJzfMEMjUqhHAbi8XC3LlzHefz6XQ6x8G0fYI1TO6v575Y\nPeHGuqCx72RdysPzC97m4YcfvuBrBAYGUlxczLXXXsvEiROZOnXqRW96Kak0N0hT+LmgjOqahmkU\nGgUiO/vRJyywXq3NAEICWv5YIU8lU6NCtHPtZV2xfvWXlStXOqVBdPCycM9ALyb30zMwrF7KQyWk\n79eydHMJW47Xra8tWrQIuPDor2PHjrz77rsXnPq0F5c+f2rz2JnGi0v3Dg1w2sDSq0sAPl5t7/+H\nkEAoRJvQlk6sqB+wO3XqxM6dO8nNzaVHjx6EhYXxzDPPOL0PgxbG9NLx3B1R9PU67kh5MFsg65SR\n1zfkszanhukzk9myap7Ta9nPCrQHwdDQUCZOnEhERESjeX6VZgs5BWVOqQp78koorWy8uHT0eQnp\nXTv4tsni0sI1CYRCtHJt6cSKlJQUfv/731NWVtZom+7du2MyhXLPDT24pw/0OLeDDj4KcAJQ2HLc\nwvvbq9mj7cvXm7c5nvfJJ5+4vJ+fnx8PPfQQ48aNaxD0TpZV8d3B0451vN15JRwocl1cWq9V6BkS\ncN65eQEE+Xp2bp4nkDVCIVqxtlKNxmKxEB0dTU5OzgXbRXXSc3dfhcn99PTsVNffo8VW0g8aSD/k\nzbd7ChzrhN26dePw4cMu7xUcHMwTTzzBs88+C4qGQycbbmApbKS4dAdfvdMhsX3CbMWlvXSygaUt\nkzVCIdqhzMxMcnNzSUlJcQqCABqNhuTkZOLi4sjMzGyxdIz6lV4KCgr47rvvWLVqVaPtjQa4K8bA\nfTFahkXUfeSoel/U6Du577U1pG87xbnKMmbPns3XL74I2N6fqyA4+o6xJP52KoYu3dmTX8b4pZvZ\ne4Hi0hGd/BqcqBBq9JbcPOEggVCIVqy1nVjx6aefMnXqVAoKCi7YTqvAzd21TO7vxbjeOnxqUx6s\nqso3R+Dd/1Viir+P+c++QRJppIwfD9imVh988EHeffddrFYr2oBOeIV0xyukO77hvQjrcy17qzXM\n/b4S2O30mj56rS03r97UZu/QAPwM8jEnLkz+hQjRirnrxApXO1STk5N57bXXLvi82BBbysO9sXpM\nAXUj2D1FFpZvN7Nip5ljJbblmIiSz/nLX60kJiYyZ+48Xnr9bY7rr+CzXA0hE+fi1SUSrY/zcsfJ\n2pnOLkZDgzqbER5UXFo0LVkjFKIVc8caoasdqsHBwQ2S3e26+CncE6tnSn89A0Lr+nC6UuHD7ZW8\nv93Mj3l105Yab3/0IZEYQrrTfeCNBFwZzZGzZtA07L9WoxAV7N/gRIVO/oYmea+ifZE1QiHaIa1W\ny8KFC0lKSiIhIaHREyuaMgiev0N1+/btxMfHO7Xz1sHsu68l1pLNqCgdutqRWLVF5f/21fD+DjOZ\neT6UajriFRJJYFR3DCHd0YdEoDOGOO5zBjhTYgWNFmtlGZqSPIbGRnLH0GvoYzISFeKPt15y80Tz\nkkAoRCuXmJhIamoq06dPJy4uzvF4ZGRkk6RO2Gt8rlu3jsWLFxMWFkZMTAzl5eX4+PhQWVmJ2Ww7\n/+6GK22nPNzVR0+g9z7AVhsz65iVjw4F8e+SXpR16IV+UHf8gyMwGnxdvqb5bD7mgoNUFx3Cr/oM\nY64fwN1jb+HGG5PaZJEA0bbJ1KgQbcTFVJa5UBtX11JTU3nooYcoLy93+ZrBwcHcdfO1dD6xnsn9\n9XTvUDc1e5IOvJ8XwWeaERwJ6I/iYmpTrammuugw1YUHMRceorrwEDcP6s3TTz7R5ivkiNZHjmGS\nQCg83IWqzwANrnl5eVFd7bpQdKABJvT14oGBflwfVlfSrEw18LllCGnWoXxv7Y1KXWC0lJ+hujbY\nVRccxFx0CPOpY6Da1gf9/Px45513mDhxYhO+ayHqyBqhEB6sseozc+fOJSkpCVVV6dOnD8OGDWPr\n1q2Ul5c7BUHFyxffLhHc2jeQSVHl3NrxGN5KDWDBoipsssawyjKUr6yDqLDqMZ8+jrkwE9/q0xT+\n/BPVhQexlp9t0C9FUUgcP55HH330V5/yIERLkRGhEG3c+TtLVVV1FLX+6KOPKCwsdGqvNYbgFRKJ\nV5fueAVHMCDcm4lBu7lTu4lgpe6X6H3WK0itGsKn+SaO5p2iuvAQNUW5VBXmota4rthi5+3tzdNP\nP80LL7wgwU+0GJkalUAoPFRGRgbx8fFkZWVx4sSJuilQrQ6vzt3Qh0TWJqVH4hUSicbbnxDOMFa7\niURtJtGao457nSWAlCOd+XC3hi37TlBzthCo+4xobDpVURSGDRvG7373O5eFroVoCTI1KoSHysvL\nQ+Nj5O1/f8uq9ZsJG/4HrtAHoQk0oWjrfsS9qeIWzVYSNd8yVJuNVrEFuCoL/DtHZfn/zjFj6SoO\nf7WO75a92uB1Jk6cyObNm53KngUEBDB+/Hj+8Y9/XPQZf0K0NhIIhWhDLFaVw6fKHefm7TpezKZd\nGrr+8SO+tkKH+L6cA+xjMeu5s1xTkkmSYTN3djqMUV+X2L7piC3f75NdZs5W2h67J6+gQZUavV5P\nWFgYK1asAGgXZyIKUZ9MjQrRSpVX1bA3v5Q9ebaAt+Xn4xwpsWBWXZcRM58+gbnwIEGU4X04k4ld\nDnNv7xoigup2dpZoOrBT05cpf/2CA2ca/ux369aNo0ePYrU6F7BetWpVqznqSYjGyNSoEG2UqqoU\nlFSxO6+YPXmljnPzck+V0/D3VAWruRJzUa4jTUFfXsjYYdey+sOlTOrvw6S+CnED7D/aGkqqVD7Z\nZWb5djOHrQGMT+rFgTOfu+zL+ac9GI1Gli1bJkFQtGsyIhSiBZktVg4UlbH7REndCeknSjhTYXbZ\n3lJ2iuqC2ty8woO2nZtnTjhy83QauDVKx1M3hRLXqRiDzjZatFhVvjxQw/vbzaw9ABXV1gajvAvx\n8fFh5syZPP/88zL1KdoMGREK0coUV5jrDomt/TOnoIxqS8OApNUo9Aj2Izo0AEPlSbZ8+RmZ/5eC\ntaLYqZ1GowHVyjUm2ykPk2L0BPtpgBJA4VCFH9uVPjz65tfklzX8JfdCSfQAnTt3ZurUqTz77LMS\nAIXHkEAoxGWyWlWOnTnH7rxidtdObe7JK+H42XMu2wcYdETXHhTbNyyQXl38KMzZzhdr0lk5fyX5\n+flO7f39/SkrK+OPDyQxQPszv/E6QN+Qeqc8VOtY9mM5+/2u5R/pm7j55hiXQRBoEAR1Oh3+/v4s\nWbJE0h6Ex5JAKMSvUGm2sK92A4t9lLcnr5SyqhqX7a/o4ON0bl4fkxGT0YuNGzeSl5fL3i9ymPHO\nOxw5csTxHHvg69ChA1VlZxgbWcXk/r6M7PZfFFRAS2WNyqZTHViw7gQ/14RxMPc0t9ziS2pqKr//\n/e8v+B68vb257bbbeOyxxzAYDAwdOpTw8PAWO+FeiNZGAqEQjSgqrXKs4dkD38GiMqwuBlteOg29\nuvjTp/Zk9GiTkd4mI4E+ttMZ7Cc8/OmFpXz11VeUlDS+fHGuooL4CC1Pjwzgxs41+HvZd4mqcGUc\nP1p6MWLq36mwVmI21wC2IBoQEEB8fLzTWmBiYiJpaWnMnDkTLy8vhg8f7lTqrLS0FGi5E+6FaI0k\nEAqPV2OxcuhkbW5e7Qhv94kSTpa5LiPWyc/LcSq6Peh1D/ZDr9U0aFtdXc0f/vAHVq5cSWVl5QX7\nMe2+23h65BWw42PC/VXgLKBwpEzHv34o4+M9sPPYajYtXUpxFYBtg42iKKiqyqpVq1i1ahVQt9uz\nY8eOpKWlMW7cuBY94V6ItqRV7BpVFOVxYAYQCmwHpqqquuUinyu7RsVFK600O3Lz7GkK+/JLqapp\nuIFFUSCys59tlFcb+PqajAQHGFAU17l8dtXV1YwePZoNGza4vG4PXh19FP4Q14mxEecYfEXd2tzZ\nSpWPd5lJzdGzbl+p43G9Xo+iKI1ueNFoNEyYMIEVK1ag1WrdcsK9EO7SZneNKooyEfgr8AjwPTAN\n+FJRlKtUVS284JOFaISqqhw/e84xurNPbR45XeGyva+X1rGBpY8pkGhTAFeFBuDrdXE/Ivapz4yM\nDNavX8/mzZtp7JdMvQbGxfhyd28rt/fS4aWtBrTUWFW+PqIl/qm36T70Ac6UVuLl5Ryg7Qfk1hcT\nE8Pw4cPp0aMHjz32mFOps5Y+4V6ItsjtI0JFUb4HflBV9Ynav2uAo8BiVVX/4qK9ATDUeygAOCYj\nQs9VVWMhp6DMKS9vT14JJZWuN7CYAr2dRnl9TEau7OiLRnPhUV5j0tLSeOSRRygqKmq0jZeXF/07\n1zC5vxf3xOjo5Fs3MvtfnoXUHB3vbCmhsFxlw4YNZGRk8OKLLzZ6P71ez0033cRnn32Gj4/PRfXx\n/PMIIyMjWbBggSTLi3ajTY4IFUXxAgYC8+yPqapqVRRlHXBdI09LBma3QPdEK3S6vNoxrWkPfPsL\ny6hxsYNFp1GICvF37Na0r+d18Lv84tD2095Xr17N66+/DtiC0/kjtq5Ghfv66ZkW34UQTV1O4IlS\nKyt2mnl/u5nsQiuJiYkUlqcBttMk/vOf/7i856Xm+SUmJjJ27FipEyqEC24dESqKEgYcB+JUVc2q\n9/h8YJiqqoNdPEdGhB7AalXJPVVum9rMK64NfKXkl7jecBLoo2+QphAV4o+XruEGlsvlanQFEB4e\nzvHjxxlyTQy9LPuY3E9PfKQWTe16YoVZ5bM9Nezzvpq5K79z2n3avXt3Dh482OC1IiIi+O1vf0vP\nnj0leAnxC9rkiPBSqKpaBTi28/3SpgXR+lVU24pL1x/l7csvpaLa4rJ9RCdfpx2bfcKMmAK9m/zf\nQv11P6vVSseOHTly5AiLFy/m9ttvZ8aMGTz++OMAdOoYRLQhn3kJ3kyIPYa3pm66csOhGjaWdeWd\nTYUcLihFp9vSIAWjfhA0GAxUVVXx4osvSoUXIVqAuwPhScACdDnv8S5AfsPmoi2zF5c+fy3vkMvi\n0mDQaehtMtLHFOCUm+dvaP5/tqmpqTz88MOcPXvW5fUdO3YQFRVFdGdbqbNpNwXhXW3f2GLl51MW\n3t9u5sMdZg4Xq8A+QkNDAaipcb12aRcWFiZrd0K0ILcGQlVVqxVF+REYAaSDY7PMCGCJO/smLo+r\n4tJ78ko5Xe56239wgOG8qc0AIjr5oXORm9ccqqurefPNNzlw4ADZ2dlkZGQ02razr8L4K4qYVP5P\nFj3uX3uD05h1/vwz6zTph7zJOmqlrMz5vZ5fOq2+oKAgHnjgAcaOHSvTn0K0MHePCMGWOrFcUZSt\nwBZs6RN+wDK39kpctOJzZqe8vF8qLt29s1+DhPTgAIOLOzc/i8XCvffey6effnrB0xkMOoXbe2qZ\n3E/PbT116LW2aVizReXznBq2VPdg5j/X8/jzXYBK+zrFBQ0ePJiRI0c2qPYihGhZbg+Eqqp+rChK\nMPAStoT6bcBoVVUL3NszcT5VVTl6+uKLS/sbdLV5eXWpCr26BOCtd/8HvsViYe7cufzlL3/h3DnX\n/QcYHK7l4Wv9GNfLlvxu98NxC+uKOvPX/x7hZIUKZLPmwEgAhg0bxjfffNPoPTt06MDbb79NUlJS\nk72f5mbfJSs7TkV75PY8wssllWWaR6XZws8FzhtY9uaVUtpIcenwIB+nUV4fk5ErOvhccm5eU7N/\nkB8/fpz//ve/fPrpp1RU1CXXe3t7A1BZWUm3QFvKw+T+enp1qvuwP1ai8uGOaj7YYWZ3kZV7772X\nFStWMGTIEDZv3tzgNe2VXUaNGsXVV1+NTqdrk6M/V7tkIyIiWLhwoaxjilbFY3aNiqZXVNpwA8uB\nxopLazX0rC0ubV/Piw41Euirb/mO/wJ78Fu1ahXvvfceZWVlLttpNBqiu1/BAK/DTO7vy/CIuh+L\n8mqVnZZIjEN/T+ydTzh9T+wJ9L1798ZsNvPjjz863bdr165tPlikpaWRlJTEmDFjSElJcVSleeWV\nV0hKSiI1NbVNvz8hQEaEHsViVTl0soxdtTl59vW8olLXxaU7+nk1mNrsEezvsrh0a/Pxxx/z6KOP\ncubMmUbbaBS4ubuWyf28GBetw1dvG71aVZUNhyx8dsiH974/hSGgIwcPHiQoKKjRe0VGRjJ//nw6\nd+7cbqYPpU6paGtkRCiclFXVsLfemXm7T5Swr6CUSnMjxaU7+RF9XgWWLsZfLi7dGiUkJLB69epG\nrw8w6ZkUo+G+fl6Y/Ove396TFpZvN/PRzhqOFFsB29Rp+enTXHedc6Ej+0nvRqORtLS0NjfdeTEy\nMzPJzc0lJSXFKQiCbRSdnJxMXFwcmZmZcpahaNMkELZxqqpyoriSPfV2bO7OK+HwKdfFpX30Wnqf\nN8rr/SuKS7c29ZPeAfbu3cvq1asd63MAN998Mzuy1nNPjJ4p/fVcbaoLWCcrrHy6F979sYKtJxrf\nNbpnz54GjymKwrJlyxgxYkTTvqlWwn5GYUxMjMvr9sflLEPR1rXNTz8PVV1jJaew1HGiwu68Yvbk\nlVJ8ruGJBAChRu/aYFd3okK3Tn5oW8kGlktVv87ne++95zLpXafTocPCnVfpWHhzGabr/NHVvu9q\ni8qan2t4f7uZz3NqcDFI/kXh4eHtPundfkZhdna2nGUo2jVZI2ylztiLS9fbxHKgqAyz5QLFpeuN\n8qJNRjo2QXHp1qD+qG/37t1kZGRw+vRpx3WDwUBycjKFhYW8+eabxHXVMqW/nrv66gnyrgv6m4/Z\ngt/Hu2o4fe7X/bv39fVlwoQJjBw5kvDw8Da//ncxZI1QtDWXukYogdDNrFaVw6crGpyokFfsuri0\n0VvXIBm9Zxd/DLr2+UH0S6XOvLy8iI2N5czBn3gxMZrrfA/To2PdB/axUviurCv/2FzCdz8X/eIp\n8XYhISEMHTqU6OjoNpny0FTq7xpt7CzD9jwqFm2LBMI2EAjPVVvYm++8gWXvBYpLX9nR1xHsok0B\n9A0PJKwZiku3FvVz/YqKili3bh1r165ttL3RAE+O6MrIkFMM7VY3y19arZK6y8z7O8xsLfKmrLyc\nQYMGsXXr1gu+fnx8PA899JDHjPgulpxlKNoKCYStKBCqqkphaZVTXt7uvBJyT5a7zM0z6DT0Dg2o\ny8ur3cAS4N36cvOaS2NHG9nZN79cf91gjEU/Mrm/noTeOrx19pQHWHewhg931vDfo97kn677GTAa\njVit1kbzCIOCgnj77beZMGFCk7+v9kIqy4i2QNIn3MRssXKwqLxeYWlb8DvVSHHpzv6GehtYbNOb\nkZ1brrh0a2GxWFi/fj0ffPABP//8M1u2bGn0g9XLy4voDjVM7m9g2oiTaCp8Hdd2FVo4G3Ebd72c\nyolSlTFjxpC/fY3T80tKGv489OzZk4kTJ3r0tOevodVqJUVCtFsSCH+FkkpzgzSFnwvKqK5puO1Q\no0CPYH+nUV60KYCQAG839Lz1sNf4nDdvXoP1OovFgre3N3PmzGHGjBmE+itMitUzuZ+e/qG137eK\nIsx6I29uLGL5NjM/5Vt5+ulunCi1DbWfeuop9uzZw4EDB1y+focOHXjrrbeYOHFis75PIUTbIVOj\nLqiqyrEz5xpMbR47c+Hi0o46m2Gtp7i0O52/5pebm8vy5cudNr4EBARQWlrq+LuPDl685zfEWLK5\npYfWkepRVaPy7301/FDdg2eXbSCoU7DjOZ07d+bkyZN07NiR6OhoNm3aBMADDzyAr68viqIwePBg\nunbtKlN6QrRjMjV6GY6dqeC7/afqAl9+CaWVjReXrn9mXh9TYKsqLt1a/NKan9FoZO3atTz44IOU\nlZZyw5VapgzQM6GPHqNhL/Z/mt8dtac8mDlbCYqyi+/uTABs1U2sVisnT54E4PTp02zatImQkBDe\neust2cghhLgoEgiBrAOnmLlqh9Nj9uLS9dMU+phaZ3Hp1sa+5X7gwIEcPnyYUaNG8cMPP2CxWByj\nwZKSEsqP7mDpxCuJLMkjskPdGmnuWSuf7tPyXcWVpH+70+neYWFhfPfddwBO5wded911jBgxQtb8\nhBC/mkyNAvvyS3lpzS6iQ2tHemFtp7h0a1F/GvTpp59m4MCB7Nq1i9jYWKZNm8aIESN44YUXeP0v\nLzCxr56pNwbT11i3i7OkSuXf+zW8vaWMjUcsjBo9mi+//BKDwXDB3D/Z8SmEsJOp0ctwVWgAK37X\nsISUaKj+NvqQkBAA1qxZw8qVK8nPz3e0++677zhz5gwpKSnkHtzPmF46ft9hE/9vegAGnQKUYbGq\nnOl0DbNW7uDDrWewag1UVdlyKgMCAujfvz/btm1z2Y/g4GCeeOIJnn32WRn9CSEui4wIxS+qX9vz\n/IBXX2hoKImJibz55pv0798PTUE2Hz1zB93L/4eXudjRbkeBhZV7NSzbWsZf315Bbm4us2bNumAf\nbrvtNkaOHElwcLAkvAshXJIRoWgWrja9dOjQgTNnznD99dezdetWqqqqmDt3Lps3byb9g7eYEefF\nM6Mr6GTxh7MbADh5TmFjSRgnw2/m4aWLHffatGmTI4fQflpEfTL1KYRobhIIRaPsm15uu+02SkpK\n6Nq1Kw8++CDPPPMMACNHjmTTpk3cfGMcRevfIH3yILg6AI0CWAqprIH0vWai736Rg0oE4ydMJCTk\nNEFBQZSWlmKxWHjjjTccr6fRaEhKSqJXr14AsvGlHZIKNaI1kqlR4XD++t/vfvc7OnfuzLFjxxpM\nh0b3voqrDCe5s1s59w8yorPUbWjJPFzDBztq2KvtTeYPO+jfvz+PPfYYS5cu5aeffnK6j5+fH337\n9uWuu+5i6tSpeHm1jxMzREOuZhciIiJYuHChpLqIJnGpU6OoqtqmvwAjoBYXF6vi0q1atUqNiIhQ\nAacvRVHUa665RgXUvLw8dckLT6ov32RQc5/0V9XZRsfX/qn+6o6/T1JLj+xQAXXatGmqyWRqcD9A\nNZlM6rRp09QNGzaoNTU17n7rogWsWrVKVRRFveOOO9SsrCy1tLRUzcrKUu+44w5VURR11apV7u6i\naF9++jIAABWWSURBVAeKi4vtnzNG9VfEERkRCqe8P1ejvwfvScCQs5ZXJ11NQPFex+NnK1X+fUDD\n5opI3vp8Gxs2bMBgMBAXF8f69etZtGgRW7duZf78+Zw6dUo2ungoOddQtBTZLCMuSv1DbgGGDh3K\n9OnTGThwIFu3buWOO+5g1qxZZGZ8zab3ZvPgQB9Gd1+PVy8fKN6Lqmj5+oiGL04E8sa6g0y8dwrL\nU5fj6+uLRqNhzpw5mEwmFi1axNq1a+W8OkFmZia5ubmkpKQ4BUGwrQsnJycTFxdHZmamFPYWbiGB\n0IOkpaXx6KOPUlhY2OBacXExd9wxhvQ3/oxm58cMqF7BjLvrTnnYXgjv/VTJKVM8Vw0cysJ3nwMg\nKysLRVHw8vJi2LBhjva7du2SICgAyMvLAyAmJsbldfvj9nZCtDQJhB7CPv2pqio33HADL774IgCP\nPfYY5Xk/c1/fcqYP3IXmnREAeAN5pVY+3Gnm/e1msgutzJo1i4+XLeODz7503Pfnn38G4OzZs4SG\nhnL33XczduxYmf4UDiaTCYDs7GyGDGlYuCI7O9upnRAtTdYIPYB9jaagoIARI0awevVqNOYK2PN/\nnM54g6AzO9DUnnqv6nxQosdgib2Lq0b/Hr3Bh7179zrdr0uXLvj5+VFRUcHy5cs5deqUbIUXjZI1\nQtFSZI1QNMq+RqNR4NU/3Iom/VHY828wV9ARQFHIyK1h+XYzD76WwtCbb0ML/GX+AiZMmIDBYKCq\nqooZM2bQrVs3/vOf/zjW/2655RY3vzvR2mm1WhYuXEhSUhIJCQkkJycTExNDdnY28+bNY82aNaSm\npjqKKkieoWhpEgg9QHnuj8wbYeC+fnqu2Jpcd6FjD+h/D2kH9Ix/8SkAum3cwoDBQ8nOzub9998H\nsKep8NprrwEQGRkp63/iV0lMTCQ1NZXp06cTFxfneLz+vyXJMxTuIlOj7VX5ScheBdtT4ERdEnuN\n3h9d/4nQ/x64YhAoCllZWU4fTv+/vbsPsquu7zj+/u6V0AJJFCWwgUaiiwbNmoAdjBnBLBRTbMYs\ny1YhYweYToMM4DQT7TRhRMNYxJogTshklNry4ExQ4prKokJVtq6wUkGBrIiCJDzUJICUJDwlcPPr\nH/fu5mbJJrubvXvuzXm/Zu4k95xzb745MPvJ7+n8+kyePJnJkydz3333sWzZMk444QT/la4DMliL\nr28Me968eSxdurS/xXjVVVf1txgNQ+3PSLtGDcI6MaQuo9d3wO/vKIXfo3fCrtLmwqnhTfx4Y/Dv\nv3qFV6fM4bvrbusfp9m1axfz58/npz/9KYcddhiFQoEtW7b0f+XUqVNZvny5P4RUNY4harQ4RngQ\n22eX0dlnw9P3lcKv97vw6gu7Pzj5JJhxHjH9HLbf8TO+fUM76cEfcNppp/GFL3wBgGXLlvHzn/8c\ngJtvvpn58+c7RqMx5TpDZc0grHGVXUZr1qzp7zK6fvkV/OraBZz5yNsZv7PiSTDjJ8P7Pg4zzoVJ\nJ/Yf7hujufjii7n77rs588wz+89NmjSJ1atX97f6/GGjseQ6Q2XNIKxhxWKRxYsXM2/evFKX0c4X\n4bfrmPXILcyafi9wKOzcTDrkMOLEj5XCb+pp0LD3FlxbWxvz58/f48ky7vCgrLnOUFlzjLCGdXV1\nccbpLfR+/zpO3PFr+G0nvP5K+WzwwpEz+PQ3e/iHFf/JqX91Vqa1SiPlGKFGy0jHCBv2f4kyseVh\njn7wOp5adAQn3rcU1t9aCsG3vQvOuAIW9VK4sJObH3qNp599Yf/fJ9WovnWGnZ2dtLa20tPTw/bt\n2+np6aG1tZXOzk6WL19uCKpq7BqtJS8+Wwq8B9fA5oc4EWB8A68dMoFDZn4CZp4Hk0+G8lNgent6\nALuMVP+Gss5Qqha7RrP22qvw+x/Cg7fAo/8FqVg63nAI6YS5XLS6i+fechJrv/d9u4x00PPJMjoQ\nLp+oJynBU/eWlzx8D3Zs3X3u2L8sTXqZfg5x2JH89aEdQ3o0lXQwKBQKQ5q1bGBqNBmEY+n5DfDQ\nd0oB+H8bdh+fcBzM+AS871w46l17fMQuI2lPPopNo82u0Wp7dSv8Zl2p6/PJe3YfH3cEvGd+qfX3\n9g9Bw77nLfkvYAkfxaZ98hFrtRSExdfh8btKLb9HbofXXy2fCHjHnNJzPk+cB+MOz7BIqb64zEL7\nU3NjhBFxOfA3wExgZ0rpzXu5ZgqwGmgBXgRuBJaklF6vVl1VtXl9qeX30HfgpYpd4I+aVgq/930c\nJkzOrj6pju3rUWwpJebOncttt93GypUrueyyywxDDVk1xwjHAbcCPcDfDzwZEQXgdmAzMBtoBG4C\nXgOWVrGu0bV9S3nJwy2wZf3u44e9FZr/thSAjTP6lzxIGpnBHsU2cMxw0aJFfO1rX3PMUENWtQX1\nKaXPp5S+Cqwf5JKPAO8BPplSeiCl9EPgc8AlETGuWnWNitdegfVr4VvtcM00uPPyUggWxpXG/c67\nBRb/Ds76MkyeaQhKo6DyUWx9+sYMm5ub+cY3vgHAqlWraG5upr29nY6OjkxqVX2p+hhhRFwAXDuw\nazQirgQ+llKaWXFsKvA4cHJK6dfsRUQcChxacWg88PSBjBEOaSLKrl3w1C9K436/WQc7Krqf/+ID\npUkv7z0b/vwtI6pB0r4NHCNMKfW/7+jooK2trX+MMCIcM8yhmhsjHIJjgC0Djm2pODeYJcDnR6uI\n/U7F/tMf4KFvl7o+X3hi9wffPKW03GHGufDWd45WOZIG0fcotr51tXPnzmXjxo0sXbqUtra2N6yr\ndfsmDdWwukYj4uqISPt5TatWsWVfAiZWvI4b6RdVdqtUPt9w1sxp3HHVefzpX98PK0+G//5yKQTH\njYeT/g4u+AF8+kE4/XJDUBpDfetq169fz6WXXgrAwoUL6e3tfcPSCbdv0lANt0W4ArhhP9c8PsTv\n2gycMuDY0RXn9iqltAPY0fc+Rjj+9oYtjlIRHvsJs55cwwdOvp+Y8Wfw8mOkaCDeeXpp0su7Pwrj\nDhvRnydpdPRtJ7Zy5UoWLVrEqlWruOiii97Q/en2TRqqLMcIzwI6gcaU0jPlYwuBrwCTyoE3lO8f\n0TrCrq4uWlpaSi3Aw5+GH/0zvPRs//mXjjieKzoe4ZzPr2H2XGeeSbXGdYUaqOa2YYqIKRExE5gC\nFCJiZvl1RPmSO4GHgZsjYkZEzAW+CKwaaggeiD2mYh8xqRSChx8Fsy6Bi7rZtbCba3p28sTzVS9F\n0gi4fZNGSzUny1wJnF/xvm8WaAvQlVIqRsQ8Sgvqe4CXKC2ov6KKNfXbY1fsUz4En/wuTJ0DhdIt\ncYsjqfb5LF6Nhtw+Ys1uFeng4bN4BfW5fCJTA6diu8WRVL+Gun2TtDe5DUKwW0WSlOOu0Up2q0hS\n/bNr9ADYrSJJ+VW15ROSJNUDg1CSlGsGoSQp1wxCSVKuGYSSpFwzCCVJuWYQSpJyzSCUJOWaQShJ\nyjWDUJKUawahJCnXDEJJUq4ZhJKkXDMIJUm5ZhBKknLNIJQk5ZpBKEnKNYNQkpRrBqEkKdcMQklS\nrhmEkqRcMwglSblmEEqScs0glCTlmkEoSco1g1CSlGtvyroASTrYFItFuru72bRpE42NjZx66qkU\nCoWsy9IgbBFK0ijq6OigqamJlpYWFixYQEtLC01NTXR0dGRdmgZhEErSKOno6KC9vZ3m5mZ6enrY\nvn07PT09NDc3097ebhjWqEgpZV3DAYmICcDWrVu3MmHChKzLkZRTxWKRpqYmmpubWbduHQ0Nu9sZ\nu3btorW1ld7eXh599FG7Satk27ZtTJw4EWBiSmnbUD9ni1CSRkF3dzcbN25k6dKle4QgQENDA0uW\nLGHDhg10d3dnVKEGYxBK0ijYtGkTANOnT9/r+b7jfdepdhiEkjQKGhsbAejt7d3r+b7jfdepdjhG\nKEmjwDHC7DlGKEkZKhQKrFixgs7OTlpbW/eYNdra2kpnZyfLly83BGuQC+olaZS0tbWxdu1aFi9e\nzOzZs/uPT506lbVr19LW1pZhdRqMXaOSNMp8skw2Rto1WpUWYUQcD3wOOB04Bvgj8C3gX1JKOyuu\nmwKsBlqAF4EbgSUppderUZckjYVCocCcOXOyLkNDVK2u0WmUxh8vAh4DpgPXA4cDnwGIiAJwO7AZ\nmA00AjcBrwFLq1SXJEl7GLOu0Yj4LHBxSukd5fdnAZ3A5JTSlvKxTwFfBo6qbDnu53vtGpUk1cWs\n0YnA8xXvPwis7wvBsjuACcB7B/uSiDg0Iib0vYDxValWkpQLYxKEEdEEXAZ8veLwMcCWAZduqTg3\nmCXA1orX06NUpiQph4YVhBFxdUSk/bymDfjMscCPgFtTStePQs1fotS67HsdNwrfKUk1oVgs0tXV\nxZo1a+jq6qJYLGZd0kFvuJNlVgA37Oeax/t+ExGTgbuAe4CFA67bDJwy4NjRFef2KqW0A9hR8Wfs\npxxJqg8dHR0sXryYjRs39h87/vjjWbFihWsQq2hYLcKU0rMppUf289oJ/S3BLuB+4MKU0q4BX9cD\nNEfEpIpjZwLbgIdH/DeSpDrkXobZqcqs0YoQfAI4H+hv26eUNpevKQAPUFpj+E+UxgVvBv4tpTTk\n5RPOGpVU73xO6eiotVmjZwJNwBmUJrNsqngBkFIqAvMohWQPpQX3NwFXVKkmSapJ7mWYraosqE8p\n3cD+xxJJKT0BfLQaNUhSvXAvw2y5+4QkZcy9DLPlQ7clKWOOEY6OWhsjlCQNkXsZZsv9CCWpBriX\nYXbsGpWkGuJehiNXU/sRSpJGxr0Mx55jhJKkXDMIJUm5ZhBKknLNIJQk5ZpBKEnKNYNQkpRrBqEk\nKdcMQklSrhmEkqRcMwglSblmEEqScs0glCTlmkEoSco1g1CSlGsGoSQp1wxCSVKuGYSSpFwzCCVJ\nuWYQSpJyzSCUJOWaQShJyjWDUJKUawahJCnXDEJJUq4ZhJKkXDMIJUm5ZhBKknLNIJQk5ZpBKEnK\nNYNQkpRrBqEkKdcMQklSrhmEkqRcMwglSblmEEqScs0glCTlWtWCMCK+HxFPRsSrEbEpIm6OiMkD\nrpkSEbdHxMsR8UxEfCUi3lStmiRJo6tYLNLV1cWaNWvo6uqiWCxmXdKwVbNFeBfwceDdwDnAO4G1\nfScjogDcDowDZgPnAxcAV1axJknSKOno6KCpqYmWlhYWLFhAS0sLTU1NdHR0ZF3asFQtCFNKX00p\n/SKl9ERK6R7gamBWRBxSvuQjwHuAT6aUHkgp/RD4HHBJRIyrVl2SpAPX0dFBe3s7zc3N9PT0sH37\ndnp6emhubqa9vb2uwjBSStX/QyKOBFYDx6aUPlQ+diXwsZTSzIrrpgKPAyenlH49yHcdChxacWg8\n8PTWrVuZMGFCtf4KkqSyYrFIU1MTzc3NrFu3joaG3W2qXbt20draSm9vL48++iiFQmHM6tq2bRsT\nJ04EmJhS2jbUz1V1skxEfDkiXgL+BEwB5lecPgbYMuAjWyrODWYJsLXi9fToVCtJGoru7m42btzI\n0qVL9whBgIaGBpYsWcKGDRvo7u7OqMLhGVYQRsTVEZH285pW8ZGvACdR6gYtAjdFRBxgzV8CJla8\njjvA75MkDcOmTZsAmD59+l7P9x3vu67WDXeG5grghv1c83jfb1JKzwHPAb+PiN8CTwGzgB5gM3DK\ngM8eXf5182BfnlLaAezoe3/guSpJGo7GxkYAent7mTVr1hvO9/b27nFdrRuTMUIoLZUAngBaUkpd\nEXEW0Ak0ppSeKV+zkFIrclI58IbyvROArY4RStLYcIxwCCLiAxFxaUTMjIi3R8TpwBrgD5RagwB3\nAg8DN0fEjIiYC3wRWDXUEJQkjb1CocCKFSvo7OyktbV1j1mjra2tdHZ2snz58jENwQNRrckyLwNt\nwE+A3wHfBB4CPtwXcimlIjCP0thhD/At4CbgiirVJEkaJW1tbaxdu5b169cze/ZsJkyYwOzZs+nt\n7WXt2rW0tbVlXeKQjVnXaLXYNSpJ2SkWi3R3d7Np0yYaGxs59dRTM2sJjrRr1MeZSZJGrFAoMGfO\nnKzLOCA+dFuSlGsGoSQp1wxCSVKuGYSSpFwzCCVJuWYQSpJyzSCUJOWaQShJyjWDUJKUawfNk2W2\nbRvy03QkSQehkebAwfCs0WNxl3pJ0m7HpZT+d6gXHwxBGMBkYHvWtVTReEphfxwH999zpLw/g/Pe\nDM57s2/1en/GA39Mwwi3uu8aLf9lh5z89aiU9QBsH84T1fPC+zM4783gvDf7Vsf3Z9i1OllGkpRr\nBqEkKdcMwvqwA1hW/lVv5P0ZnPdmcN6bfcvN/an7yTKSJB0IW4SSpFwzCCVJuWYQSpJyzSCUJOWa\nQShJyjWDsI5ExPER8c2I2BARr0TEHyJiWUSMy7q2WhERl0fEPRHxckS8kHU9WYqISyJiY0S8GhH3\nRsQpWddUCyLitIi4LSL+GBEpIlqzrqlWRMSSiPhlRGyPiGciYl1EvDvruqrNIKwv0yj9N7sIeC+w\nCPgUcFWWRdWYccCtwOqsC8lSRHwCuIbSOrCTgQeBOyJiUqaF1YbDKd2PS7IupAZ9GFgFzALOBA4B\n7oyIwzOtqspcR1jnIuKzwMUppXdkXUstiYgLgGtTSm/OupYsRMS9wC9TSpeW3zcATwErU0pXZ1pc\nDYmIBJydUlqXdS21KCKOAp4BPpxS+lnW9VSLLcL6NxF4PusiVDvKXeXvB37cdyyltKv8/oNZ1aW6\nNLH860H9M8YgrGMR0QRcBnw961pUU94GFIAtA45vAY4Z+3JUj8q9CNcCd6eUerOup5oMwhoQEVeX\nB+339Zo24DPHAj8Cbk0pXZ9N5WNjJPdH0gFbBUwHzs26kGqr+/0IDxIrgBv2c83jfb+JiMnAXcA9\nwMLqlVUzhnV/xHNAETh6wPGjgc1jX47qTURcB8wDTkspPZ11PdVmENaAlNKzwLNDubbcErwLuB+4\nsDz2c1Abzv0RpJR2RsT9wBnAOujv5joDuC7L2lTborQb70rgbGBOSmlDxiWNCYOwjpRDsAt4AvgM\ncFTfLtIpJf+lD0TEFOBIYApQiIiZ5VOPpZRezK6yMXcNcGNE3Af8D/CPlJYN/EemVdWAiDgCaKo4\nNLX8/8nzKaUnMyqrVqwCFgDzge0R0TemvDWl9Ep2ZVWXyyfqSHlJwF5/kKWUYmyrqU0RcQNw/l5O\ntaSUusa2mmxFxKXAZylNkHkA+HRK6d5sq8peRMyh1Ksy0I0ppQvGtpraUl5OsjcXppRuGMtaxpJB\nKEnKNWeNSpJyzSCUJOWaQShJyjWDUJKUawahJCnXDEJJUq4ZhJKkXDMIJUm5ZhBKknLNIJQk5ZpB\nKEnKtf8HQMgfi7yTk9oAAAAASUVORK5CYII=\n", | |
"text/plain": [ | |
"<matplotlib.figure.Figure at 0x21b25086ba8>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"fig, ax = plt.subplots(figsize=(5,5), dpi=100)\n", | |
"\n", | |
"ax.plot(x, y, 'ko', ms=6, mfc='none')\n", | |
"ax.plot(xnew, y_scipy, '-', label='OLS')\n", | |
"ax.plot(xnew, y_sieg, '-', label='Siegel')\n", | |
"ax.legend()" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": { | |
"collapsed": true | |
}, | |
"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.1" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 2 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment