Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Linear regression example
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linear regression with simple basis extension"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear regression on fake BMI data"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# fake BMI data\n",
"\n",
"import numpy.random as random\n",
"\n",
"def generate_fake_bmi_data(N=100,noise=2,seed=2):\n",
" random.seed(seed)\n",
" ages=random.uniform(low=15,high=80,size=N)\n",
" heights=random.normal(1.75,0.15,size=N)\n",
" bmis=random.normal(30,5,size=N)\n",
" weights=heights*heights*bmis+noise*random.normal(0,1,N)\n",
" return (ages,heights,bmis,weights)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def compute_weights(x_train,y_train):\n",
" #w=np.dot(np.dot(np.linalg.inv(np.dot(x_train,x_train.transpose())),x_train),y_train)\n",
" w=np.dot(np.linalg.pinv(x_train.transpose()),y_train)\n",
" return w\n",
"\n",
"def compute_ridge_weights(x_train,y_train,l):\n",
" w=np.dot(np.dot(np.linalg.inv(l*np.eye(x_train.shape[0])+np.dot(x_train,x_train.transpose())),x_train),y_train)\n",
" #w=np.dot(np.linalg.pinv(x_train.transpose()),y_train)\n",
" return w\n",
"\n",
"def compute_mse(w,x_test,y_test):\n",
" y=np.dot(x_test.transpose(),w)\n",
" mse=np.dot((y-y_test).transpose(),(y-y_test))/y.size\n",
" if len(mse.shape)==2:\n",
" mse=mse[0][0]\n",
" \n",
" return mse,y"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"# generate higher order features by padding products of features up to max_order\n",
"\n",
"def gen_high_order_features(x1,x2,x3,max_order=2):\n",
"# assert(xs.shape[0] == 4) # 'Only work with 3 basics features'\n",
" \n",
" from sympy.abc import x\n",
" from sympy.abc import y\n",
" from sympy.abc import z\n",
" \n",
" mat=[np.ones(x1.size)]\n",
" p=sympy.poly(x+y+z)\n",
" \n",
" for order in range(1,max_order+1):\n",
" for mon in p.monoms():\n",
" mat.append(x1**mon[0]*x2**mon[1]*x3**mon[2])\n",
" p=p*p\n",
" return np.array(mat)\n",
"\n",
"def bmi_train(ages,heights,bmis,weights,order):\n",
" N=ages.size\n",
" x_train=gen_high_order_features(ages,heights,bmis,order)\n",
" y_train=weights.reshape(N,1)\n",
" ws=compute_weights(x_train,y_train)\n",
" mse,y=compute_mse(ws,x_train,y_train)\n",
" \n",
" return (ws,mse,y)\n",
"\n",
"def bmi_test(agest,heightst,bmist,weightst,ws):\n",
" N_test=agest.size\n",
" x_test=gen_high_order_features(agest,heightst,bmist,order)\n",
" y_test=weightst.reshape(N_test,1)\n",
" \n",
" mse_test,y=compute_mse(ws,x_test,y_test)\n",
" return (mse_test,y_test)"
]
},
{
"cell_type": "code",
"execution_count": 571,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[ 6.62858971 1.00510313 0.3199968 ]\n",
"[ 8.64919034 4.7262191 7.70046893]\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7f9debd33990>"
]
},
"execution_count": 571,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEPCAYAAABFpK+YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FOX2wPHvCR0hSO/SroI0BRSCKEZRBAsCUq0XRUBF\nkijWqz/g2q4KJkFFERHligqKNEE0F8kFS+ggJXhBVEoIWJBOIMn5/TG7MYRUks3sZs/nefZhdnZ2\n5mSynJ28Z973FVXFGGNM8AhxOwBjjDHFyxK/McYEGUv8xhgTZCzxG2NMkLHEb4wxQcYSvzHGBBmf\nJ34RiRCRjZ7HKF8fzxhjTO58mvhFpBVwD3AJcDFwo4g09eUxjTHG5M7XV/wXAitUNUVV04BlQF8f\nH9MYY0wufJ34NwFXiEhVEakIXA809PExjTHG5KK0L3euqltF5EUgDjgCrAPSfHlMY4wxuZPiHKtH\nRJ4Ddqnqm1nW24BBxhhTQKoqZ/O+4rirp6bn3/OAPsAH2W2nqvYogseYMWNcj6EkPex82vn010dh\n+LSpx2O2iFQDTgH3q+qhYjimMcaYHPg88atqV18fwxhjTP5Zz90SJjw83O0QShQ7n0XLzqd/KNbi\nbo5BiKg/xGGMMYFCRNCzLO4WRxu/MaYING7cmF9++cXtMEwxa9SoET///HOR7tOu+I0JEJ4rPLfD\nMMUsp997Ya74/aaNf+O+jW6HYIwxQcFvEn+36d14YOED/H7sd7dDMcaYEs1vEn/iA4mICBe+fiGv\nr3yd1PRUt0MyxpgSyW8Sf/WK1Xnt+tdYcucSZifOpt3kdnz101duh2WMcUF6ejqVK1dm9+7dRbqt\ncfhlcVdV+TTxU0bHjaZ93faMv3Y8Tao2cTFCY9znz8XdypUrI+LUGY8ePUq5cuUoVaoUIsLkyZMZ\nPHiwyxEGLl8Ud/0y8XsdP3WcCd9NIDohmvsvuZ/HL3+cc8qe40KExrjPnxN/Zk2bNmXq1KlcddVV\nOW6TlpZGqVKlijEq38vuZyroz5nd9iX6rp7sVChTgae6PsWGERv48cCPtHi9BR9s/CAgPvzGBKvs\nBhF7+umnGTRoELfeeitVqlRhxowZJCQk0LlzZ6pWrUr9+vWJiIggLc0ZtT0tLY2QkBB27twJwB13\n3EFERATXX389oaGhdOnSJaNPQ0G2Bfj8889p3rw5VatWZdSoUVx++eVMnz49x5/l+eef529/+xu1\natXi1ltv5eDBgwD8+OOPhISE8O6779KoUSOuu+66bNcBzJkzh9atW1OtWjWuueYa/ve//2Uco2HD\nhowfP562bdtSqVKlovgV5M3tEeY8HxDNj+W/LNf2k9trl6lddE3Smny9x5iSIr//T9zWuHFjXbJk\nyWnrnnrqKS1XrpwuXLhQVVVPnDihq1ev1pUrV2p6err+9NNP2rx5c3399ddVVTU1NVVDQkL0l19+\nUVXV22+/XWvWrKlr167V1NRUHThwoN5xxx0F3nbfvn1auXJlXbBggaampuorr7yiZcuW1ffeey/b\nn2X8+PF6+eWX6969ezUlJUXvvffejH1t375dRUTvvvtuPX78uJ44cSLbdYmJiVqpUiWNj4/X1NRU\nff7557V58+aampqqqqoNGjTQSy65RJOSkvTEiRNnxJDT792z/uxy7tm+sSgfBflAp6al6pQ1U7T2\ny7V16Lyhuu/Ivny/15hAlp//J1D4R2HllPi7deuW6/vGjx+vAwYMUFUnmYvIacn8vvvuy9h2/vz5\n2qZNmwJv+84772jXrl1PO27dunVzTPznn3++Llu2LOP5zp07tXz58qrqJP6QkBDdvXt3xuvZrRsz\nZozedtttGc/T09O1Tp06+s0336iqk/jff//9HM+LLxK/Xzf1ZKdUSCmGth/K1pFbCS0XSqtJrXjl\nu1c4mXbS7dCMcV1RpH5fadjw9FlXf/jhB2688Ubq1q1LlSpVGDNmDL/99luO769Tp07GcsWKFTly\n5EiBt01KSjojjgYNGuS4n507d3LTTTdRrVo1qlWrRtu2bQkJCWH//v0Z29SvX/+M92Vel5SURKNG\njTKeiwgNGjRgz549+YrBFwIu8XudW/5cJlw3geVDlhO3I462b7Rl8fbFbodljMmB964fr+HDh9Om\nTRt27NjBwYMHGTdunLcFwGfq1q3Lrl27TluXOQFn1bBhQ+Li4vjjjz/4448/OHDgAEePHqVWrVr5\nPma9evVOqzGoKrt37z4t2Wc9N74WsInfq0WNFiy6dRHju4/nwc8f5MYPbmTb79vcDssYk4fDhw9T\npUoVKlSoQGJiIpMnT/b5MW+88UbWrVvHwoULSUtLIyYmJte/MoYPH84TTzyR8WWxf/9+FixYkPF6\ndl9UWdcNGDCA+fPns2zZMlJTU3nppZcIDQ2lY8eORfRTFVxxTL0YJSKbROR7EZkhImV9cAxuvOBG\nNt23ia6NutJ5amcejXuUQyk22ZcxxS2/V68TJkzg3XffJTQ0lPvuu49BgwbluJ+89pnfbWvVqsXM\nmTOJioqiRo0a/PTTT7Rr145y5cplu/3DDz9Mz5496datG1WqVOHyyy9n9erVuR4r67qWLVvy3nvv\nMWLECGrVqsWXX37J/PnzM27bLO6rffDxffwiUg/4GmihqidFZCawUFWnZ9lOizKO5CPJPLHkCb7Y\n/gXPd3ueOy+6kxAJ+D9uTJALlPv4A0l6ejr16tVj9uzZdOnSxe1wshWo9/GXAs4RkdJARSDJ1wes\nU6kO026extxBc3lz9Zt0ntqZFbtX+PqwxpgA8MUXX3Dw4EFSUlL45z//SdmyZV1tdnGDTxO/qiYB\nE4CdwB7gT1X9jy+PmVnH+h359p5vGXnpSPrO6sudc+4k6bDPv3eMMX7s66+/pmnTptSuXZu4uDjm\nzp1LmTJl3A6rWPm6qedcYDbQHzgIfAJ8rKofZNlOx4wZk/E8PDy8yOfmPJxymOeXP8+UtVMYfdlo\nosKiKFc6+3Y9Y/yRNfUEJ+/vPT4+nvj4+Iz1nrug/G+sHhHpB1ynqvd6nt8BdFLVkVm2K9I2/txs\n/2M7D3/5MJv3b2ZC9wn0at7LleKKMQVliT84BdwgbSLSEZgKXAqkANOAVar6epbtii3xe33545dE\nLo6kQWgDYnrE0LJmy2I9vjEFZYk/OAVccVdVV+I076wDNgACvOXLY+ZX92bd2TBiAzecfwNXvnsl\nkYsjOXD8gNthGWOMz/n1sMzF5dejv/L00qeZu3Uu48LHMbT9UEqFlKwhY03gsyv+4BRwTT35DsLl\nxO+1Pnk9oz4fxeGTh4ntEUvXRl3dDsmYDJb4g5Ml/mKgqszaPItH4h6hc8POvHzty5xX5Ty3wzLG\nEj9w8uRJqlevzrZt204biK0kC7g2/kAkIgxsPZCtI7fSonoL2k1ux7j4cRw7dczt0IzxW5UrVyY0\nNJTQ0FBKlSpFxYoVM9Z9+OGHZ73fzp0788EHf939XbZsWQ4fPhw0ST87qempzNo8q1D7sMSfg4pl\nKjLuqnGsHbaWzb9upuXrLfl488dBf8VlTHYOHz7MoUOHOHToEI0aNWLhwoUZ60rifLvemcLyWlfQ\nfeTmzxN/8vI3L9M0timvrXytQO/NyhJ/Hhqd24hZ/Wfxbu93eXb5s1z13lVsSN7gdljG+C3vZB+Z\npaen88wzz9CsWTNq1arFHXfcwaFDziCKx44dY/DgwVSvXp2qVavSuXNnDh48yOjRo1m1ahVDhw4l\nNDSURx55hJSUFEJCQkhKcnrgDx48mKioKHr06EFoaChXXHHFacMuL1y4kAsuuIBq1aoRFRV1xl8Q\n+Y3xhx9+oEyZMrz99tucd9553HDDDdmuA5g9ezatWrWiWrVqdO/ene3bt2cco27dukyYMIHWrVtT\npUqVfJ/TkYtG0iS2CRv2bWDOwDksG7Is3+/N1tnO4FKUDwJkSrlTaaf0jVVvaK2Xa+mIBSP016O/\nuh2SCSKB8v8kuxm4/vWvf2nXrl01OTlZU1JSdMiQIXr33XerqmpsbKz2799fU1JSNC0tTVevXq3H\njh1TVdWwsDD94IMPMvZz4sQJDQkJ0T179qiq6qBBg7R27dq6fv16TU1N1X79+umQIUNUVTUpKUkr\nVaqkixYt0tTUVH3ppZe0bNmyOmPGjGzjzi5G7762bt2qIqL33ntvxpSK2a37/vvvtXLlyrps2TI9\ndeqUPvPMM9qyZUtNS0tTVdU6depox44dNTk5OdtpFrMD6JP/eVJ3H9x9xno925x7tm8sykegfKC9\nfj/2uz646EGt+VJNnZgwUU+mnnQ7JBME8vP/hLEU+lFY2SX+Jk2a6LfffpvxfMeOHVqxYkVVVZ00\naZKGh4frpk2bzthXWFjYaYn6xIkTKiKnJf4HH3ww4/VPP/1U27Vrp6qqb731ll599dUZr6Wnp2ut\nWrVyTPy5xbh161YNCQnR5OTkjNezW/ePf/xD77rrroznaWlpWrNmTV2xYoWqOon/o48+yvb4Ocnp\n916YxF+6cH8vBKdqFaoxsedEhnUYRuTiSN5c8yaxPWK5puk1bodmgpyO8c8a1K5du7j++uszhkdR\nT1PQH3/8wT333ENycjL9+vXj6NGj3HHHHTz77LP5Hkolv9Msiki20yTmJ0aAkJAQateufdp7sq7L\nOs1iSEgI9evXd3WaxexYG38htK7Vmrg74nju6ucYtmAYfWb2YceBHW6HZYzfadCgAV999dUZUxhW\nq1aNsmXLMm7cOBITE1m2bBkff/wxH330EVC4SUqyTrOoqrlOs5hbjDnFknVd1mkW09PT2bNnj6vT\nLGbHEn8hiQi9W/RmywNbuLTepVw65VKeXPIkR07mPBG0McFm+PDhPPbYY+zevRtwpjD87LPPAFiy\nZAmJiYmoKpUqVaJ06dIZs1PVrl2bHTvO7mKqV69erFy5ksWLF5OWlsaECRP4888/zypGyN80iwMH\nDmTOnDl8/fXXpKam8sILL1CjRg06dOhwVj+Dr1jiLyLlS5fnySue5PsR37Pr0C5avNaC979/327/\nNEEnuyvaxx57jGuvvZarr746YwrDdevWAc5k5zfffDOhoaG0bduWG2+8kQEDBgAQFRXFe++9R/Xq\n1Xn88cfP2H9uV8916tThww8/5MEHH6RmzZokJSXRpk2bHKdZzC3GnI6VdV2bNm2YOnUqw4YNo1at\nWixdupR58+YREhKSZ7zFyXru+si3u74lYnEEZULKMLHnRC6pd4nbIZkAZz13CyctLY06derw2Wef\n0alTJ7fDyTfruRtALmt4GSuGrmBo+6Hc9OFN3D3vbpKPJLsdljFBZfHixRw6dIgTJ04wduxYzjnn\nHL9rdnGDJX4fCpEQ7m53N1sf2Eq1CtVoPak1478dz8m0k26HZkxQWLZsGU2aNKFOnTosXbqUOXPm\nULq03cxoTT3F6IfffuChLx9i2+/biL4umhsuuMHtkEwAsaae4GSjc5YQi7YtIuqLKJpVbUb0ddE0\nr9Hc7ZBMALDEH5wCro1fRC4QkXUistbz70ERGeXLYwaC68+/no33baRbk250eacLo78czcETB90O\nyxgTJIrtil9EQoDdOJOt78ryWlBd8We278g+nlzyJIu2L+K5q5/j7xf/nRCx0os5k13xB6eAbuoR\nke7A06p6RTavBW3i91qdtJpRn4/iZNpJJvacyGUNL3M7JONnGjdufFqvUBMcGjVqxM8//3zG+kBJ\n/FOBNao6KZvXgj7xg9ML8IONH/DYfx4jvHE4L17zIvVDcx5bxBjjO1t/20pMQgwzN8+kd4veRHaK\n5KI6F7kdVobCJP5iua9JRMoAvYDHc9pm7NixGcvh4eGEh4f7PC5/IyLc1vY2bm5xMy8sf4G2b7bl\nobCHePiyhylfurzb4RlT4qkqcTviiEmIYc3eNYzoMILEBxKpU8n9Gb/i4+OJj48vkn0VyxW/iPQC\n7lfVHjm8blf82dhxYAcPf/kwG5I3MKH7BHq36O03Xb6NKUmOnzrOjI0ziEmIQUSICovi1ja3+vUF\nl9839YjIh8BiVX0vh9ct8efiPzv+Q+TiSOpUqkNsj1ha1WrldkjGlAjJR5KZtGoSk9dM5pJ6lxAV\nFkW3Jt0C4gLLrxO/iFQEfgGaqurhHLaxxJ+H1PRU3lj1Bs8se4aBrQYy7qpxVKtQze2wjAlI65PX\nE50Qzfwf5jO49WAiOkUEXH8av078+QrCEn++/XbsN/5v6f8xO3E2Y68cy7AOwygVUsrtsIzxe2np\naSzctpDohGi2/b6NkR1HMqzDsIC9gLLEH4Q2JG8gYnEEB04cILZHLOGNw90OyRi/dOTkEaatm0bs\niliqVahGVFgU/Vr2o0ypMm6HViiW+IOUqvLJlk94JO4RLq1/KeOvHU+jcxvl/UZjgsDOgzt5dcWr\nTFs/jfDG4USFRXFZw8sCov0+P/x2yAbjWyJC/1b9SXwgkTa12tD+rfaMWTqGY6eOuR2aMa75btd3\nDPh4AO0mtyNd01l17yo+GfAJXc7rUmKSfmHZFX8JsvPgTh6Ne5Rvd33LS9e+xMBWA+2DboJCanoq\ns7fMJjohmv1H9xPRKYIh7YYQWi7U7dB8pkQ09aSnK5ajisbyX5YzavEoKpetTGyPWNrVbed2SMb4\nxIHjB5iydgqvrXyNJlWbENkpkl7NewXFDQ8loqnnzjvhxAm3oygZrmh0BavvXc3tbW+n54yeDF8w\nnF+P/up2WMYUmW2/b2PkopE0m9iMjfs3MmfgHP779//S58I+QZH0C8tvEv/Jk3DVVZBssxMWiVIh\npRjWYRiJDyRSoUwFWk5qSUxCDKfSTrkdmjFnRVVZ+tNSen3Yiy7vdKFKuSpsun8T/+7zbzrUs+kU\nC8Kvmnr++U945x2YNw8uvtjtqEqWLb9uIXJxJLsP7SamRwzdm3V3OyRj8iUlNYUPN31ITEIMJ9NO\nEhkWye1tb6dimYpuh+aqEtHG741j1ix44AGYMgV693Y5sBJGVVnwvwU89MVDtKrVigndJ/C3an9z\nOyxjsrX/6H7eXP0mb6x+g7a12xIVFkX3Zt1tvgqPEpX4AVavhj594P774fHHsaJvEUtJTSE6IZrx\n345naPuh/OOKf1C5XGW3wzIGgE37NxGTEMPsxNn0b9mfiE4RNj5VNkpc4gfYs8e54m/Rwrn6L++/\ng+QFrKTDSTyx5An+s+M/vNDtBW5ve7tdTRlXpGs6i7cvJjohms37N3P/pfczvMNwap5T0+3Q/FaJ\nTPwAx47BkCGwcyfMmQN13B8Su0RK2J3AqM9HESIhTOw5kY71O7odkgkSR08e5d/f/5uYhBgqlqlI\nVFgUA1oNoFzpcm6H5vdKbOIHUMWKvsUgXdOZvmE6Ty55ku7NuvNCtxeoW7mu22GZEmrPoT28tvI1\n3l73Nl0adiEqLIqujbpah8MCKBH38edEBMaMgfHj4dprYe5ctyMqmUIkhL9f/He2jtxK7XNq0+aN\nNrz0zUukpKa4HZopQVYnrea2T2+jzRttOHbqGN/d8x1zB83lysZXWtIvRn5/xZ+ZFX2Lz7bft/HQ\nlw+x9betvNL9FW684Eb7j2nOSlp6GnO3ziU6IZpdh3YxquMo7ml/D+eWP9ft0AJaiW7qycpb9G3e\nHN5+24q+vrZ4+2IiF0fS+NzGxPSIoUWNFm6HZALEoZRDTF07lYkrJ1Kvcj2iwqLo3aI3pUOKZarv\nEs+vm3pEpIqIfCwiiSKyWUQ6FWZ/9evDf/8Lp05ZT9/i0ONvPdh430aua3YdV0y7goe+eIg/T/zp\ndljGj+04sMO5WIhpzMqklczsN5Nv7v6Gfi37WdL3E8XRxh8LLFLVC4GLgMTC7rBiRfjoI+jRAzp1\ngvXrCx2jyUWZUmWI6hzF5vs3czjlMC1ea8GUNVNIS09zOzTjJ1SV5b8sp+/MvnSc0pHypcuzYcQG\nPrzlQ7tLzA/5tKlHREKBdaraLI/tznpY5o8/dtr833rLaf83vrd271pGfT6KY6eOMbHnRC4/73K3\nQzIuOZl2klmbZxGTEMOhlENEhkVy10V3cU7Zc9wOrcTz2zZ+EbkIeAvYgnO1vxqIUNXjWbYr1Hj8\n3qLvfffBE09Y0bc4qCofbfqIR//zKFecdwUvXvMiDas0dDssU0x+P/Y7k9dM5vVVr3NhjQuJDIvk\n+vOvtw6AxcifE38HIAHorKqrRSQGOKiqY7Jsp2PG/LUqPDyc8PDwAh3Lir7uOHryKC9+8yKvr3qd\nyE6RjL5sNBXKVHA7LOMjib8mErsilpmbZ9KnRR8iwyJpW7ut22EFhfj4eOLj4zOejxs3zm8Tf23g\nO1Vt6nl+OfCYqt6UZbsimYHLevq65+c/f2b0l6NZs3cN468dT98L+9rtnyWEqhK3I47ohGjW7V3H\niEtGcN8l91G7Um23QwtqfnvFDyAi/wXuVdX/icgYoKKqPpZlmyKbetF6+rrrq5++ImJxBDUr1iS2\nRyxtardxOyRzlo6fOs6MjTOISYghREKICoticJvBlC9tf077A39P/BcBbwNlgB3AEFU9mGWbIp9z\n14q+7klNT+WtNW8xNn4s/Vv2559X/ZPqFau7HZbJp72H9zJp1STeWvsWl9a7lKiwKK5ucrX9Bedn\n/Drx5ysIH022vmaN0+5vRV93/H7sd8bEj2HW5lmMuXIMwy8Zbvdx+7H1yeuJTohmwQ8LGNx6MBFh\nEVxQ/QK3wzI5sMSfi6QkuPlmK/q6aeO+jUQsjuDXY78S2yOWq5tc7XZIxiMtPY3P/vcZ0QnR/Hjg\nR0ZeOpJ7O9xLtQrV3A7N5MESfx68Rd9ffnEGebOib/FTVT5N/JTRcaNpX7c9468dT5OqTdwOK2gd\nOXmEaeumEbsiluoVqxMVFsUtF95CmVJl3A7N5JNfD9ngD7w9fa+/3nr6ukVEuKXlLWy5fwvt6rTj\nkimX8PRXT3P05FG3QwsqOw/u5JEvH6FxTGOW7VzG9D7TSbgngUGtB1nSDyJBccWfmbfoO3ky9O1b\nLIc02dh9aDePxj3K8p3LefGaFxncerAVD33ou13fEZ0QzZKfljDk4iGM7DiSxuc2djssUwjW1FNA\n3qLviBHw5JNW9HXT1zu/JmJxBBVKVyC2Rywd6nVwO6QSIzU9ldlbZhOdEM2vx34lolMEQy4eYvMr\nlxCW+M9CUpKT/C+4wIq+bktLT2Pa+mk89dVT3HjBjTzf7XlqnVPL7bAC1oHjB5iydgqvrXyNplWb\nEhkWyU0X3ESpkFJuh2aKkLXxn4V69ZzhnVNTITzchnd2U6mQUgxtP5StI7dSpVwVWk1qxSvfvcLJ\ntJNuhxZQtv2+jZGLRtJsYjM27d/E3EFzif97PL1b9Lakb04TtFf8XqrwzDPOVf+8edCunSthmEy2\n/raVqC+i+OnAT0RfF03P83u6HZLfUlWW/ryUmIQYEnYnMKzDMO6/9H7qVa7ndmjGx6yppwhY0de/\nqCoLty0k6osomldvTvR10Zxf/Xy3w/IbKakpfLjpQ2ISYjiVforITpHc3vZ2GyAviFjiLyJW9PU/\nKakpxK6I5aVvXuLudnfzVNenCC0X6nZYrtl/dD9vrn6TN1a/wUW1LyIqLIruzbrbHVFByBJ/EfIW\nfc8/32n+qWAXUH4h+UgyTyx5gi+2f8FzVz/HXRffFVRjv2/ct5GYhBg+3fop/Vv2JzIskpY1W7od\nlnGRJf4idvy409P355+tp6+/WblnJaM+H0W6pjOx50TCGoS5HZLPpGs6i7cvJjohms37N/PApQ8w\n/JLh1KhYw+3QjB+wxO8DqvDsszBlihV9/U26pjPj+xk8vuRxujXpxr+u+VeJKmYePXmU6RumE7si\nloplKhIVFsXA1gMpW6qs26EZP2KJ34c++cQZ3dOKvv7ncMphnl/+PG+tfYvRnUcT1TkqoMeK33No\nD6+tfI23173N5eddTlRYFFecd4W135tsWeL3MSv6+rftf2xn9Jej2bh/I690f4VezXsFVLJctWcV\n0QnRLN6+mDva3sGoTqNoVq2Z22EZP+fXiV9EfgYOAunAKVXtmM02fp34wYq+gSDuxzgiFkfQILQB\nMT1i/Lr4mZaextytc4lOiGb3od2M6jSKe9rdQ5XyVdwOzQQIf0/8O4AOqnogl238PvGDU/S9+274\n6Scr+vqrU2mnmLRqEs8uf5ZbW9/K2PCxVK1Q1e2wMhw8cZCp66by6spXqVe5HlFhUfRu0dsmqDEF\n5u9DNkgxHcfnKlSADz6AG26Ajh1h3Tq3IzJZlSlVhoiwCLbcv4WUtBQufP1CJq+eTFp6mqtx7Tiw\ng8jFkTSJbcLqpNXM7DeTb+7+hn4t+1nSN8WuuK74/wTSgLdUdUo22wTEFX9mVvQNDOuT1zPq81Ec\nPnmY2B6xdG3UtdiOraos37mc6IRolv+ynKHthzKy40gahDYothhMyeXvTT11VXWviNQE4oCRqvp1\nlm0CLvGDFX0Dhaoya/MsHol7hM4NO/PytS9zXpXzfHa8k2knmbV5FtEJ0Rw5eYSIThHcddFdnFP2\nHJ8d0wQfv078px1MZAxwWFVfybJex4wZk/E8PDyc8PDwYourMKzoGziOnTrGi1+/yGurXmNUx1E8\n0uURKpapWGT7/+3Yb0xePZlJqydxYY0LiQqLouf5PYOqh7Hxnfj4eOLj4zOejxs3zj8Tv4hUBEJU\n9YiInAN8CYxT1S+zbBeQV/xe3qLvjh1O0bduXbcjMrn55c9feCTuEVbsWcHL175M/5b9C3X7Z+Kv\nicQkxDBryyz6tOhDZFgkbWu3LcKIjTmT317xi0gTYA6gQGlghqr+K5vtAjrxg/X0DUTxP8cTsTiC\nquWrEtsjlovqXJTv96oqcTviiE6IZt3eddx3yX2MuGQEtSvV9mHExvzFbxN/voMoAYnfy4q+gSU1\nPZW3177NmPgx9G3Rl2eufibXsXCOnzrOjI0ziEmIIURCiAqLYnCbwQHdY9gEJp/dzikit2da7pLl\ntZFnc8CSrl8/+OILiIx0/gIoId9nJVbpkNKMuGQEiQ8kUqZUGS58/UImrpjIqbRTp2239/Benv7q\naRrHNmYBy1N2AAAUBklEQVTeD/OY2HMiG0ZsYEi7IZb0TcDJ9YpfRNaqavusy9k9L1QQJeiK32vv\nXrj5Zvjb32DqVCv6BopN+zcRuTiSvUf2EtsjlhoVaxCdEM2CHxYwuPVgIsIiuKD6BW6HaYzvmnpE\nZJ2qtsu6nN3zwiiJiR+s6BuoVJV5P8zjoS8e4lT6KR7s+CBD2w+lWoVqbodmTAZfJn674i8kK/oG\nLm9vX5uo3PgjXyb+Y8B2nGEXmnmW8TxvqqpF0iOlJCd+r9mznY5eb74Jt9zidjTGmEBXmMSf1yAh\nF57NTs2ZbrkFmjRxOnslJsI//mE9fY0x7ijQ7ZwiUh3oCuxU1TVFFkQQXPF7WdHXGFMUfHk752ci\n0tqzXBfYBNwN/FtEIs/mgMGubl3473+dtv/wcOeLwBhjilNeg4g0UdVNnuUhQJyq3gR0wvkCMGfB\nO7zzTTdBp06wdq3bERljgkleiT9zL5ZuwCIAVT2MM6OWOUsi8NRTEB0N113nFH+NMaY45FXc3SUi\nDwK7gfbAYgARqQCU8XFsQcGKvsaY4pbXFf89QCvg78BAVf3Tsz4MmObDuIJK+/awYgUsWAC33eZ0\n/DLGGF+xQdr8iLen748/Op29rKevMSYnvuzANT+3N6tqr7M5aDbHscTvoQrPPQdvveUM89C+SPpG\nG2NKGl8m/l+BXcCHwAqcHrsZVPW/Z3PQbI5jiT8L6+lrjMmNLxN/KeBaYDDQFlgIfKiqm8/mYLkc\nxxJ/NtaudYq+997r3AFkRV9jjFexTMQiIuVwvgBexpk+8bUCBBgCrAZ2Z9c8ZIk/Z3v3Osm/aVN4\n5x3r6WuMcfis565n5+VEpC/wPvAAMBFnOsWCiAC2FDw8U7cuxMc7V/tXXmk9fY0xhZfXkA3Tge9w\n7uEfp6qXquozqronvwcQkQbA9cDbhYo0iFWoADNmOGP8WE9fY0xh5dXGnw4c9TzNvKEAqqqheR5A\n5GPgOaAK8LA19RSOt+j7xhvONI/GmODks2GZVTXPpqDciMgNwD5VXS8i4WS5K8gU3C23OO39N9/s\n9PS1oq8xpqDyGrKhsLoAvUTkeqACUFlEpqvqnVk3HDt2bMZyeHg44eHhPg4tcLVr5/T07dMHtmyx\noq8xwSA+Pp74+Pgi2Vex9dwVkSuxpp4idfw43HMPbN/udPaqV8/tiIwxxcWnd/UY/2VFX2PM2bCx\nekoIK/oaE1x8OeeuCRBW9DXG5Jdd8Zcwe/c6Rd8mTazoa0xJZm38JkPdurB06V89fZOS3I7IGONv\nLPGXQJmLvmFhVvQ1xpzOmnpKuE8/heHDrehrTEljxV2To759nfb+m292Ons9/bQVfY0JdnbFHySs\n6GtMyWLFXZMn7/DOISFW9DUm2FniDyLly8P77zsTu3TqBGvWuB2RMcYN1tQTpLxF30mToH9/t6Mx\nxhSUFXdNgXmLvr17Oz19rehrTPCwK/4gl5zsJP/GjWHaNCv6GhMorLhrzlqdOk7Rt1Qp6NrVir7G\nBANL/Caj6NunjxV9jQkG1tRjTjNnDgwbZkVfY/ydFXdNkenTx2nvt6KvMSWXT6/4RaQcsAwoi/Ml\n84mqjstmO7vi9zNW9DXGv/ltcVdVU4CrVLUdcDHQU0Q6+vKYpmh4i76lS1vR15iSxufFXVU95lks\nh3PVb5f2AaJ8efj3v/8q+q5e7XZExpii4PPELyIhIrIOSAbiVHWVr49pio4IPPkkTJwIPXvCrFlu\nR2SMKSyfF3dVNR1oJyKhwFwRaamqW7JuN3bs2Izl8PBwwsPDfR2aKQDvyJ7eOX3/7/+s6GtMcYqP\njyc+Pr5I9lWst3OKyNPAUVV9Jct6K+4GCCv6GuMf/La4KyI1RKSKZ7kCcC2w1ZfHNL5lRV9jAp+v\n2/jrAktFZD2wAvhCVRf5+JjGx6zoa0xgs567plC8PX1ffx0GDHA7GmOCR2Gaeizxm0Jbv94p+t59\ntxV9jSkulviN67xF30aNnKJvxYpuR2RMyea3xV0TPLxF3zJlnDl99+xxOyJjTE4s8Zsi4y369u0L\nYWFW9DXGX1lTj/EJK/oa41vWxm/8krfoO2QIjBljRV9jipIlfuO3kpOd+/3PO8+KvsYUJSvuGr9V\npw4sXQply1rR1xh/YYnf+Fz58jB9uhV9jfEX1tRjipUVfY0pGtbGbwKKFX2NKTxL/CbgWNHXmMKx\n4q4JOJmLvl27WtHXmOJkid+4xlv07dfPir7GFCdr6jF+Ye5cuPdeK/oak1/Wxm9KhA0boFcvp+j7\nf/8HIfb3qDE58ts2fhFpICJfichmEdkoIqN8eTwT2C66CFasgC++gEGD4NgxtyMypmTy9TVVKvCQ\nqrYCOgMPiEgLHx/TBDBv0bdcOSv6GuMrPk38qpqsqus9y0eARKC+L49pAp8VfY3xrWJrRRWRxsDF\nOJOuG5MrEXj8cXj1VejZE2bOdDsiY0qO0sVxEBGpBHwCRHiu/M8wduzYjOXw8HDCw8OLIzTj53r3\nhiZNnJ6+iYlW9DXBKz4+nvj4+CLZl8/v6hGR0sBnwOeqGpvDNnZXj8nVvn1OT98GDeDdd62nrzF+\ne1ePxzvAlpySvjH5Ubs2fPWV0/5vRV9jCsfXt3N2AW4DrhaRdSKyVkR6+PKYpuQqXx7ee88p+nbq\nBKtWuR2RMYHJOnCZgDRvntPT99VXYeBAt6MxpvhZz10TlDZscIq+f/+7FX1N8LHEb4KWFX1NsPL3\n4q4xPmNFX2MKzhK/CXjeom///lb0NSY/rKnHlChW9DXBwtr4jcnEW/S96y5nTl8r+pqSyBK/MVl4\ni7716zvNQFb0NSWNFXeNycJb9K1QwYq+xmRlid+UWFb0NSZ71tRjgsK8eTB0qFP0HTTI7WiMKTxr\n4zcmH6zoa0oSS/zG5NO+fdC3L9SrZ0VfE9isuGtMPtWuDUuWOEXfK66A3bvdjsiY4meJ3wQdb9F3\nwABnTt+VK92OyJjiZU09JqjNn+8UfSdOtKKvCSzWxm9MIXz/PfTqBXfeCWPHWtHXBAa/beMXkaki\nsk9EvvflcYwpjLZtYcUKp+1/wAA4etTtiIzxLV9f20wDrvPxMYwpNG9P33POcXr6WtHXlGQ+Tfyq\n+jVwwJfHMKaolCvnTOYycKAVfU3JVtrtAIzxJyLw6KPQogXccIMzsXvDhs59/5kfVas62xoTiPwm\n8Y8dOzZjOTw8nPDwcNdiMaZXL1i+HOLiICkJli51/vU+jh8/88sgu0flyvYFYYpGfHw88fHxRbIv\nn9/VIyKNgAWq2jaXbeyuHhNQjh2DvXtP/zLI+vCOCJrXl0Pduk5twZiC8OvbOUWkMU7ib5PLNpb4\nTYl0+HDuXw7eR7ly+fuCKFfO7Z/I+Au/Tfwi8gEQDlQH9gFjVHVaNttZ4jdBSxX+/DPvL4fkZKhU\nKe8viDp1oEwZt38q42t+m/jzHYQlfmPylJ4Ov/+e9xfE/v1QrVreXxC1akGpUm7/VOZsWeI3xmRI\nS3OSf15fEH/84ST/vL4gqle33sz+yBK/MabATp1ymo+SknIvVB8+7DQf5fUFce65dgdTcbLEb4zx\nmRMn/vqCyO2RkpL/W1xN4VniN8a47ujR/N3iGhKSvzuYbJKc3FniN8YEBNX83+JaocKZXwZ2i+tf\nLPEbY0oUVThwIH+3uIaG5v0XRO3aJe8WV0v8xpiglJ4Ov/2W9xfEr786dyfl9QVRs2bg3OJqid8Y\nY3KR31tcDxzI/y2ubt/BZInfGGOKwMmTsG9f3l8QR45kX3PI+qhSxXdfEJb4jTGmGJ04kfcdTElJ\nTl+J/NziWqlSwWOwxG+MMX7oyJH83eJaunT+bnGtUOGvfVviN8aYAKUKhw7l/dfD3r1O3wbvF0Fc\nnCV+Y4wp0VSd8ZW8XwQ9eljiN8aYoFKYph4bc88YY4KMzxO/iPQQka0i8j8ReczXxzPGGJM7nyZ+\nEQkBXgOuA1oBg0WkhS+PGeyKajJm47DzWbTsfPoHX1/xdwS2qeovqnoK+Ai42cfHDGr2H6to2fks\nWnY+/YOvE399YFem57s964wxxrjEirvGGBNkfHo7p4iEAWNVtYfn+eOAquqLWbazezmNMaaA/PI+\nfhEpBfwAdAP2AiuBwaqa6LODGmOMyVVpX+5cVdNEZCTwJU6z0lRL+sYY4y6/6LlrjDGm+BRbcVdE\nporIPhH5PpdtJorINhFZLyIXF1dsgSavcykiV4rInyKy1vN4qrhjDCQi0kBEvhKRzSKyUURG5bCd\nfT7zkJ9zaZ/P/BORciKyQkTWec7nmBy2K9hnU1WL5QFcDlwMfJ/D6z2BhZ7lTkBCccUWaI98nMsr\ngfluxxkoD6AOcLFnuRJOXapFlm3s81l059I+nwU7pxU9/5YCEoCOWV4v8Gez2K74VfVr4EAum9wM\nTPdsuwKoIiK1iyO2QJOPcwng8sRwgUNVk1V1vWf5CJDImf1N7POZD/k8l2Cfz3xT1WOexXI4ddms\n7fMF/mz60338WTt77cE6exVGZ8+ffQtFpKXbwQQKEWmM89fUiiwv2eezgHI5l2Cfz3wTkRARWQck\nA3GquirLJgX+bPr0rh7jmjXAeap6TER6AnOBC1yOye+JSCXgEyDCc7VqzlIe59I+nwWgqulAOxEJ\nBeaKSEtV3VKYffrTFf8eoGGm5w0860wBqeoR75+Hqvo5UEZEqrkcll8TkdI4ierfqjovm03s85lP\neZ1L+3yeHVU9BCwFemR5qcCfzeJO/ELObXvzgTsho8fvn6q6r7gCC0A5nsvM7Xsi0hHntt0/iiuw\nAPUOsEVVY3N43T6f+ZfrubTPZ/6JSA0RqeJZrgBcC2zNslmBP5vF1tQjIh8A4UB1EdkJjAHK4gzh\n8JaqLhKR60VkO3AUGFJcsQWavM4l0E9E7gNOAceBgW7FGghEpAtwG7DR05aqwJNAI+zzWSD5OZfY\n57Mg6gLveYa4DwFmej6LwynEZ9M6cBljTJDxpzZ+Y4wxxcASvzHGBBlL/MYYE2Qs8RtjTJCxxG+M\nMUHGEr8xxgQZS/xBTkTSRWR6puelRORXEZl/lvu7SUQeLboIz55n+N8FbseRGxEpKyJxnuGJ+2d5\nrblnON41ItLkLPYdISLliy5aU1JY4jdHgdYiUs7z/FpOH/CpQFR1gaq+VCSRFY0i6aji6UDjC+1x\nOuK0V9WPs7zWG/hYVTuo6k9nse9IoGJB3uCZLtWUcJb4DcAi4AbP8mDgQ+8LInKpiHzruer8WkTO\n96yPFJGpnuU2IvK9iJQXkbtE5FXP+mkiMklEvhOR7Z4r8KkiskVE3sl0jMOZlm8RkWkFeX9mItJD\nRBJFZDXQN9P6ip73Jnh+ll6e9RVEZKaIbBKRTz2vt/fGJSLjPT1Qw0SkvYjEi8gqEfncO/SAiDT1\nPF8lIv8VkTMGHBORqiIyR0Q2eM5naxGpCfwbuNRzxd8k0/Y9cRL3fSKyxLPuNnEm5VgrIm+IiHjW\nTxKRlZJpog4ReRCoByzN9P7czvMbIpIAvJjNubrJs13LTMdfLyLNsvsdmADg9iQD9nD3ARwCWgMf\n44z3vQ7oimeiDJzJNEI8y92ATzzLAsTjXJWuAsI86+8CJnqWpwEfeJZ7AQeBlp7nq4G23hgyxXML\n8E5B3p/pveWAnUBTz/OZmX6O54BbPctVcCYIqQA8DLzhWd8KOAm09zxPB27xLJcGvgGqe54PwJlD\nGuA/QDPPckdgSTbneSLwtGf5KmCdZznHSUlwhuJ4yLPcAmdMllKe568Dt3uWz/X8G4IziFdrz/Md\nQNXMv+tczvP8TK/ldK4mAoMznY9ybn9+7XF2DxuW2aCqm8QZO30wsJDTB387F5juudJXPOM7qaqK\nyBDge+BNVU3IYffeNvaNQLL+NZzsZqCx5/25TcqRn/d7tQB2qOoOz/P3gXs9y92Bm0TkEc/zssB5\nOLOZxXh+ps0isjHT/lKBTz3LzXG+IOM8V9ohQJKInANcBnzsvQIHymTzc1yO5y8QVV0qItXEGbo4\nv7rhNAut8hynPOAdiGuQiNyL87upA7QENpH7oIhZZW5myulcfQf8Q0QaAHNUdXsB4jd+xBK/8ZoP\nvIwz+FuNTOufAb5S1b4i0gjnitLrAuAwTpNCTlI8/6ZnWvY+937+MrfDZy1G5uf9meWU6ATn6n3b\naSsl17x4Qj2Xt573b1LVLlneXxk4oKrtc9sRZ9YaCjoDlQDvqeo/shy/Mc5fLR1U9ZCn+SY/Bd2s\n2xzN8vyMcwX84GkOuhFYJCLDVDU+n/EbP2Jt/MabgN4Bxqnq5iyvV+Gvsb0zRv0TZ6jYWJxmoeoi\ncksBjpVVsjh3sIQAfc7i/V5bgUaZ2soHZ3rtCyBj4m/5a0Lqb/CMDinOTFBtcjjeD0BNcYa9RURK\nizMhxmHgJxHpl2nfbbOJbTlwu+f1cOBXLdhkL0twRrWs6dlHVRE5DwgFjgCHPTWHnpnec8jzuld+\nz3O250pEmqjqT6r6KjAPyO7nNAHAEr9RAFXdo6qvZfP6S8C/RGQNp39eXgFe9fy5PxR4QURqZHlv\n1qtczWH5CZwmpq+BpLN4P56fIQUYhnM1upq/mkLA+culjDhF6I3APz3rJwE1RGSTZ90mnFrCacdQ\n1VNAP5zi53qcWkhnz8u3A/d4Cp6bcOoRWY0DOojIBuB5nFpIvqlqIvAU8KVnH18CdVT1e2A9zty2\n7+OcQ68pwGJvcZf8n+dn+etcec8LwABPEXwdTj1kOiYg2bDMJqh5rn7LqGqKiDQF4oDmqprqcmjG\n+Iy18ZtgVxHnlkdvQfY+S/qmpLMrfmOMCTLWxm+MMUHGEr8xxgQZS/zGGBNkLPEbY0yQscRvjDFB\nxhK/McYEmf8H+kz9Eqe5awYAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f9deef82710>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"\n",
"N=30\n",
"N_test=20\n",
"max_order=3\n",
"\n",
"ages,heights,bmis,weights=generate_fake_bmi_data(N,2,1)\n",
"agest,heightst,bmist,weightst=generate_fake_bmi_data(N_test,2,5) # no noise when testing\n",
"#x_train=np.concatenate((np.ones(N),heights,bmis,ages),axis=0).reshape(4,N)\n",
"\n",
"mse_train_list=[]\n",
"mse_test_list=[]\n",
"for order in range(1,max_order+1):\n",
" (ws,mse,y_train)=bmi_train(ages,heights,bmis,weights,order)\n",
" (mse_test,y_test)=bmi_test(agest,heightst,bmist,weightst,ws)\n",
" mse_train_list.append(mse)\n",
" mse_test_list.append(mse_test)\n",
"\n",
"print(np.array(mse_train_list))\n",
"print(np.array(mse_test_list))\n",
"\n",
"line1,=plt.plot(range(1,max_order+1),mse_train_list,label='Training error')\n",
"line2,=plt.plot(range(1,max_order+1),mse_test_list,label='Testing error')\n",
"plt.xlabel('Maximum degree of features')\n",
"plt.ylabel('MSE')\n",
"plt.legend(handles=[line1,line2])\n",
"#plt.savefig(\"bmi0.pdf\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear regression for quadratic fitting"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# general polynomial data\n",
"import numpy as np\n",
"\n",
"def gen_poly_data(N,isRandom=True,seed=1):\n",
" if isRandom:\n",
" random.seed(seed)\n",
" x=random.uniform(0,10,N)\n",
" y=(x-3)**2+random.normal(0,2,N)\n",
" else:\n",
" x=np.linspace(0,10,N)\n",
" y=(x-3)**2\n",
" return (x,y)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def gen_poly_features(x,max_order=3):\n",
" mat=[np.ones(x.size)]\n",
" p=x\n",
" for order in range(1,max_order+1):\n",
" mat.append(p)\n",
" p=p*x\n",
" return np.array(mat)"
]
},
{
"cell_type": "code",
"execution_count": 580,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"64.0703368074\n",
"58.1449568138\n",
"1.76608000787\n",
"2.24223953871\n",
"1.76602144107\n",
"2.24943481913\n",
"1.74206072203\n",
"2.80766073752\n",
"0.734274250627\n",
"44.2056154586\n",
"0.575730627164\n",
"228.218092695\n",
"0.575420075101\n",
"194.306626484\n",
"0.55113926434\n",
"624.354036092\n",
"1.96999204674e-10\n",
"41349.6855465\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEPCAYAAABLIROyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt81NWd//HXZ3JPZr4h4RYFxctau8jNskopWKMUL2jF\n7Vq8t3Wtpe3WC9X+tH2sD1C3uq3arm59bN2VVt2HxWvxRr1WI4u23sEbYi2Vi8gtCUlIQkIm5/fH\ndwghJJmZ3L7fZN7PxyMPJt98Z+bDEOY955zvOcecc4iISOaJBF2AiIgEQwEgIpKhFAAiIhlKASAi\nkqEUACIiGUoBICKSofo1AMxssZltMbN32h0rMbNnzWyNmT1jZsX9WYOIiHSuv1sAvwVO7nDsGuB5\n59yRwAvAj/u5BhER6YT190QwMxsHPOGcm5T4/kPgeOfcFjMrAyqcc5/v1yJERGQ/QYwBjHLObQFw\nzm0GRgVQg4hIxgvDILDWohARCUB2AM+5xcxGt+sC2trViWamcBAR6QHnnCU7ZyBaAJb42uNx4FuJ\n298EHuvuzs65If11+E3H8pVv/SnpeQsXLgy81rB86bVI/7VYtXkVR91xVOD1huG1yISvVPX3ZaC/\nA14BPmdm683sIuDfgdlmtgaYlfg+Y0VzPeqaa4MuQ4a4qsYqhhcOD7oMCZl+7QJyzp3XxY++0p/P\nO5jEcmNUKQCkn1U2VFJaUBp0GRIyYRgEzmjF+R71LXVJzysvL+//YgYJvRZ7pfpaVDVWMbxgaLcA\n9HuRPgVAwIYVeDTGk7cA9Mu9l16LvVJ9LSobh34LQL8X6QviKiBpp7TIo7FVXUDSO4cccgjr1q1L\net7N3DwA1chAGTduHJ988kmP768ACFhpNEaTbQ+6DBnk1q1bl9bVHzI0mCW90rNb6gIK2IioRzN1\n6P+uiAw0BUDASgo9IgW1NDYGXYmIZBoFQMC8PI+solp27gy6EhHJNAqAgMXyYkTyFQAiqWhtbSUW\ni7Fx48Y+PTdTKQAC5uV5WF4ddcmnAogMOrFYDM/z8DyPrKwsCgsL244tWbIk7ceLRCLU1dUxduzY\nPj03U+kqoIB5eR4uTy0AGZrq2n2yOeyww1i8eDEnnHBCl+fH43GysrIGorQB09nfKd2/Z3+9LmoB\nBMzL82jNUQDI0NfZQmXXXnst55xzDueddx7FxcXcd999/PnPf2b69OmUlJQwZswYLr/8cuLxOOC/\nEUYiEdavXw/AhRdeyOWXX86cOXPwPI8ZM2a0zYdI51yAp556iiOPPJKSkhIuu+wyZs6cyb333tvl\n3+XGG2/k7/7u7xg1ahTnnXceNTU1APz1r38lEolw9913M27cOE4++eROjwEsXbqUCRMmUFpayle+\n8hU++uijtuc46KCDuOWWW5g0aRLRaLQv/gn2owAIWCw3RjyrVl1AkrEeffRRLrjgAmpqajj77LPJ\nycnh9ttvp6qqipdffplnnnmGO++8s+38jte+L1myhJ/+9KdUV1dz0EEHce2116Z97tatWzn77LO5\n9dZb2b59O4ceeiivv/56lzX/4he/4KmnnmLFihVs3LiRaDTKpZdeus85//d//8eaNWtYtmxZp8c+\n/PBDvvGNb3DHHXewbds2Zs2axRlnnNEWdgAPPPAAzzzzDDt27EjjFU2dAiBg+dn5OGuluq4p6FJk\nCDPrm6/+MHPmTObMmQNAXl4eU6dO5ZhjjsHMOOSQQ7jkkkt46aWX2s7v2Io466yzOProo8nKyuL8\n889n5cqVaZ+7bNkyjj76aE4//XSysrJYsGABw4d3vXbSnXfeyY033khZWRm5ublce+21PPTQQ20/\nNzOuv/568vPzycvL6/TY/fffz9y5czn++OPJysrimmuuoaamhldffbXtca644goOOOCAtsfoaxoD\nCJiZkes8KuvqgP75RxYJ80TDgw46aJ/v16xZw5VXXsmbb75JQ0MD8XicadOmdXn/srKyttuFhYXs\n7KY/tatzN23atF8d3Q0er1+/nq9+9atEIv5naOcckUiErVv37m81ZsyY/e7X/timTZsYN25c2/dm\nxtixY/n0009TqqEvqAUQAvl4VO7UekCSmTp208yfP5+JEyeydu1aampquO666/p9mYsDDjiADRs2\n7HOs/RtxRwcddBDPPfccVVVVVFVVUV1dTX19PaNGpb7F+YEHHrjPGIRzjo0bN+7zpt/bpR6SUQCE\nQH4kRnWDBgFEwL9yqLi4mIKCAlavXr1P/39/Of3003n77bdZtmwZ8Xic//iP/2D79q7X6Jo/fz4/\n/vGP20Jj69atPPHEE20/7yywOh6bN28ejz/+OMuXL6elpYWf//zneJ7Hscce20d/q+QUACFQmOWx\no1EtABnaUv00e+utt3L33XfjeR7f+973OOecc7p8nGSPmeq5o0aN4oEHHmDBggWMGDGCv/3tbxx9\n9NFd9r1feeWVnHrqqcyaNYvi4mJmzpzJG2+80e1zdTw2fvx47rnnHr773e8yatQonn32WR5//PG2\nyz37+9M/gIV5BUEzc2Gur69M/Pc5lG34F56747SgS5FBysy0Gmgfam1t5cADD+SRRx5hxowZQZfT\npa7+3RPHQ7EpvCQRy/Woa1ILQCRIzzzzDDU1NTQ1NXH99deTm5s7oN0xQVAAhEAsL8bO3RoDEAnS\nihUrOOywwxg9ejTPPfccjz76KDk5OUGX1a/UBRQCZy++ktdeOIC/3XdV0KXIIKUuoMykLqAhoKRQ\n20KKyMBTAIRASVGMXU4BICIDSwEQAsOjHk1oDEBEBpYCIARGeh67I2oBiMjAUgCEwAjPI55VS2tr\n0JWISCZRAITAsPwYkYJa6uuDrkRkcGpubiYWi7F58+agSxlUFAAh4OV5WEGdNoWRIaevt4TcY/r0\n6fzud79r+z43N5e6urp9VvuU5LQcdAj4+wJrVzAZetLdEnKwC/P2j51RCyAE9mwLqV3BZCjrbEvI\n1tZWbrjhBg4//HBGjRrFhRdeSG2tf0FEQ0MD5557LsOHD6ekpITp06dTU1PDVVddxeuvv863v/1t\nPM/jRz/6EU1NTUQiETZt2gTAueeey4IFCzjllFPwPI/jjjtun+Wely1bxuc+9zlKS0tZsGDBfi2K\nVGtcs2YNOTk53HXXXRx88MGcdtppnR4DeOSRRzjqqKMoLS3lpJNO4uOPP257jgMOOIBbb72VCRMm\nUFxc3HcvehIKgBCI5cVozamjrk4zOSWz3HzzzTz//PO88sorbNy4kZycHBYsWADAXXfdRTwe57PP\nPqOyspJf/epX5Obmcsstt3DMMcewePFiamtrufnmm4HOt3/82c9+RnV1NWVlZSxcuBCAzz77jHPO\nOYfbbruNbdu2ceCBB/LWW2+lVeMVV1zR9vN4PM5rr73GRx99xGOPPdbpsXfffZeLLrqIX//612zd\nupUvf/nLzJ07l9Z2V348+OCD/PGPf6SysrJvXtwUqAsoBLIj2URa86isbQCKgi5HhiC7rm+WFnYL\n+/ZDyp133sl9993H6NGjAX+T+AkTJrB48WJycnLYtm0bf/nLXzjqqKOYOnXqvrV0aE10tt7+5MmT\nATjvvPO44YYbAHjyySc59thjOfXUUwG46qqruOWWW9Ku8Te/+Q3gB88NN9xAfn5+2306HnvggQf4\n2te+xnHHHQfAT37yE26//XbeeOONtgXnfvjDH7Y9x0BRAIRETqvH1tpaFADSH/r6jbuvbNiwgTlz\n5rR9et/zJl5VVcXFF1/M5s2bOeuss6ivr+fCCy/k3/7t31JeJz/V7R/NrNPtG1OpESASiez3xt3x\nWMftHyORCGPGjBnQ7R87oy6gkMjDo7JOk8Eks4wdO5YXXnhhv60VS0tLyc3N5brrrmP16tUsX76c\nhx56iPvvvx/o3WYpHbd/dM51u/1jdzV2VUvHYx23f2xtbeXTTz8d0O0fO6MACIk8i1FVr1FgySzz\n58/n6quvZuPGjYC/teKTTz4JwB//+EdWr16Nc45oNEp2dnbb1TGjR49m7dq1PXrOM844g9dee42n\nn36aeDzOrbfeyo4dO3pUI6S2/ePZZ5/N0qVLWbFiBS0tLdx0002MGDFiv26tgaYACImCiEd1g1oA\nMnR19gn36quvZvbs2Zx44oltWyu+/fbbgL8p+9y5c/E8j0mTJnH66aczb948ABYsWMA999zD8OHD\nueaaa/Z7/O4+TZeVlbFkyRIuvfRSRo4cyaZNm5g4cWKX2z92V2NXz9Xx2MSJE1m8eDHf+c53GDVq\nFC+++CKPPfYYkUgkab39KbD9AMxsAXAx0Aq8C1zknGvucE5G7AcAMP76uXx+10X8/sYzgy5FBiHt\nB9Bz8XicsrIynnzySaZNmxZ0OWkZlPsBmNmBwKXAF5xzk/AHo8/p/l5DW1GOR622hRQZEE8//TS1\ntbXs2rWLRYsWUVRUFHh3TBCC7ALKAorMLBsoBDYFWEvgYrkxdjZrDEBkICxfvpxDDz2UsrIyXnzx\nRZYuXUp2duZdFBlkF9BlwE+BBuBZ59yFnZyTMV1AZ95xDR+/W8x7v/5x0KXIIKQuoMzU2y6gQCLP\nzIYBc4FxQA3wsJmd55zbby72okWL2m6Xl5dTXl4+QFUOrGH5Hg1xdQGJSPoqKiqoqKhI+36BtADM\n7CzgZOfcJYnvLwSmOed+0OG8jGkBLFjyK+5/fjWfLb4j6FJkEFILIDMNykFgYD3wRTPLN//6p1nA\n6oBqCYXSaIxdTmMAIjJwAukCcs69ZmYPA28DuxN//ncQtYTFiKhHs6kLSHpm3LhxgV1LLsFpv7xE\nTwQ27O2cuw64LqjnDxvtCyy98cknn+x37PvLvs/4keP5wbE/2P8OImgmcGiMGhajJUsBIH2nqrGK\n4QXDgy5DQkwBEBIjPQ+XU0dLS9CVyFBR2VhJaUFp0GVIiCkAQqI438PytTG89J2qxiqGF6oFIF1T\nAISEl+dBnraFlL5T2aAWgHRPARASRTlFuKxd1NTGgy5FhgiNAUgyCoCQMDMi8Shba9QEkN5rjjfT\n2NLotyxFuqAACJGcuMeWHboSSHqvurGakvwSzQ2QbikAQiTHeWzXtpDSByobKzUALEkpAEIklxiV\nO9UFJL2nAWBJhQIgRArMo6peLQDpPQ0ASyoUACFSEPHYoX2BpQ9oEpikQgEQIkXZHjt2KQCk99QC\nkFQoAEKkKCdGXZPGAKT3NAYgqVAAhIiX51HXrBaA9J6WgZBUKABCpDjPo75FASC9pzEASYUCIESG\nFWhfYOkbGgOQVCgAQqSkKMauVo0BSO9pIpikQgEQIsOjHk2oBSC9p0FgSYUCIESGxzyatC+w9AF1\nAUkqFAAhMrrY07aQ0muNuxuJuziFOYVBlyIhpwAIkZHFMeJZGgOQ3tnz6V8rgUoyCoAQKSvxaM1R\nC0B6R5eASqoUACEyMuZvC9ncHHQlMphpEpikSgEQIvk5eeCMypqmoEuRQUxXAEmqFAAhE2mJ8VmV\nuoGk53QFkKRKARAyWS3aFlJ6p7KxUgEgKVEAhEx23GNbrQJAeq6qsUpdQJISBUDI5BJje60uBZWe\nq2zQMhCSGgVAyOThUaltIaUXdBmopEoBEDIF5lGtAJBe0CCwpEoBEDIFWR47GhUA0nNqAUiqFAAh\nU5QTo2aXxgCk5zQRTFKlAAiZWI62hZSec85pIpikTAEQMl6ex04FgPRQ/e56crJyyM/OD7oUGQQU\nACFTnK99gaXnKhs0CUxSF1gAmFmxmT1kZqvN7H0zmxZULWEyrDBGo7aFlB7SJDBJR3aAz30b8Afn\n3NfNLBvQ7hVAaaFHo1MLQHpGewFLOgIJADPzgOOcc98CcM61gDbDhcS2kHoppIc0ACzpCKoL6FBg\nu5n91szeMrP/NrOCgGoJlZExj92mLiDpGU0Ck3QEFQDZwBeAO5xzXwAagGsCqiVURhTHtC+w9Jgm\ngUk6ghoD2AhscM69kfj+YeDqzk5ctGhR2+3y8nLKy8v7u7ZAlQ3ziGcrAKRnqhqrGBMbE3QZMsAq\nKiqoqKhI+37mnOv7alJ5YrOXgEuccx+Z2UKg0Dl3dYdzXFD1BWV7ZZyRt+fSuqhFm3pL2r756Dcp\nH1fORUdfFHQpEiAzwzmX9A0kyKuALgPuM7McYC2g31ig2MuClgJ2NtcTy4sGXY4MMloGQtIRWAA4\n51YBxwT1/GGVkwM0x9haU0tslAJA0qOJYJIOzQQOochubQspPaOJYJIOBUAIZbd4bK1RAEj6NBFM\n0qEACKEc57GtRnMBJD2trpXqxmpK8kuCLkUGCQVACOW6GNt3qgUg6altqqUot4icrJygS5FBQgEQ\nQnnmUaUAkDRpGQhJlwIghAojHlUNCgBJj5aBkHQpAEKoICtGTaPGACQ9WgZC0qUACKFotkdtk1oA\nkh5NApN0KQBCKJanAJD0aRKYpEsBEEJenraFlPRpEpikq9sAMLML2t2e0eFnP+ivojJdcX6MhhaN\nAUh6KhvVApD0JGsB/LDd7f/s8LN/7uNaJKGkyKOhVS0ASY9aAJKuZAFgXdzu7HvpI6VFHk3aF1jS\npGUgJF3JAsB1cbuz76WPjIh6NGlbSEmTJoJJupItB/15M3sH/9P+4YnbJL4/rF8ry2AjvRi7TS0A\nSY8mgkm6kgXA3w9IFbKPUcWe9gWWtGkimKSr2wBwzq1r/72ZDQe+DKx3zr3Zn4VlshHFhbhIEy2t\nLWRHgty0TQaLeGucuqY6huUPC7oUGUSSXQb6pJlNSNw+AHgP/+qf/zWzKwagvozkeYbt9qhr0jiA\npKZ6VzXF+cVkRbKCLkUGkWSDwIc6595L3L4IeM4591VgGroMtN9Eo0BTTLOBJWW6BFR6IlkA7G53\nexbwBwDnXB3Q2l9FZbqiInC7PHY0KgAkNVoGQnoiWQfzBjO7FNgIfAF4GsDMCgDtOtFPsrL8fYG3\n1tbCAUFXI4OBWgDSE8laABcDRwHfAs52zu1IHP8i8Nt+rCvjZce1LaSkTpPApCeSXQW0FfhuJ8df\nBF7sr6IEclpjbK9TF5CkprKhktJ8tQAkPd0GgJk93t3PnXNn9G05skcengJAUqa9AKQnko0BTAc2\nAEuAV9H6PwMmzzyq6hUAkprKxkrGjxwfdBkyyCQLgDJgNnAucB6wDFjinHu/vwvLdIURjx3aFlJS\npGUgpCe6HQR2zsWdc087576JP/D7MVChvQD6X2F2jBpdBiop0iCw9ETSdQbMLA84Db8VcAhwO7C0\nf8uSaI5HTdOaoMuQQUKXgUpPJBsEvheYgD8B7Lp2s4Kln8VyPbY1qwUgqdFEMOmJZPMALgCOAC4H\nXjGz2sRXnZnWK+5Pxfke9doWUlKkFoD0RLJ5ANo0PiDDCmI07FbGSnLN8WYaWxrx8rygS5FBRm/w\nIVVS6NGofYElBVWNVZTkl2Cmq7QlPQqAkCqNejQ5dQFJcpoEJj2lAAipEbEYzRpmkRRoL2DpKQVA\nSI30YuyO1OKcC7oUCTlNApOeUgCEVImXBy6LXS27gi5FQk6TwKSnAg0AM4uY2VvJFp3LRLGYvydA\nXbPGAaR7VY1VWglUeiToFsDlwAcB1xBK0SjYbm0LKclVNqgFID0TWACY2VhgDnBXUDWEWTQK7PIU\nAJKUJoFJTwXZAvgl8CNAo5ydiMWgtVEBIMlVNmoZCOmZQALAzE4DtjjnVuLvMaAZLB0UFPgBoCWh\nJZnKRl0GKj2TdDXQfjIDOMPM5gAFQMzM7nXOfaPjiYsWLWq7XV5eTnl5+UDVGKhIBLJbY2yvVQtA\nuqeJYFJRUUFFRUXa97OgrzM3s+OBKzvbXtLMXND1Bano7O/yr9+ezI9nfy/oUiTExv5iLK9c/AoH\nFx8cdCkSEmaGcy5pz0rQVwFJN/LwqNS2kJKEJoJJTwXVBdTGOfcS8FLQdYRRvnlU12sMQLrWuLuR\nVtdKYU5h0KXIIKQWQIgVZsWo1raQ0o09l4BqJVDpCQVAiBVme9TsUgBI17QMhPSGAiDEYjkeO5vU\nBSRd0yQw6Q0FQIh5eR512hdYuqG9gKU3FAAh5uXHqG9RAEjX1AKQ3lAAhNiwAo8GbQsp3dAyENIb\nCoAQGx712NWqMQDpmnYDk95QAIRYaZFHE2oBSNe0DIT0hgIgxIbHorRYPa2uNehSJKTUBSS9oQAI\nsWIvQlZrITubdwZdioSUBoGlNxQAIRaNQlaLR53mAkgXNBFMekMBEGLaFlKSUQtAekMBEGKxGFiz\ndgWTzjnndBWQ9IoCIMSiUXDaF1i6UL+7npysHPKz84MuRQYpBUCIRaP+tpB1zRoDkP1pGQjpLQVA\niMVi0FKvMQDpnPYClt5SAIRYXp7fBVSlXcGkE5oEJr2lAAgxM8h1HpU7FQCyP3UBSW8pAEIuzzyq\ntC2kdEKXgEpvKQBCriASo1pdQNIJLQMhvaUACLnCLI8d2hdYOqEWgPSWAiDkojketVoKQtpxzrG2\nei2rt6/WILD0SnbQBUj3/ABQCyBT7WrZxXtb32PV5lWs3LySlVtW8s6Wd/DyPCaPnsyXDvpS0CXK\nIKYACDkvL8YmbQuZEbbWb93njX7V5lX8tfqvHFF6BFPKpjClbApnfv5MJpdNZkThiKDLlSFAARBy\nwwo8GhQAQ0q8Nc5fqv7S9ma/aov/Z2NLI5NHT2ZK2RRmHzabq6ZfxfiR48nLzgu6ZBmiFAAhN6zQ\nozGuMYDBamfzTt7d8q7/qT7xZv/e1vcYVTSq7VP9/KnzmVI2hYOLD8bMgi5ZMogCIORKizx2ObUA\nws45x6a6Tfu80a/cvJKNtRsZP3I8U8qmMHn0ZC6YdAGTRk+iOL846JJFFABhVxLNJ16/m+Z4M7lZ\nuUGXI8Du+G4+3P7hPm/0KzevJGKRtk/1c4+cy8LjF3LkiCPJjui/mYSTfjNDLhYzcnf6u4Lpkr+B\nt2PXjv366j/c/iEHFx/c9qn+yulXMqVsCmXRMnXhyKCiAAi5aBSy4/6S0AqA/uOc45Mdn+z3qX57\nw3YmjZ7ElLIpTBszjflT5zNh1ASKcouCLlmk1xQAIReLQaRFcwH60q6WXby/9f193uzf2fIO0dxo\n26f68yeez82zb+bw0sOJmOZLytCkAAi5aBSsWXsC9NS2+m37fapvf2395NGTmXvkXF1bLxlJARBy\n0ShYk1oAycRb43xc9fE+b/Srtqyivrm+7Y1+1qGzuHL6lbq2XiRBARBysRi07vIHgcW359r69m/2\n7a+tnzx6MvOnzmdy2WTGFY/TwKxIF8w5F3QNXTIzF+b6BsL69TD+/81n7MyXKIuW4XDseU063gZ/\nMLOr20HeJ9n907lPXVMd40eOb5s1O6Vsiq6tF2nHzHDOJf3kowAIuaoqOHTiFh57eTUAhv9vamb7\n3N7zs65up3qfzs7ryX36qp7O7jO8YDg5WTk9eTlFMkKoA8DMxgL3AqOBVuB/nHO3d3JexgdAczMU\nFfl/qidDRFIR9gAoA8qccyvNLAq8Ccx1zn3Y4byMDwDwN4evqYH8/KArEZHBINUACOQCZ+fcZufc\nysTtncBqYEwQtQwG0Sjs3Bl0FSIy1AQ+w8XMDgGmAK8GW0l4xWIKABHpe4FeBpro/nkYuDzREtjP\nokWL2m6Xl5dTXl4+ILWFSTQKdboKVES6UFFRQUVFRdr3C+wqIDPLBp4EnnLO3dbFORoDAL74Rfjl\nL2H69KArEZHBINRjAAm/AT7o6s1f9lIXkIj0h0ACwMxmAOcDJ5rZ22b2lpmdEkQtg4G6gESkPwQy\nBuCcexnICuK5B6NoFNat82cFiz8fIhLp2ZfmUojspZnAg8Att8Dt+02Ty1zOQWtrz74AsrJ6HiDd\nfQXxuFOnwplnwnBtFSHthHoiWKoUANKXnOtdeHT2FY/37eOl85i7d8Py5fDss/ClL8G8eX4YlJQE\n/UpL0BQAIhli50548kl48EF4/nk47jg/DObOhWHDgq5OgqAAEMlAdXXwxBN+GLzwAhx/vB8GZ5wB\nxVosNWMoAEQyXG0tPP64HwYVFXDCCX4YfPWr4HlBVyf9SQEgIm127NgbBsuXw6xZfhicfro/z0SG\nFgWAiHSquhoee8wPgxUrYPZsPwxOO82/5FgGPwWAiCRVVQWPPuqHwZ/+BCed5IfBnDn+PhQyOCkA\nRCQtlZWwdKkfBq++Cqec4ofBqadCYWHQ1Uk6FAAi0mPbtu0Ngzfe8ENg3jw/FAoKgq5OklEAiEif\n2LoVfv97PwzeessfK5g3D04+WbvUhZUCQET63ObNe8Ng1Sr/KqJ58/yxg7y8oKuTPRQAItKvPvsM\nHnnED4N33/Unm82b519VlJsbdHWZTQEgIgPm00/3hsEHH/jLUMyb5883UBgMPAWAiARi40Z4+GE/\nDNas8ReomzcPTjwRcnKCri4zKABEJHDr1+8Ng48/hn/8Rz8MTjgBsgPdkXxoUwCISKisW7c3DNau\nha99zQ+D449XGPQ1BYCIhNbf/rY3DNatg3/6Jz8MvvxlfwMc6R0FgIgMCmvXwkMP+WHw6ad7w2Dm\nTIVBTykARGTQ+fjjvWGweTOcdZYfBjNm+FtgSmoUACIyqH300d4w2L59bxhMn64wSEYBICJDxocf\n7g2D6mr4+tf9MJg2TWHQGQWAiAxJH3zgB8GDD/r7IX/96/7X+PH+fgYKBAWAiAxxzsH77/stg0ce\n8a8mamjwl66OxfxtL5P92d3PiooGb5goAEQk47S2+q2Cujp/T+SOf3Z2rKs/Gxv9EOhNiLQPE0v6\ndtx3FAAiIr0Qj+8bJumER8f77Nrld0/1RZgUFiYPk1QDQPPvREQ6kZUFxcX+V2+1tPhhkiw0Kivh\nk0+6D5ymJj8MuguRVCkARET6WXY2DBvmf/VWS0vylkeq1AUkIjLEpNoFNEjHuEVEpLcUACIiGUoB\nICKSoRQAIiIZSgEgIpKhFAAiIhkqsAAws1PM7EMz+8jMrg6qDhGRTBVIAJhZBPgVcDJwFHCumX0+\niFoGi4qKiqBLCA29FnvptdhLr0X6gmoBHAv8xTm3zjm3G7gfmBtQLYOCfrn30muxl16LvfRapC+o\nABgDbGg9fApWAAAGDElEQVT3/cbEMRERGSAaBBYRyVCBrAVkZl8EFjnnTkl8fw3gnHM/63CeFgIS\nEemB0O4HYGZZwBpgFvAZ8BpwrnNu9YAXIyKSoQJZDto5FzezHwDP4ndDLdabv4jIwAr1ctAiItJ/\nQjkIbGaLzWyLmb0TdC1BMrOxZvaCmb1vZu+a2WVB1xQUM8szs1fN7O3Ea7Ew6JqCZmYRM3vLzB4P\nupYgmdknZrYq8bvxWtD1BMnMis3sITNbnXjfmNbt+WFsAZjZTGAncK9zblLQ9QTFzMqAMufcSjOL\nAm8Cc51zHwZcWiDMrNA515AYQ3oZuMw5l7H/4c1sATAV8JxzZwRdT1DMbC0w1TlXHXQtQTOzu4GX\nnHO/NbNsoNA51+UeYaFsATjnVgAZ/4/pnNvsnFuZuL0TWE0Gz5dwzjUkbubhj1+F79PLADGzscAc\n4K6gawkBI6TvZQPJzDzgOOfcbwGccy3dvfmDXrRBw8wOAaYArwZbSXASXR5vA5uB55xzrwddU4B+\nCfyIDA7BdhzwnJm9bmaXBF1MgA4FtpvZbxNdg/9tZgXd3UEBMAgkun8eBi5PtAQyknOu1Tl3NDAW\nmGZm44OuKQhmdhqwJdE6tMRXJpvhnPsCfovoXxJdyJkoG/gCcEfi9WgArunuDgqAkEv04z0M/K9z\n7rGg6wmDRLP2ReCUoGsJyAzgjETf9xLgBDO7N+CaAuOc+yzx5zZgKf5aY5loI7DBOfdG4vuH8QOh\nS2EOAH2y8f0G+MA5d1vQhQTJzEaYWXHidgEwG8jIwXDn3E+ccwc75w4DzgFecM59I+i6gmBmhYkW\nMmZWBJwEvBdsVcFwzm0BNpjZ5xKHZgEfdHefQCaCJWNmvwPKgeFmth5YuGdgI5OY2QzgfODdRN+3\nA37inHs62MoCcQBwT2Ip8QjwgHPuDwHXJMEbDSxNLBuTDdznnHs24JqCdBlwn5nlAGuBi7o7OZSX\ngYqISP8LcxeQiIj0IwWAiEiGUgCIiGQoBYCISIZSAIiIZCgFgIhIhlIAyJBgZq3tZ8OaWZaZbeuP\npZLN7EUz63aGpchgoACQoaIemGBmeYnvZwMbAqwnqcSy1iKBUQDIUPIH4LTE7XPx18kBwMyOMbNX\nzOxNM1thZkckjl9hZosTtycmNpvJb/+gZpZvZksSG2z8Hshv97PZicd9w8weMLPCxPE5iU05Xjez\n28zsicTxhWZ2r5mtAO5NrHD688RmNyvbr2ZpZleZ2WuJ4xm/AY70PQWADBUOuB84N9EKmMS+S2ev\nBmY656YCC4GbEsdvAw43szPx1126xDm3q8Njfw+od84dlbjvPwCY2XDgX4FZzrl/wN+w54eJ5/81\ncLJz7hhgJPsu2/z3wInOufOBi4Edzrlp+IuYfcfMxpnZbOAI59yxwNHAP2TwKpfST0K5FpBITzjn\n3kvsm3AusIx9FxMchv+J+wj8N+PsxH2cmV0EvAP82jn3504e+sv4QYFz7l0zW5U4/kVgPPCymRmQ\nA/wJ+DzwV+fc+sR5S4D269Q/7pxrTtw+CZhoZl9PfO8BRySOzzaztxJ/j6LE8RVpvSgi3VAAyFDz\nOHAz/mKCI9odvwF/1cyvmdk4/OWk9/gcUAccmOJzWLs/n018kt/7Q7PJdL+SbX2Hx7rUOfdch8c4\nBbjJOfc/KdYkkjZ1AclQsecN9zfAdc659zv8vBj4NHG7bYXExBLTt+F/yh9uZv/UyWMvx1+VFTOb\ngN+9BPBnYIaZHZ74WWGihbEGONTMDk6cd3Y3dT8DfD+x7wNmdkRiHOEZ4J8TSxxjZgea2cjuXgCR\ndCkAZKhwAM65T51zv+rk5z8H/t3M3mTf3/tfAP/pnPsY+DZwk5mN6HDf/wKiZvY+sAh4I/Fc24Fv\nAUsS3UKvAEcmxhC+DzxjZq8DtUBNF3Xfhb9m+1tm9i7+2EFWokXwO+BPZvYO8BAQTfXFEEmFloMW\n6QdmVuScq0/cvgP4KNM39ZHwUQtApH9cYmZvJ1oNHnBn0AWJdKQWgIhIhlILQEQkQykAREQylAJA\nRCRDKQBERDKUAkBEJEMpAEREMtT/B73ki9MfeUKRAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f9debea6950>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"import matplotlib.pyplot as plt\n",
"\n",
"max_order=9\n",
"mse_train_list=[]\n",
"mse_test_list=[]\n",
"for order in range(1,max_order+1):\n",
" x,y=gen_poly_data(10,True,3)\n",
" x_train = gen_poly_features(x,order)\n",
" y_train = y\n",
"\n",
" xt,yt=gen_poly_data(100,False)\n",
" x_test = gen_poly_features(xt,order)\n",
" y_test = yt\n",
"\n",
" ws=compute_weights(x_train,y_train)\n",
" mse,_=compute_mse(ws,x_train,y_train)\n",
" print(mse)\n",
" mse_train_list.append(mse)\n",
"\n",
" mse_test,ye=compute_mse(ws,x_test,y_test)\n",
" print(mse_test)\n",
" mse_test_list.append(mse_test)\n",
"\n",
" line1,=plt.plot(xt,ye,label=\"Estimated curve\")\n",
" plt.plot(x,y,'x')\n",
" line2,=plt.plot(xt,yt,'r',label=\"Original curve\")\n",
" plt.legend(handles=[line1,line2])\n",
"# plt.savefig('curve_fit'+str(order)+'.pdf') # save figure\n",
" plt.cla()\n",
"\n",
"p_order=6\n",
"line1,=plt.plot(range(1,p_order+1),mse_train_list[0:p_order],label=\"Training error\")\n",
"line2,=plt.plot(range(1,p_order+1),mse_test_list[0:p_order],label=\"Testing error\")\n",
"plt.legend(handles=[line1,line2])\n",
"plt.xlabel('Max degree')\n",
"plt.ylabel('MSE')\n",
"plt.ylim([0,10])\n",
"# plt.savefig('curve_fit_err.pdf') # save figure"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"#### Incorporate regularization"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/phsamuel/anaconda3/envs/conda2/lib/python2.7/site-packages/sklearn/linear_model/coordinate_descent.py:484: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
" ConvergenceWarning)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Wd0VFUXgOH3THovlNBDL6JUBUGRAIrSi4KABRURBUUQ\nlaJowAKIqOhnQ5Ci9N6lCBGkSAelhJ4AKRDSSSOZ8/24SQghlCQTJiH7WeuumdvPjLjnZN9TlNYa\nIYQQ9z6TtQsghBDi7pCAL4QQxYQEfCGEKCYk4AshRDEhAV8IIYoJCfhCCFFMSMAXhY5S6oxSqnUB\nXLelUupcLo7vq5TaaulyCGEtEvBFcZPbjid3taOKUspLKTVfKRWhlLqolPpNKeV6N8sg7l0S8IUo\nXD4DPABfoBpQBvC3ZoHEvUMCvijUlFIPKaW2K6WilFIXlFLfKaVss+w3K6XeUEodV0rFKKXGKqWq\nKqW2KaWilVLzsh5vnKJGKqUuKaVOK6X6ZNnhrZRakX6dnRgBN+uJ3yilgtP371ZKPVoAH7kysExr\nfUVrHQcsBeoWwH1EMSQBXxR2acAQwBtoBrQGBmY7pi3QEHgYeB/4GegDVAQeAHpnObZM+rXKAS8B\nU5RSNdL3/QAkAD5AP+CVbPfZBdQDvIA5wEKllH1OhVZKDU//kYpMf836PvIWn/d7oJNSylMp5QU8\nDay5xfFC3DEJ+KJQ01rv01rv0oZgYArQMtthE9JrxEeB/4D1Wuug9BryWowfg8xLAqO11le11luA\n1UBPpZQJ6J6+L0lrfRiYma0sc7TW0Vprs9b6a8ABqHWTck/QWntprb3TX7O+977FR94H2AOXgUtA\nKvDjnXxXQtyOBHxRqCmlaiilViqlQpVS0Rg57pLZDruY5X0iEJ5tPetDzyitdVKW9SCM2n4pwBY4\nn21f1rK8q5Q6klFjB9xzKEt+LQQCAZf0658GZlv4HqKYkoAvCrsfgaNANa21J/ABoPJxPS+llFOW\n9UpACNdq0xWz7QNAKdUCeA94JqPGDsTerCzpzwnilFKx2ZY4pVTsLcpXH/g5/a+MBOAnoF1ePqgQ\n2UnAF4WdGxCrtU5QStUG3sjn9RQwRilllx7EOwALtNZmYDHgr5RyUkrdB/TNcp4rcBW4rJSyV0p9\nlF62HGmtx2mt3bTW7tkWN621+y3Ktwt4VSnlmP7DNAA4lL+PLIRBAr4ojLK2fX8XeC69VvwzMO8W\nx+a0nl0oEIVRq/8NGKC1PpG+7y2MIB4K/Jq+ZFiXvhwHzmA83L3jTly58ApQBSO1dA6j1U7fW50g\nxJ1SlpgARSk1FKNVgxn4F3gZIwc5H6M98Vmgp9Y6Jt83E0IIkSf5ruErpcph1Iwaaa3rYTz46g2M\nADZqrWsBm4CR+b2XEEKIvLNUSscGcEnv4OIEXAC6cK1Z20ygq4XuJYQQIg/yHfC11iHAJCAYI9DH\naK03Aj5a6/D0Y8KA0vm9lxBCiLyzRErHE6M274vRntlFKfUcuX+YJoQQogDZ3v6Q23ocOK21jgRQ\nSi0FmgPhSikfrXW4UqoM13eOyaSUkh8CIYTIA611rvqkWCKHHww8nN5uWAFtgCPACoyxSsBoVrb8\nZhfQWsuiNR9//LHVy1BYFvku5LuQ7+LWS17ku4avtd6llFoE7MfomLIfY7wTN2CBUuoVjC7qPfN7\nLyGEEHlniZQOWusxwJhsmyMx0j1CCCEKAelpW4j4+flZuwiFhnwX18h3cY18F/ljkZ62+SqAUtra\nZRBCiKJGKYXO5UNbi6R0hBA5q1y5MkFBQbc/UIib8PX15ezZsxa5ltTwhShA6bUwaxdDFGE3+zeU\nlxq+5PCFEKKYkIAvhBDFhAR8IYQoJiTgCyFEMSEBX4hirEqVKmzatMnaxRB3iQR8IYQoJiTgCyGu\nEx0dTadOnShdujQlSpSgU6dOXLhwIXP/jBkzqFatGu7u7lSrVo25c+cCcOrUKfz8/PD09KR06dL0\n7t0785zt27fTpEkTvLy8aNq0KTt27Ljrn0tIwBdCZGM2m3nllVc4d+4cwcHBODs78+abbwKQkJDA\n22+/zbp164iNjWX79u00aNAAgNGjR/Pkk08SHR3N+fPneeuttwCIioqiY8eODBkyhMuXLzN06FA6\ndOhAVFSU1T5jcSUBXwgrU8oyi6V4e3vTrVs3HBwccHFxYeTIkWzZsiVzv42NDf/++y9JSUn4+PhQ\np04dAOzs7AgKCuLChQvY29vTvHlzAFavXk3NmjXp06cPJpOJXr16Ubt2bVauXGm5Qos7IgFfCCvT\n2jKLpSQmJjJgwAAqV66Mp6cnLVu2JDo6Gq01zs7OzJ8/nx9//JGyZcvSqVMnAgMDAZg4cSJms5km\nTZrwwAMPMH36dABCQkLw9fW97h6+vr7XpYnE3SEBXwhxnUmTJnHixAl2795NdHR0Zu0+o3v/E088\nwfr16wkLC6NWrVr0798fgNKlSzNlyhQuXLjATz/9xMCBAzl9+jTlypW7YSyY4OBgypcvf1c/V1F3\n4vIJZh+azadbPuXVFa/m6RoS8IUo5lJSUkhOTiY5OZmkpCSioqJwcnLC3d2dyMhI/P39M4+9ePEi\nK1asICEhATs7O1xdXbGxsQFg0aJFmbV2T09PTCYTJpOJ9u3bc+LECebNm0daWhrz58/n6NGjdOzY\n0Roft0g6GHaQ5r82Z+XxlVxJuUKT8k3ydiELTLNVE2OWq33przHAYMALWA8EAusAj5ucr4W4VxX2\nf9+VK1fWJpNJm0wmrZTSJpNJ9+vXT7dq1Uq7urrqWrVq6SlTpmiTyaTT0tJ0aGiobtmypfb09NRe\nXl66VatW+ujRo1prrd9//31dvnx57ebmpqtXr66nTp2aeZ9t27bpxo0ba09PT/3ggw/q7du3W+sj\nFzmArvhVRT3v33k3bNe5jNcWHS1TKWUCzgNNgTeBy1rrL5RSwwEvrfWIHM7RliyDEIWJjJYp8ksp\nxad/fcoHj31ww3ady9EyLR3w2wKjtdYtlFLHgJZa63ClVBkgQGtdO4dzJOCLe5YEfJFfSinMZjMq\nW1OswjA88rPAnPT3PlrrcACtdRhQ2sL3EkKIYiF7sM8ri814pZSyAzoDw9M3Za/W3LSak/WhkJ+f\nn8xbKYQQ2QQEBBAQEJCva1gspaOU6gwM1Fo/lb5+FPDLktLZrLWuk8N5ktIR9yxJ6Yj8KqwzXvUG\n5mZZXwG8lP6+L7DcgvcSQgiRSxap4SulnIEgoKrWOi59mzewAKiYvq+n1jo6h3Olhi/uWVLDF/ll\nyRq+TGIuRAGSgC/yq7CmdIQQQhRiEvCFEDf1xhtv8Nlnn910v8lk4vTp03exRIXX7b6rwkBSOkIU\noMKe0qlcuTIXL17E1tYWV1dXnnzySb7//nucnZ3v6HwbGxtOnDhB1apVC7ikxZekdIQQFqGUYvXq\n1cTGxnLgwAH279/PuHHj7vj8u/VjZun7FOYf4YIkAV+IYi4j+JUuXZonn3ySAwcOZO57+eWX+eij\njzLXJ06cSLly5ahQoQLTp0+/rgdoZGQknTp1wsPDg6ZNmzJ69GhatGiRuf/YsWO0bduWEiVKUKdO\nHRYuXHjTMrVq1YoPP/yQRx99FBcXF86cOUNsbCz9+vWjXLlyVKxYkdGjR2eW3Ww2M2zYMEqVKkW1\natX4/vvvMZlMmM3mPF3vVtM1Dh06FB8fHzw8PKhfvz5HjhzJ8bv65ZdfqFGjBiVLlqRr166EhoZm\n7jOZTPz888/UrFkTb2/vzBnFCpoEfCEEAOfPn2ft2rXUqFEjx/1//PEHX331FX/++ScnTpxg48aN\n1+0fOHAgbm5uXLx4kRkzZjBz5szMH4SEhATatm3L888/T0REBPPmzWPQoEEcO3bspuX5/fffmTp1\nKnFxcVSqVIm+ffvi4ODA6dOn2b9/Pxs2bGDq1KkATJkyhXXr1nHo0CH27dvHsmXLbhiOIDfXu9l0\njevXr+fvv//m5MmTxMTEsGDBAkqUKHFD2Tdt2sSoUaNYtGgRoaGhVKpUiV69el13zOrVq9m7dy8H\nDx5kwYIFrF+//lb/eSzCYkMrCCHyRo2xzDgp+uO8pSm6du0KQHx8PG3atLluqJOsFi5cyMsvv5w5\npaG/v3/mBOZms5klS5Zw5MgRHBwcqFOnDn379uWvv/4CYNWqVVSpUoUXX3wRgPr169O9e3cWLlzI\n6NGjc7zfSy+9RO3axniLERERrF27lpiYGBwcHHB0dGTIkCH88ssv9O/fn4ULF/L2229TtmxZAEaM\nGMGmTZvyfL2s0zWWL18+c7pGOzs74uLiOHLkCE2aNKFWrVo5ln3OnDn069eP+vXrAzBu3Di8vLwI\nDg6mUqVKAIwcORI3Nzfc3Nxo1aoVBw4coG3btrf6T5VvEvCFsLK8BmpLWb58Oa1atWLr1q306dOH\niIgI3N3dbzguJCSEBx98MHM967SFly5dIi0tjQoVKmRuq1ixYub7oKAgdu7cibe3N2CkkdLS0njh\nhRduWq7s51+9ejUzoGeM754RPENCQq47Puv7vFxv4sSJfPjhhzRp0gRvb2/eeecdXn75ZVq1asWb\nb77JoEGDCA4Opnv37nz55Ze4urre8F01btw4c93FxYUSJUpw4cKFzHv4+Phk7nd2diY+Pv6m34Wl\nSMAXopjLyFu3aNGCvn37MmzYMJYuXXrDcWXLluXcuXOZ60FBQZlpk1KlSmFra8v58+epXr06wHXH\nVqxYET8/P9atW3fH5cqakqlYsSKOjo5cvnw5x5Ejy5Yty/nz5zPXg4OD83W9jOkaAbZt28bjjz9O\ny5YtqVq1Km+++SZvvvkmERER9OjRg4kTJzJmzJjrzi9XrhxBQUGZ61euXOHy5cvX/SBag+TwhRCZ\nhgwZwoYNG/j3339v2NezZ09mzJjB0aNHSUhIYOzYsZn7TCYT3bt3x9/fn8TERI4dO8asWbMy93fs\n2JHjx4/z+++/k5qaytWrV9mzZ88tc/hZlSlThrZt2zJ06FDi4uLQWnP69OnM+XZ79uzJ5MmTCQkJ\nITo6mi+++CJf17vZdI179uxh165dpKam4uTkhKOjIybTjWG0d+/eTJ8+nUOHDpGcnMyoUaN4+OGH\nc/zL426SgC9EMZa9dluyZEn69u17XTDP8NRTTzFkyBBat25NzZo1adOmzXX7v/vuO6Kjoylbtix9\n+/alT58+ODg4AODq6sr69euZN28e5cqVo1y5cowYMYKUlJQ7KhfArFmzSElJ4b777sPb25sePXoQ\nFhYGQP/+/Wnbti316tWjcePGdOjQAVtb28xgnNvr7d69m6ZNm+Lu7k7Xrl359ttvqVy5MrGxsfTv\n3x9vb2+qVKlCyZIlee+99264dps2bfjkk0/o3r075cuX58yZM8ybN++mn89S493fjnS8EqIAFfaO\nVwVpxIgRhIeHM3369Lt+7z/++IM33niDM2fO3PV7W5p0vBJCFDqBgYGZqaBdu3Yxbdo0unfvflfu\nnZSUxNq1a0lLS+PChQuMGTPmrt27KJEavhAFqDjV8Pfs2UPv3r0JDQ3Fx8eHAQMG8P7779+Veycm\nJtKyZUsCAwNxcnKiY8eOfPPNNze0nimKZHhkIYqI4hTwRcGQlI4QQohcs0jAV0p5KKUWKqWOKqUO\nK6WaKqW8lFLrlVKBSql1SikPS9xLCCFE3liqhj8ZWJM+SXl94BgwAtiota4FbAJGWuheQggh8iDf\nOXyllDuwX2tdLdv2Y0BLrXW4UqoMEKC1rp3D+ZLDF/csyeGL/CpsOfwqQIRSarpSap9Sakr6pOY+\nWutwAK11GFDaAvcSQgiRR5YYS8cWaAQM0lrvUUp9jZHOyf6TdNNqTtbR+fz8/PDz87NAsYQQBWXc\nuHGcOXMmc7wZSx17OyaTiZMnTxbLGbYCAgIICAjI1zUskdLxAXZoraumrz+KEfCrAX5ZUjqb03P8\n2c+XlI64ZxWFlM6MGTP46quvOHXqFB4eHnTt2pVx48bh4VH42lkUxykVC1VKJz1tc04pVTN9Uxvg\nMLACeCl9W19geX7vJYSwrEmTJjFy5EgmTZpEbGwsO3fuJCgoiCeeeILU1NQcz0lLS7vLpbzmbvx4\nZsySdS+yVCudwcBspdQBjFY6nwMTgCeUUoEYPwLjLXQvIYQFxMXF4e/vz//+9z+eeOIJbGxsqFSp\nEgsWLODs2bP8/vvvAIwZM4YePXrwwgsv4OnpycyZMxkzZsx1Y9nPmjWLypUrU6pUKT799FOqVKmS\nOQFJ1mODgoIwmUzMmjULX19fSpcuzeeff555nd27d9O8eXO8vLwoX748b7311k1/eLKLiorilVde\noXz58pQoUSJzaIWZM2deN9UiGKmh06dPA8bUhAMHDqRDhw64ubnx5ZdfUrZs2et+XJYuXZo5mYnW\nmvHjx1O9enVKlSpFr169iI6OztV3by0WCfha64Na64e01g201t211jFa60it9eNa61pa67Za66Lx\njQhRTGzfvp3k5GS6det23XYXFxfat2/Phg0bMretWLGCnj17Eh0dTZ8+fYBrIzweOXKEQYMGMXfu\nXEJDQ4mJiSEkJOS6a2YfDXLbtm2Z0ySOHTuWwMBAwEjZfPPNN0RGRrJjxw42bdrEDz/8cEef5/nn\nnycxMZGjR49y8eJFhg4detP7Z1+fO3cuo0ePJi4ujrfffhtXV9frZsyaO3cuzz//PADffvstK1as\nYOvWrYSEhODl5cXAgQPvqIzWJj1thbA2pSyz5FJERAQlS5bMcTz3smXLEhERkbnerFkzOnXqBICj\no+N1xy5evJjOnTvTrFkzbG1tcxxa+fqPq/D398fe3p569epRv359Dh48CECjRo1o0qQJSikqVarE\na6+9ljlN4q2EhYWxbt06fv75Z9zd3bGxsbmhVp9V9tRQly5dePjhhwFwcHCgV69ezJkzBzD+Elqz\nZk3mROY///wzn332GWXLlsXOzo6PPvqIRYsWFYlUkMx4JYS1WemhbsmSJYmIiMBsNt8Q9ENDQylZ\nsmTm+q0m7sg+vaCTk1OOE3tndbPp/U6cOME777zDnj17SExMJDU19bqpAm/m3LlzeHt75zg1453I\n/vn69OnDI488wk8//cSSJUto3Lhx5mxVQUFBdOvWLfM701pjZ2dHeHh45pSJhZXU8IUoppo1a4aD\ngwNLliy5bnt8fDxr167l8ccfz9x2qwk6sk8vmJiYyOXLl/NUpjfeeIM6depw6tQpoqOj+eyzz+7o\nQW3FihWJjIwkNjb2hn0uLi4kJCRkrmdMcpJV9s9Xp04dfH19WbNmDXPnzs1MYwFUqlSJtWvXEhkZ\nSWRkJFFRUVy5cqXQB3uQgC9EseXu7s5HH33EW2+9xbp160hNTeXs2bM8++yzVKpUKTNnfTvPPPMM\nK1euZOfOnVy9evW6fjU5uVUAj4uLw93dHWdnZ44dO8aPP/54R2UoU6YM7dq1Y+DAgURHR5OamsrW\nrVsBqF+/PocPH86cbnDMmDF3NMNUnz59mDx5Mlu3bqVHjx6Z2wcMGMCoUaMy5829dOkSK1asuKNy\nWpsEfCGKsffee4/PP/+cd999Fw8PD5o1a4avry8bN27Ezs7ujq5x33338d133/Hss89Srlw53N3d\nKV26dOb0htnd6gHql19+yezZs3F3d2fAgAH06tXrludm9dtvv2Fra0vt2rXx8fFh8uTJANSoUYOP\nPvqINm3aULNmzVvm9rPq1asXW7ZsoU2bNnh7e2duf/vtt+nSpQtt27bFw8OD5s2bs2vXrju6prXJ\nePhCFKCi0PHK0q5cuYKnpycnT57E19fX2sUp8gpVxyshhFi1ahWJiYlcuXKFYcOGUa9ePQn2hZAE\nfCFEvi1fvpxy5cpRoUIFTp06xbx586xdJJEDSekIUYCKY0pHWJakdIQQQuSaBHwhhCgmJOALIUQx\nIUMrCFGAfH1976iTjxA3Y8nWTvLQVgghCsDba9+mpHNJRrccXSDXl4e2QghRCGitWRa4jG51ut3+\n4DzI68CcEvCFEMLC9oXuw97Gnrql6lr82leuwDPP5O1ciwR8pdRZpdRBpdR+pdSu9G1eSqn1SqlA\npdQ6pVThmyBTCCEKwNJjS+leu7vFn9+EhMBjj0EeR4G2WA3fjDFheUOtdZP0bSOAjVrrWsAmYKSF\n7iWEEIXa0mNLLZ7O2b8fHn7YqN1Pn563a1gq4KscrtUFmJn+fibQ1UL3EkKIQiswIpCoxCialG9y\n+4Pv0Lp10LYtTJoEI0fmaYIzwHIBXwMblFK7lVKvpm/z0VqHA2itw4DSFrqXEEIUWkuPLaVb7W6Y\nlGXC66xZ8OKLsGwZZBmWP08s1Q7/Ea11qFKqFLBeKRWI8SOQ1U3bXmadMMHPzw8/Pz8LFUsIIe6u\nJUeX8Hmbz/N9Ha1hwgT46ScICIDw8AD8/QPydU2Lt8NXSn0MxAOvYuT1w5VSZYDNWus6ORwv7fCF\nEPeE87Hnqf9TfcKGhWFnc2cTyOTEbIZhw+DPP2HtWihf/sYDlI3N3W+Hr5RyVkq5pr93AdoC/wIr\ngJfSD+sLLM/vvYQQojBbdmwZHWt2zFewT02Ffv3gn3/gr79yCPZaw8CBebq2JVI6PsBSpZROv95s\nrfV6pdQeYIFS6hUgCOhpgXsJIUShteToEgY3HZzn85OSoHdvSEyEDRvAxSXbAVrDkCFw4ECeri9D\nKwghhAVEJERQ7dtqhA4LxdnOOdfnX7kCXbuClxf8/jvY22c7QGsYPtzI8/z5J8rLS4ZWEEIIa1gZ\nuJLHqz6ep2AfGwtPPQUVKsDcuTkEe4CPPzbaZ65fD56eeSqjBHwhhLCAxUcX83Sdp3N9XmQktGkD\n9erBtGlgY5PDQWPHwuLFRp6nRIk8l1ECvhBC5FNMUgxbgrbQsWbHXJ136RK0amUs//sfmHKKyJ99\nZlT7N22C0vnrziQBXwgh8mnl8ZX4VfbD3eHOB7kJDzcCfefORnv7HHvPTphg9LzatAl8fPJdTgn4\nQgiRT4uOLOKZ++58CMvQUPDzg5494ZNPbhLsJ06EqVONYF+2rEXKKQFfCCHyIS45jk1nNtGpZqc7\nOv7CBSPYv/ACfPTRTQ6aOBGmTIHNm3NoiJ93MsWhEELkw+oTq3m00qN4OXnd9tgLF4w0Tr9+RgvL\nHGUN9hUqWLSsEvCFECIf7jSdExICrVvfJth/8QX88kuBBHuQlI4QQuTZlZQrbDi9gS61utzyuJAQ\no2b/yiu3CPaff27k7AMCCiTYg9TwhRAiz9aeXEvT8k0p4XzztvGhoUbN/uWXbxHsx441ml7+9ZfF\nHtDmRAK+EELk0e3SOeHhRqeqF16AESNyOEBr8PeHRYuMmr0Fml7eiqR0hBAiD66kXGHtybV0rZ3z\nZH6XLhnB/tln4YMPcjhAa2P6qqVLjZx9AQd7kBq+EELkyarjq3i4wsOUdrmx9+vly/D448ZgaDk2\nvdQa3nkHtmwxgn0+hkvIDQn4QgiRB3P/m0vv+3vfsD0qCp54whgMLcdOVWYzvPkm7NtnjHyZx4HQ\n8kKGRxZCiFyKTorG9xtfgocE4+Hokbk9NtYI9s2awddf5xDs09Lg1Vfh5ElYswbc3PJcBqWU9YZH\nVkqZlFL7lFIr0te9lFLrlVKBSql1SimP211DCCGKgqVHl9K6Suvrgn18PLRrB40b3yTYX70KffoY\nva/++CNfwT6vLPnQ9m3gSJb1EcBGrXUtYBMw0oL3EkIIq5l3eN516ZyEBOjUCWrXNka9vCHYJyXB\n008brytW5DCV1d1hkYCvlKoAtAemZtncBZiZ/n4mkPOjbCGEKEIuXrnIP+f/yRwKOSkJunUz+kpN\nmZLDEMfx8dCxIzg7G80vHR3vfqHTWaqG/zXwHpA1Ge+jtQ4H0FqHAfkbyFkIIQqBRUcW0aFmB5zt\nnElJgR49wMMDpk/PYfKSyEgjqV+lCsyeDXZ5n9zcEvId8JVSHYBwrfUB4FYPEOTJrBCi6Fm9GqKj\nM1fn/jeXFyp1Im3Fanr3NoL87Nlgm73NY1iYMSxm8+ZG1T/HqazuLks0y3wE6KyUag84AW5Kqd+A\nMKWUj9Y6XClVBrh4swv4+/tnvvfz88PPz88CxRJCCAt45BGj59RnnxGsYjkf/B+P/xPAgIjxJCTA\nsmU5VNzPnoW2bY0uth9+eJMB73MnICCAgICAfF3Dos0ylVItgWFa685KqS+Ay1rrCUqp4YCX1vqG\nzsXSLFMIUehFR8MHH/B9S2fqzlzDQs9tHAvzZNUqcHLKduyRI/Dkk/DeezB4cIEVKS/NMgsy4HsD\nC4CKQBDQU2sdncM5EvCFEIWePnMGVbUqA15azJGT3fnjjxwa2/zzD3TpApMmwXPPFWh5rB7w80IC\nvhCi0IuOJuztfrT12Mfo+e14cvfnuFfK1kN2wwYjyE+fDh06FHiRrNrxSggh7jmrV0NQEHrUB/Qs\n7UHExVdpu2E47gOfv+5BLvPmwfPPw+LFdyXY55UEfCGEuJlHHoGBA5lkP5jtdstZN7IzHj9/Ad9/\nD9u2Gcd89x28+y5s3AgtWli3vLchKR0hhLiFb/yj8Zz2LGsGxbPgXAP47DNjwDOtYfRoWLAA1q+H\nypXvarnyktKR0TKFEOImvv4avv/dkyZvpbBg+HY4M9sI9levwuuvw6FD8PffULpo9CuVlI4QQuTg\nhx/g229h4ZxAWi/dTkLgYZg40Rj8rGtXY+7CzZuLTLAHCfhCCHGDqVNh/HgIWBaNmvAK+958Guea\n98HQodCokVHLX74cXF2tXdRckZSOEEJkMXOmMc3s5s1QKfBvnm4aztePjjfGsH/qKaP3rJ+f1cfF\nyQsJ+EIIkW72bBg1ypiIqkYN2GzrQlKwI4+G2EL3FjBmDLz2mrWLmWcS8IUQApg/3xgNYeNGY1x7\ngJ/2/sSEK81RXbrAjBnQvr1Vy5hfEvCFEMXe4sUwZIjRuvK++4xt4fHhVJm5knZ7PY3pCB980LqF\ntABphy+EKNaWLIGBA41ZBxs0SN+YlsaeZ1tQdvdRym85AL6+Vi1jTmRoBSGEyIVly4xgv3ZtlmAf\nH4/u2pV4wQSrAAAgAElEQVSkf/dzcd3SQhns80oCvhCiWFq+HAYMMLI1DRumbzx/Hlq04LxjCu++\nXYeGtf2sWUSLk4AvhCh2li83GtusXm00qwdg3z5o1gx69+bN7g68+vBAq5axIEjAF0IUKxnB/rrn\nsMuWGZOWTJ7M6f7PsO3cdnrd38uq5SwI0kpHCFFsLFt2LY3TuDHGAGjjxxujX65ZAw89xDdrB9O/\nUX9c7YtWL9o7ke+Ar5RyALYA9unXW6S1HqOU8gLmA77AWYwZr2Lyez8hhMiLxYuvPaBt1AhISoL+\n/eHoUWOmqvLliUyM5LdDv3F44GFrF7dA5Dulo7VOBlpprRsCDYB2SqkmwAhgo9a6FrAJGJnfewkh\nRF7Mnw+DBhlNLxs1AsLCoHVrSE6GLVugfHkAftrzE11qdaGcWznrFriAWCSHr7VOSH/rgFHL10AX\nYGb69plAV0vcSwghcmP27Gudqho2BPbsgYceMsbFmTcPnJ0BSE5N5n+7/sewZsOsW+ACZJGAr5Qy\nKaX2A2HABq31bsBHax0OoLUOA4rOGKJCiHvCzJnXhkuoVw8jwLdrB5Mnw0cfgelaCJzz7xwe8HmA\nB3wesF6BC5hFHtpqrc1AQ6WUO7BUKVUXo5Z/3WGWuJcQQtyJn34yJqfatAlq10iDER8YuZ2NG6F+\n/euO1Vozacckvn7yayuV9u6waCsdrXWsUioAeAoIV0r5aK3DlVJlgIs3O8/f3z/zvZ+fH35+fpYs\nlhCimPnmG6MSHxAA1byjoENvSEmB3buhZMkbjl97ci02Jhser/r43S/sHQoICCAgICBf18j3WDpK\nqZLAVa11jFLKCVgHjAdaApFa6wlKqeGAl9Z6RA7ny1g6QgiLGTcOfv3VGOK4Utxh6NIFOnY0ZqvK\nYQx7rTVNpjbh/ebv06NuDyuUOG+sNadtWWCmUsqE8UxgvtZ6jVJqJ7BAKfUKEAT0tMC9hBAiR1rD\nBx8Ybe3/+gvKbVtotMP88kvo2/em5608vpKUtBSevu/pu1ha65DRMoUQRZ7ZDIMHw44dsG51KiW/\nGgULFxqN7zPHTsjhPG2m0c+N8Pfzp2vtotWQ0Fo1fCGEsJrUVHj1VTh1CjYvuIT7872N1jd79kCJ\nErc8d8nRJdiabOlSq8tdKq11yVg6QogiKykJevSA0FBYP3Yn7q0aQ5MmRnfa2wT7NHMaHwd8zNhW\nY1EqVxXlIktq+EKIIik2Frp2hVIlNQvbfY/ts2Nh6lTo3PmOzp9/eD7uDu60q96ugEtaeEjAF0IU\nOZcuGf2nHqkXx9eJAzBNPwzbt0P16nd0flJqEqM3j2ZKxynFpnYPktIRQhQxZ8/Co4/Ci43+45vt\nD2FydYGdO+842AN8teMr6vnUo03VNgVX0EJIavhCiCLj0CFo3x5+9ZtJ26Xv3rbJZU7Ox55n0o5J\n7O6/u4BKWXhJwBdCFAl//QV9n7nC1vsGUWXfLti8Ge6/P9fXeX/D+wx8cCBVvaoWQCkLN0npCCEK\nvUWL4KNu/3LY6UGqVMEYIiEPwX5L0Bb+Dv6bEY/e0Om/WJCAL4Qo1L7+SrO7/8/8qVvj8skImDED\nXFxyfZ1UcyqD1w7my7Zf4mKf+/PvBZLSEUIUSmYzfDgoCr/Z/RlU8SS2i7dC7dp5vt4X276glEsp\netxXdMbLsTQJ+EKIQichAT5tv43BO5/Ds29n7Cf/Do6Oeb7egbADfLPzG/a+trdYNcPMTgK+EKJQ\nCTufysqmn/J+5E+4zJ6C3dN31pHqZpJTk3lh6Qt82fZLKnpUtFApiyYJ+EKIQiPwjzPEdX2e1hVc\n8Di1H1WubL6v+dHmj6jhXYMX6r1ggRIWbfLQVghhfVpzcNgsSrRvgnq6O9WO/2GRYL81aCuzDs3i\n544/F+tUTgap4QshrEpfjuREm9dxOHyECzM20vjF+rc/6Q6cjz1Pr8W9mNZ5GqVcSlnkmkWd1PCF\nEFZzddU6onzr88+F8jj9t4f6Fgr2iVcT6TqvK281eYv2Ndpb5Jr3AktMcVgBmAX4AGbgF631t0op\nL2A+4AucBXpqrWNyOF8mQBGiuImPJ2HQe8TOX8MPjX/l/XVtcHW1zKW11jy/9Hm01szuPvueTeXk\nZQIUS9TwU4F3tNZ1gWbAIKVUbWAEsFFrXQvYBIy0wL2EEEXd33+TVKcBqxclMmPoIfy3Wi7Yg9He\nPjAikGmdp92zwT6v8h3wtdZhWusD6e/jgaNABaALMDP9sJlA0Zo/TAhhWYmJMGwYCZ168lrMl9jN\nnsGIcR6YLJhY/nnPz/yw5weW9VqGk52T5S58j7DoQ1ulVGWgAbAT8NFah4Pxo6CUKm3JewkhipAd\nOzC/9DJ7zQ150/sQv64oSd26lr3FlL1T+Pzvz9n04iYquFew7MXvERYL+EopV2AR8LbWOl4plT0x\nf9NEvb+/f+Z7Pz8//Pz8LFUsIYQ1JSTA6NGk/T6HD92/49j9z7B+Bnh4WPY2v+z9hU+3fMrmvpup\n5l3NshcvJAICAggICMjXNfL90BZAKWULrALWaq0np287CvhprcOVUmWAzVrrOjmcKw9thbgX/fUX\nvPoqoRWb8PjhybwwtCTDh4Ml0+paa77Y9gXf7/6eTX03Ud37zidBKery8tDWUjX8X4EjGcE+3Qrg\nJWAC0BdYbqF7CSEKs5gYGDECvXIlcx79kfe3duL3edCqlWVvk3A1gX4r+nEy8iTb+22XNM4dyPfj\nEqXUI8BzQGul1H6l1D6l1FMYgf4JpVQg0AYYn997CSEKuWXLoG5dEhI0Xav9x9TwTuzda/lgfzb6\nLI/8+gh2Jju2vLRFgv0dskhKJ18FkJSOEEXfhQsweDD8+y+7X/uFLl+1pF8/8PcHGxvL3SbNnMZP\ne37C/y9/Rj06iiEPDym2TS+tmdIRQhRHaWnw44/g70/aa28wutLvzPrGidmzLV+r/zf8X15b9Rq2\nJlu2vLSFOqVueCQobkMCvhAib/bvh9dfBwcHzv62lZ4f16F0aThwAEqWtNxtjl46yri/x7H25Fo+\nbfUp/Rv3x6RkVJi8kG9NCJE7MTFG+uapp9CvDWBKnwAeerEOL70EK1daJtibtZktQVt4ZsEztJzR\nklolanHirRMMeHCABPt8kBq+EOL2Vq+G5s1h7Vp47z3o0IFLK3bw4+CjLLtqYssWqJPPDEuaOY19\noftYeGQhc/+bi7eTNy83eJmZXWcW2zloLU0CvhDi9jw94f77jer7okUsP16HyFYfoAd8xs4JYG+f\nu8uZtZlzMec4GnGUA2EH2Bq8lW3B2yjvXp6utbqy9rm13F/6/oL5LMWYtNIRQtxcVJTR1GbuXBg+\nnIR/j/Nu9Js02/o13jOGUqO5A8mpyaSkpZCclkxyajJJqUkkpSZx5eoVYpNjiUmKISY5hrD4MELj\nQwmJC+FM1Bk8HD2oU7IO95e+nxaVWvCY72Mybn0u5KWVjgR8IYo5szYTHBPMqchTBMUEERwTzIWo\nYBqu3kuvBUdZ/4ATY9rac0rFU+5yMme/1TQa4U1EaRccbB1wsHHAwdYBext7HG0dMxdnO2c8HDxw\nd3DHw8EDH1cfyrmVo6xrWSp7VsbD0cLjKxQz0ixTCHFLKWkpHAw7yD8X/mF3yG4OXzzMsYhjeDh6\nUMO7Br6evjx22sybP/8Fbm4Ez59C+crNqOrvgcshMyvrfApn3mffxIkw/DMj1SOKDKnhC3EPSzWn\nsjdkLxtPb2TD6Q3sDtlNNa9qNC3flCblm/CAzwPUKVnHqG2fPGk8kD1wACZMwPx0D6ZOU3z4IQx+\nMZoR8R9gOz49yEdHwwcfwGcS9K1FUjpCCJJSk9hwagNLji1hZeBKyrqV5YmqT/B41cdpUakFbg5u\n158QGQmffgqzZsG778KQIfx30pEBA8BshilT4IHg1fDII9cH9+ho2LYNOnS4ux9QABLwhSi2Mtqt\nzzgwg2XHltGgTAO61+lO19pdqeRRKeeTkpLgu+/giy/gmWfA3584Zx8++QSmT4dPPoHXXsOiE5QI\ny5EcvhDFTHh8OL/s+4Vp+6fhZu/GSw1eYvzj4ynjWubmJ6WlwZw5MHo01K8PW7eia9Vm/nyjgt+6\nNfz7L5S5xSVE0SQBX4giaNeFXXy36ztWHV9Fj/t6sLjnYhqWaXjrgcS0hjVrYORIcHWF336DFi04\ncACGtjYyO/PmwaOP3r3PIe4uSekIUURordlwegPj/h7H6ajTvNXkLV5p+AreTt63P3nLFuMha2Qk\nfP45dO7MxUvGA9nly42m9v37g61UAYsMSekIcQ/SWrPy+ErG/jWWxNRERjwygl7398LOxu72J+/e\nDR9+CCdOwJgx0KcPCck2TB4PkybBiy9CYKA0tCkuJOALUUhl1Og/3PQhSalJjPEbQ5faXe5s8LB9\n+4xq+759RsB/5RXSbOyZNQs++ggefhh27IAaNQr8Y4hCxCIBXyk1DegIhGut66Vv8wLmA77AWaCn\n1jrGEvcT4l63J2QP7214j9C4UMb4jaFH3R53Fuj374exY+Gff2DECFiwAO3gyNKlxjNab29YsACa\nNSv4zyAKH0tNYv4oEA/MyhLwJwCXtdZfKKWGA15a6xE5nCs5fCHSBccEM+rPUWw6s4kxfmN4ueHL\n2JruoF62e7fRjnLvXqOpzeuvox2dWL/eqOCnphp9pNq1s+wk4sJ68pLDt0gLW63130BUts1dgJnp\n72cCXS1xLyHuRQlXE/h488c0/LkhVb2qEvhmIP0b9791sNfaeBj71FPQvTs8+SScOoUeMpQ//nKi\nWTMYOhTef9/4HWjfXoJ9cVeQOfzSWutwAK11mFKqdAHeS4giSWvN4qOLGbZ+GM0qNOPAgANU9Kh4\nu5OM5pXjxkF4OAwfDsuXY7ZzYNUqoyYfH2/k6p95xrJzyoqi7W4+tL1p3sbf3z/zvZ+fH35+fneh\nOEJY14nLJxi0ZhCh8aHM7DoTv8p+tz4hJcVoKD9xohHFR46EZ54hVdswbx6MH2+MSz9yJDz9tPSQ\nvdcEBAQQEBCQr2tYrB2+UsoXWJklh38U8NNahyulygCbtdY3zIkjOXxR3CSlJjHh7wl8t+s7RrUY\nxeCmg2+duomOhqlTYfJkqFXLyNE88QSxcSpzc5UqRqBv21bSNsWF1XL4GfdPXzKsAF5Kf98XWH7T\nM+fMMWovQtzjNp/ZTL0f63Ho4iH2D9jPO83euXmwP3MGhgyBqlWNESyXL4eNGzlToy3vvqeoUgV2\n7YJFiyAgwEjhS7AXt2KRgK+UmgNsB2oqpYKVUi8D44EnlFKBQJv09ZxNm2ZUUT77DC5dskSRhChU\nIhMj6be8H32X9eXLtl+yuOfinHP1WsOmTdC1Kzz0EDg4wKFD6N9+Z8PlRnTubGxWymhiP2+esS7E\nnSg8QyscOmT8bbpkCXTuDIMGQZMmVi2bEPmltWbhkYUM+WMIz9z3DJ+2/hR3B/cbD4yPh99/h++/\nN8YkHjwYnn+eiEQXZswwhih2cDA2P/ccODvf9Y8iCpl7Y3jky5fh11/hhx+MCZPfeAOefRZcZNZ6\nUbSExIUwcPVATkSeYGqnqTSrmENvpyNH4McfYfZs8PODQYMw+7Xmz02KX3+FtWuN+s/rrxudpSRl\nIzLcGwE/Q1oarFtn/M+wfbtRrenfHx544O4XUohc0Fozbf80Rv45koEPDmRUi1E42DpcOyAx0Ui8\nT5kCp05Bv37w2mucSKrIb7/BzJlQogS8/LLxz977DsZGE8XPvRXwswoKMlopTJ8OFSrAq68atX43\nt1ufJ0R+rc7dTE9nos7Qf2V/7t8dxKtvTef+mlnGGt6yxZhsZMcOI/E+YAAXm3RkwVI7fv8dzp6F\nXr3gpZegQYMC/2SiiMtLwEdrbdXFKMIdunpV61WrtO7SRWsPD61ffFHrzZu1Tku782sIkRtRUVoP\nHGi85rSeLs2cpr/d+a0uMaGEnvD3BH318iXjuJMntf72W63r1dPa1VXr4cP15b1n9JQpWrdpY/wz\nfu45rdeuNf55C3Gn0mNnruJt0ajh5+TiRWN87+3bISEBnn/eWHx8ZJ5NYVkZE3a/957R6SnbxN2B\nEYH0W9EPgGmdp1HLrbLRE3baNNiwAdq1I/GKmd/bTGfu+hLs3Ws0oXz2WWO4AycnK30uUaTduymd\nm4mOhlGjjL+Dly83HnxpbbRd7tsXypWzbGFF8XX2rNF0+MwZqFwZgFRzKpO2T2Li9omMafERbyTd\nj2nOXFiyBF2vHmdbvMCaiw8x6Od6NPA4Q/0ulenWzegcJa1sRH4VvwlQPD2N2Xsyal8JCfDEE7By\nJdx/v5EIffZZY2CpUqWsXVpRVEVHGzX7M2cya/gHks7y6rJXeCTEhtMXO+H+zXhSvUtzsN5zzGx1\nkHl/V6Di+Wgmu37Azrln2PvXRGzGfSYzjQirKto1/Aw51L5ISjLatC1YYLw2bgw9ekCXLlC2bH6L\nLYqLjHROehonMew8p55uw560YJ4544R2KMP2ir34MfJZNofUws/PGLyyXbNoKv9y7bzs1xEiv4pf\nSgdum18FjGZwf/wBCxcawf+++6BbNyP4y5Q/4lZWrzYqC3v3cvbXH3Fe/wfRTi4c5wn848fi8tB9\ntGkDbdoY/QTt7LKcl4vWPULkVvEL+NlrTXdSi0pOhs2bjR69q1aBuzt06mQszZpl+T8W+Z+2GLt6\nLoyQqWswr1iFz+GN7PdyZkmdZAJdxtOk6QBatjQCvDxwFdZS/AJ+fgOy2WzMDLFypRH8z5wxngG0\na2c8WXNxyf0PiihyUlLgyKFULizeic3GP/A9+gdlrpxip3tb9jbx4NfGy3i4bm9+ePpTPJ2l74co\nHIpfwLe00FCjd++aNfDnn0au/7HHjI5f48fDTz9JsC/CzGbjcc/hw3D4P03kjkA892zkgfCNtCSA\nSI+qhDd6CofOT5HU1o33d7xNwtUEfuzwIw+VlxHKROEiAd+S0tKM4QjXrzd+ALZvh4YNjSdyjz1m\n/GUhPX0LHa0hMhJOnoQTJ4wlMBCOHdXo4yd40vEv2jluplFcADYOtiS1eAL3bo9j/1Rr8PEhJimG\nT7Z8wsyDMxnrN5bXGr+GjUmmjBKFjwT8gpCRxhk82JhKrkYNY8LoPXugTh0j8DdvbiwVKli7tPc8\nrSEiAs6dg+Bg44+v4GAjG5exANxXLZnWJQ7QnB3cH72Vsqf+xsbFEdNjLaBVK2OpUiVzNDKzNjPj\nwAw+2PQB7au35/M2n+Pj6mPFTyrErUnAt7RbPRR2cjIC/7ZtRu1/+3ZwdDSe5DVtarw2bAgeHtb+\nFEVCYqIxFcKlS8Y0rRlLaOi1JSQELlwwHq2ULw++vulLJc39jiepFbebcud343DwH9TBg8aPc7Nm\n0KIFPPooVKqU4723BG1h2Pph2Jps+fapbyV9I4oECfiWlpuHwlrD6dPwzz/GNES7dhlj/JctazTr\na9AA6teHevWMHsD30Di3WhvdHuLjIS7u2hIbCzEx15aoKGOJjjZGwY6MNJaICEhNNfrGlSxpjI6R\nsZQte20pVw7Kl0zGOeio8d3u328sBw4YP6wPPXRtadIEXF1vWe5jEccYvnE4B8MO8nmbz+l1fy9M\nSiaCFUVDoQz4SqmngG8wZteaprWekG2/3rYt92W4WbzUWpOik0hOSyTZnMhVncRVczLJacbrVZ3M\nVXMKqTqFNJ1Kmk7FTBpmnQZodPpc6yZlylzsTHbYmeyxM9lja7LHydYZZ1sXHG2ccbVzw8XOHXsb\nuxvKpdJScTgbiOPhvTgGHsQx8CAOgQfBbCalel2Sq9clqVpdUqrWJrlyLa76VACl0JpbLmbztdfs\nS1pazktqqrFcvXr9+5QUY7l61WixmrGenGwE8axLYuK1JSEBrly5ttjZGTVvd3fj0YabmxGD3d2v\nvXp5GYunpzH8b4kSxtC/JUsasfm6/6YJCUYi/uhRYzlyxHjaevq0MeXfAw8Yf0FlLLnoSR0cE8wn\nf33CssBlDH9kOG82eRNHW8fc/QMUwsoKXcBXSpmA4xhTHIYAu4FeWutjWY7RtTutxWxKRJuSMNsk\nkmaTgNmUiNkmEbMpwXi1ScBsSjD22SRgtrlCminjfZbFlIQy22MyO2FKc0SZHTGlOWEyO6LMDpjM\nDiizPcpsh9J2oG1QZlvQJjKn5dWA0mjMaJWGVqloUwralILZlITOKI/NFcx2caTZxqK0HTYpntik\nlMAmuYTxmuSDbVIZbBLLYJtYDruECtgmlMcnIY0ayUeonnyY6imHqZISSJXkQNzM0QTbVyfIvibn\n7Ktx3qEaUbYl8U69xKoSfUmxccJkMgKjUmBjQ+Z6xnuTyXifdbG1NRYbGyMw29ld2+bgcG2bg4Ox\n2Nsbr46O116dnIzF0dEI7M7OxrqLi7HY5naQDrPZyN8EBRlNZ06fNpZTp+D4caPaX6WK8ZzkvvuM\n17p1oXZto1B5EBIXwrit45jz3xxeb/w6w5oPw9tJBpsXRVNhDPgPAx9rrdulr4/AGNJzQpZjdNvf\n2uJk64SjrSNOdk442TrhbOec+eps54yTnRMudi6Z6y7219472TrhYu+Ck60TTnZOd/3Pcq01iamJ\nRCVGEZkYyeXEy0QkRBAeH074lXDC4sMIiQvhfOx5zsWeI/FqIr6evlT2rExlj8pU965OjRI1qJPg\nQtXPf8CmfUcICzNqtps3GwHu/HmjGlyhgpHALl/eyHmUKWO8ZuRDMqrN9vZ39TtI/yKMPE5Gvubi\nxWtLWJiRhA8NNRLx588b1f6KFY3AXqWKUXOvWhVq1jTy7TaWaR1zOuo0X2z7ggWHF/BSg5cY8egI\nSruUtsi1hbCWwhjwnwae1Fq/lr7+PNBEaz04yzGFN4dfQOJT4gmKDuJM9BnORJ3hZORJTkSe4ETk\nCWLDgvlqqxMBPR7ilQ2XCX5/ALWqN6W2d00co+OvBcsLF6492QwLM2rEERFGrTk62qiye3kZ+RRX\n12tL1qq6vf31Vf6MnIpS1+eHUlONHE/GkjWfExdnJOhjY43Fycn4wfH2htKlry1lyhhJ+IylYsUC\n76a6J2QPX+34ivWn1vP6g6/zdtO3KeUig+iJe0ORHS3T398/872fnx9+fn5WK8vd4GrvSt3Sdalb\nuu4N+5JTkwk+uIXnHmzLt3OHsj1sM/8d+o5TUaeo4lmF+mXqU9+nPg3qNKBhmW45Nx3U2gjK0dHG\nkhGY4+KuT8onJxvJ+4zEfsa5Whs/AllzQhn5HgcHI5+TkcvJmqx3d7fOXxZZHq5fTbvKkqNLmLZ5\nEuUPnaXhS+/xY4cf8XCU1lKiaAsICCAgICBf17gbKR1/rfVT6es5pnSKWw3/lm4yGFxyajLHIo5x\nMPwgB8MOciD8APtD9+No60jDsg1pVKYRjcs1plHZRlR0r4i6h1oB3VZ0NJHvDOTrjiX45fRCHnSu\nzuS/3fD932/Yepe0dumEKBCFMaVjAwRiPLQNBXYBvbXWR7McIwE/Qy4Hg9NaExwTzL7QfcYSto+9\nIXtJ02k0LtvYWNJ/BHw9fO+5H4FzMedYdGQR8w7PIyYsiOm7yuMzZiJVpy6WITDEPa/QBXzIbJY5\nmWvNMsdn2y8BP4MFRufUWhMSF8Le0L3sDdnL3tC97A/bT1JqEo3KNqJhmYY0LNOQBmUaUHPnCWxa\nPFZkRgPVWvPvxX9Zc2INK4+v5FjEMbrW6kqPuj14vOrj2Aafv3FeBCHuUYUy4N+2ABLw74qw+DD2\nhe5jf+h+9ocZS+KlUCb/7caW/m2pUe0h6jv48uD3S3H+4muUl5e1i4zWmhORJ9gatJWtwVvZeHoj\nDrYOtK/eng41O9C6SmvsbdKfGdzJvAhC3EMk4ItciU2O5cjx7biOHcfcpyrw4JzNDHssiWhHqF2y\nNjVL1KRmiZpU965uNCH1rEwp51IFkhqKT4nnVOQpjkUcY3/Yfg6EHWB/2H4cbBxo4duCRys+Susq\nralZouaN98/LvAhCFHES8EXeZJkiUvv6cvHKRY5fPp65nIw6SVB0EGejz5JwNYFybuUo41qGsm5l\nKeVcCi9HL7ycvPBw8MDJzuhP4WDjgFIKszZj1maSU5OJS4kjPiWemKSYzP4JofGhnI0+S0xSDFW9\nqlKzRM3MlFODMg2o6FHx9uWXiWpEMSQBX+ReLlMh8SnxhMaFEhofSmhcKJcSLhGVGEVUUhQxSTEk\npyWTlJpEUmoScG2ICnsbe9zs3XBzcMPN3o0yrmUo41oGH1cfKntWppxbORnHRohckIAvckdSIUIU\nWRLwRe5IKkSIIksCvhBCFBN5CfiSNBVCiGJCAr4QQhQTEvCFEKKYkIAvhBDFhAR8IYQoJiTgCyFE\nMSEBXwghigkJ+EIIUUxIwBdCiGIiXwFfKfWMUuo/pVSaUqpRtn0jlVInlFJHlVJt81dMIYQQ+ZXf\nGv6/QDfgr6wblVJ1gJ5AHaAd8IO61+bXKwD5naD4XiLfxTXyXVwj30X+5Cvga60DtdYngOzBvAsw\nT2udqrU+C5wAmuTnXsWB/GO+Rr6La+S7uEa+i/wpqBx+eeBclvUL6duEEEJYie3tDlBKbQB8sm4C\nNPCB1nplQRVMCCGEZVlkeGSl1GZgmNZ6X/r6CEBrrSekr/8BfKy1/ieHc2VsZCGEyIPcDo982xp+\nLmS98QpgtlLqa4xUTnVgV04n5bbAQggh8ia/zTK7KqXOAQ8Dq5RSawG01keABcARYA0wUGY5EUII\n67L6jFdCCCHuDqv2tFVKPaWUOqaUOq6UGm7NsliTUqqCUmqTUuqwUupfpdRga5fJmpRSJqXUPqXU\nCmuXxdqUUh5KqYXpHRgPK6WaWrtM1qKUGpre0fOQUmq2Usre2mW6W5RS05RS4UqpQ1m2eSml1iul\nApVS65RSHre7jtUCvlLKBPwPeBKoC/RWStW2VnmsLBV4R2tdF2gGDCrG3wXA2xjpQAGTgTVa6zpA\nfTLSTEgAAAKJSURBVOColctjFUqpcsBbQCOtdT2M54+9rFuqu2o6RqzMagSwUWtdC9gEjLzdRaxZ\nw28CnNBaB2mtrwLzMDpsFTta6zCt9YH09/EY/1MXy34LSqkKQHtgqrXLYm1KKXeghdZ6OkB6R8ZY\nKxfLmmwAF6WULeAMhFi5PHeN1vpvICrb5i7AzPT3M4Gut7uONQN+9s5Z5ymmQS4rpVRloAFwQxPW\nYuJr4D2Mvh7FXRUgQik1PT3FNUUp5WTtQlmD1joEmPT/du6eNYooCuP4/wEFixSCEEUkARE/QCox\njRI7wVoUfKtF+zS2thY2ghYKWrgIWomIH0CFFIJlwESLtRArGwmPxb2FKLgI4lm4z6/ZmYEZDsxy\n5hwu5wJbtEHOr7Zf1kZVbtH2FFrRCCzOuiG7Zc4RSQvABLjeK/2hSDoNTHu3I37fsmM0u4AV4Lbt\nFeAbrY0fjqS9tIp2GTgILEg6VxvV3JlZJFUm/E/A0k/nh/q1IfU2dQI8sP20Op4iq8AZSZvAI+Ck\npPvFMVX6CGzbftvPJ7QPwIhOAZu2v9jeAZ4Ax4tjqjaVtB9A0gHg86wbKhP+G+CIpOW+2n6WNrA1\nqnvAe9u3qgOpYnvd9pLtw7T/wyvbF6rjqtLb9W1JR/ulNcZdzN4Cjkna03feXWO8Bexfu95nwKV+\nfBGYWSj+y0nbv2J7R9JV4AXtw3PX9mgvEABJq8B54J2kDVprtm77eW1kMQeu0abWdwObwOXieErY\nfi1pAmwA3/vvndqo/h9JD4ETwD5JW8AN4CbwWNIV4ANtS/o/PyeDVxERY8iibUTEIJLwIyIGkYQf\nETGIJPyIiEEk4UdEDCIJPyJiEEn4ERGDSMKPiBjED2/eY7K5BfIBAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7fd3468c1ed0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"from sklearn import linear_model\n",
"import matplotlib.pyplot as plt\n",
"\n",
"order = 9\n",
"\n",
"# Needed to parse for matplotlib to parse latex input\n",
"# plt.clf()\n",
"# plt.rc('text', usetex=True)\n",
"# plt.rc('font', family='serif')\n",
"\n",
"for alpha in [0.15,0.3,0.5,1,2,4,8]:\n",
" x,y=gen_poly_data(10,True,3)\n",
" x_train = gen_poly_features(x,order)\n",
" y_train = y\n",
"\n",
" xt,yt=gen_poly_data(100,False)\n",
" x_test = gen_poly_features(xt,order)\n",
" y_test = yt\n",
"\n",
" clf = linear_model.Lasso(alpha=alpha)\n",
" clf.fit(x_train.T,y_train)\n",
"\n",
"\n",
" ws=compute_ridge_weights(x_train,y_train,alpha)\n",
" mse_test,yer=compute_mse(ws,x_test,y_test)\n",
" # clfr = linear_model.Ridge(alpha=alpha)\n",
" # clfr.fit(x_train.T,y_train)\n",
"\n",
" ye = clf.predict(x_test.T)\n",
" # yer = clf.predict(x_test.T)\n",
"\n",
" plt.clf()\n",
" line1,=plt.plot(xt,ye,label=\"Lasso\")\n",
" line3,=plt.plot(xt,yer,label=\"Ridge regression\")\n",
" plt.plot(x,y,'x')\n",
" line2,=plt.plot(xt,yt,'r',label=\"Original curve\")\n",
" plt.legend(handles=[line1,line3,line2])\n",
"# plt.title(r'\\lambda = '+str(alpha)) # matplotlib appears to have bug to parse latex input\n",
" plt.title('lambda = '+str(alpha))\n",
" \n",
"# plt.savefig('regularized_fit_'+str(alpha)+'.pdf') # save figure\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 2 conda2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.