Bagging is a procedure that let's use bootstrapping to reduce the variance, and ultimately performance, of a prediction. The procedure can be defined as follows:
\n",
"We have data $D=[(X_1, Y_1),(X_2, Y_2),...,(X_N, Y_N)]$ and we want to learn:$\\:\\: E[Y|X]=\\hat{f}(X)$. Define a bootstrap sample $D^b$ as $N$ samples from $D$, sampled with replacement. Let $E^b[Y|X]=\\hat{f}^b(X)$ be the function learned from training set $D^b$.\n",
"Our bagged prediction is then the mean of all estimates of $f^b(X)$. I.e.,\n",
"
\n",
" \n",
"This is a relatively straightforward procedure to implement, which we can show on a simulated example.\n",
"\n",
"\n",
"\n",
"\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqKElEQVR4nO2df5ReVXnvP8+8MxPirRYNKBAIoRStdKWoTKPT1jCKUODWS1p6vV6iw0U0HSQs+Yeh2GXLWtwmiLY3VKDMCAHm1iX1ClXsStWCjNLOS2Uiv0NbIsUQAiUEf7QImczMc//Y75s5c+b9ffb59Z7ns9ZZ77zv2XP23ufs8z3PfvbezxFVxTAMw+h+etIugGEYhpEMJviGYRgFwQTfMAyjIJjgG4ZhFAQTfMMwjILQm3YBGnHEEUfo6tWr0y6GYRhGbtixY8dLqnpkrX2ZFvzVq1czPT2ddjEMwzByg4j8qN4+c+kYhmEUBBN8wzCMgmCCbxiGURC8CL6IbBORF0Xk8Tr7RUT+QkR2icijIvIuH/kahmEYrePLwr8NOKvB/rOBkyrbRuAvPeVrGIZhtIgXwVfV7wEvN0hyLjChjgeAw0XkaB95G4ZhGK2RlA9/JfBs4Pueym9LEJGNIjItItP79u1LpHBG91Euw5Yt7tMwDEdS8/Clxm814zKr6jgwDjAwMGCxm422KZfh9NNhZgb6++Hee2FwMO1SGUb6JGXh7wGOC3w/FtibUN5GwZicdGI/N+c+JyfTLpFhZIOkBP9uYLgyW+c9wE9V9fmE8jYKxtCQs+xLJfc5NJR2iQwjG3hx6YjIl4Eh4AgR2QP8CdAHoKo3AduBc4BdwM+BC33kaxi1GBx0bpzJSSf25s4xDIdk+RWHAwMDarF0DMMwWkdEdqjqQK19ttLWMAyjIJjgG4ZhFAQTfMMwvGPrILJJpuPhG4aRP2wdRHYxC98wDK/YOojsYoJvGIZXbB1EdjGXjmEYXrF1ENnFBN8wEqZc9i+GcRwzCoOD2SiHsRgTfMNIkDgGNG2Q1GgV8+EbRoLEMaBpg6RGq5jgG0aCxDGgaYOkRquYS8cwEiSOAU0bJDVaxYKnGYYRO1kbVO5mGgVPMwvfMIxYsUHl7GA+fMNIkCLGmLFB5exgFr5hJERRLd3qoHK13jaonB4m+IaRELUs3SIIvg0qZwcTfMNIiCJburbyNhuY4BtGQpila6SNr5eYnwVcB5SAm1X1mtD+XwT+ClhVyfPzqnqrj7wNI0+YpVsssjYdNbLgi0gJuAE4A9gDPCgid6vqzkCyS4CdqvpBETkS+BcR+ZKqzkTN3zAMI4tkcZDex7TMtcAuVX26IuB3AOeG0ijwehER4BeAl4FZD3kbhmFkkixOR/Uh+CuBZwPf91R+C3I98HZgL/AY8ClVna91MBHZKCLTIjK9b98+D8UzjPxSxHn73UIWYxz58OFLjd/C8Rp+G3gYeD9wIvD3InK/qv5syT+qjgPj4EIreCifYeSSLLoEjNbJ4iC9D8HfAxwX+H4szpIPciFwjbrAPbtE5N+AXwG+7yF/w2gLXwNpcQ/IFXXefjeRtUF6H4L/IHCSiJwAPAd8GDg/lGY3cDpwv4i8BXgb8LSHvA2jLXxZzUlY30Wet2/EQ2QfvqrOApuAbwFPAl9R1SdEZERERirJrgZ+Q0QeA+4FrlDVl6LmbRjt0ulAWtiXnsSAXNUlcPXV5s7pVpIeo/EyD19VtwPbQ7/dFPh7L3Cmj7wMIwqdWM21rPmkrO+suQQMf6QxRmPRMgtIkWd+dGI11/Olm/VtRCGNaZsWWqFg2MyP9q3metZ8nqzvrK34NNIZozHBLxg286N9sji9rh3sIZ9N0mhXJvgFw2Z+dEaWrPl2rfWJCXjtNVC1h3zWSLpdmeAXkAsucJ/Dw3bjt0varpF2rfVyGW691Yk9uFWfPh7yaZ8HozNM8AtEWCyGh9MuUb5IwjXSTEjbdclNTsJsJWqVCHzsY9HLbC6i/GKzdApEFoM55Ym4z19VSD/zGfdZaxZVu/FZgukPO8zPQz7r7SiOWWjdMrPNLPwCYf77aMRx/oIW/eQkHDgA8/Pus5b13u5AXxwDg1luR3H0PrqpR2OCXyDyPtskbXyfv7CQXHqpE3twnytW1C9HO3n7HhjMcjuKYxZaN81sM8EvGFmabZJHfJ6/sJA8/DD09Dix7+mB/fv95BMHWW1HcfQ+styjaRcTfMMIEMfsk3rHDAvJeefB/ffnX1jSnMETR+8jyz2adhHV7IacHxgY0Onp6bSLYRSERr7aTkWsmf83fNy8T3es1vfAAddLueEG2Lix/WPk+RykjYjsUNWBWvvMwjeMCvV8tVEG7Zr5f8OukbhcJUmJaHDgeX4eNm2CNWtq51mrTN00QJpFTPANo0I9X22UQbss+H+TFNGhoYVxCHDnrNb5qlembhogbZckHsom+IZRoZ6vNopoZ8H/m6SIDg46N86mTS6/Zctqn696ZcrCAzINknoom+AXgKL4RH3Us5ZLJapopz2jJWkR3bjRuXEana9GEUjTfkCmQVIPZRu07XKK4hMtSj07Jc6HfpQB7aIJez18tl8btC0wPiyHPNyYFhGyMXEOBncqVGn3fLJEUj0bE/wuJ2p3Pm3LuZWHTbkM27YtRITs7S2O7zdtijzI6pskHoBeBF9EzgKuA0rAzap6TY00Q8BWoA94SVVP85F3XknKao5qOaR5Q7f6sJmcdOUDFxHywgvjK2O9qYRZ7wHFRSsGRZHPT9aILPgiUgJuAM4A9gAPisjdqrozkOZw4EbgLFXdLSJvjppvnknaao5iOaQ5a6LVh024jHGFfa513WDht95e97Bp9p6BbhLAZgZF2j3ErJJWG/Bh4a8Fdqnq0wAicgdwLrAzkOZ84C5V3Q2gqi96yDe35KkbnOasiVYfNkmVsV5Y4Opvc3MwNga3315f2LpRABsZFHlq65CMEKfZBnwI/krg2cD3PcC7Q2neCvSJyCTweuA6VZ2odTAR2QhsBFi1apWH4mWPpKxmX403rcG1doQ8iTLWu279/QsDxs0GjfMmgLVop13laV59UkKcZhvwIfhS47fwXM9e4FTgdGA5UBaRB1T1X5f8o+o4MA5uWqaH8mWOJCzSbrEkszSTo951u/deN0to2zZ3EzcStrz7vNttV3maV5+UEKf5EPQh+HuA4wLfjwX21kjzkqq+ArwiIt8DTgGWCH5RiFvIusGSzCL1FmYNDjrffTNhy7vPu5N2laWHdiOSEuI0H4I+BP9B4CQROQF4Dvgwzmcf5OvA9SLSC/TjXD7/x0PeRh3y1JXuFloVtqR83nH0FOJuV90WWrlRXmk8BCMLvqrOisgm4Fu4aZnbVPUJERmp7L9JVZ8UkW8CjwLzuKmbj0fN26hPJ403y66EolBLUDu5LnH1FOIUxSz0bjoV4rzcO17m4avqdmB76LebQt8/B3zOR35Ga7TTeLNwsxlLBRU6uy6+XXphQeu2wcwo5One6Um7AEY2qDflMC+Uy7Bli/vMexkGB+HKK91np9el2lMolaK7XqqC9pnPuM+4znG4zCtW+L+mcbSTTq5RWu3VQisYQL59/lmwsOIqQ6fXxafrJSnLO1jmFSvgssv8ns+sXKM026tZ+AawcLNdfXW2u6S1yELvJK4yRLkuwZ5CFHz2FppRLfP+/f7PZ1auUZrt1Sx84xB5mT4XJu3eSbkMu3c7QQT/ZUj7uqQxjXDFCvfmLFV/5zPOdtLONcr7PHzDSBUfguTjJeW9vfCJTzSPpdNWYR57DLZudVHhPvUp93aRctmt9AJ45zudObxihfusqsfEBLzwwsIxn3kGXnoJzj8fPvvZ2vkFjxGqwCBlBpkEhmD8MbjlFjjsMDj5ZFeGhx5ayO/ll90T8HWvc2VeswauvRb27nXHPvzwhbrdcgsccwycffahvMsMctllzgLu6XHVr7Veod3rldSCx6hrMWJFVTO7nXrqqWoYcTM1pbp8uWqp5D6nplr/382b3f+B+9y8OXTgkRHVdetU165VHRtzv73jHapveIPqhg3u++bNC5mOjqr29KiKqPb26jwc2hTc/v7+ahQHt4m4z54et6+vb/H+8DY6urTyPT0LxwifhOAJanbsWlv12MHy9vbWTrdsmT5z/Dp9nJP1PtbpXazXf377enceK2WamlI9rX9KPy2b9bT+qYWiVs/l6KjqmWe68x2+0CMjqusXH88XUdqRT4BpraOpZuFXyMs8WsM/HQ9Kjo8z8rU7WaVHciJP8aa5l1l97U/gmyfDe94Df/ZnC3GbAb7/fWepVwP3f+lL8OUvu9/6++HSS50lXEFnZ4GF2CUKyF13wcGDi8tRPd78PBw8iKrWjHdyiLvuWrDyq5WvvnV8fn7pSQieoGB9WqV67GB5K3Vbku7AAVb96HuLf3+yst16K9x3H09NwPaZ0+lnhpmZfr46cS+D4Lpa1aBGAN/+tvus9oqGhlw9qmzbBl/4wtKeESyEXA3Oj52YgJ07XR4XXbTwHsef/AQmJ3nTfxzDn7z6Vk5jkudfPYanJkZduYK9pxUrXG+oepxgj6eaxwsvwFFHeegq1qDekyALW1IWflaezEY6hK//o2Mhq1vVWY2//MvOcjzzTGedVyzTJVZ4J1up5I4f+G0OWXTs+XoWfsA6n+3r1wP0NS5TChb+ovLUs/BFGp9DEdXNm/WZkc16ENetmqGkz4xsXtzVCm5nnunqsHnzQk8oeLy+Pvd/4Z5RX5/7rda+YJrAMcPXaq5UUl22bOHchvMPlqNWHsuWdSRGmIXfmLwu+DAiUvGFD77wAnvf9jIHn9vHf1l2kNeN/NDtP+ww52z92tcWLO9du5YcpqE1vSRxwMIH56SuWvi/93uLLPy9Gy5n81+fyCWzWxERll3+KU787EZYv76uD/9Lu4cYG4OP6ARv4QVE4L3vhSN+9kxtH354LmQtH37Y6fxYez78x1jDDy++lqPn9/K90hAf+sThHD88tNSH/9BDzI19kZIu7kUcOr+VEc7jgblb+5mbmaGnv98dq7o/aOEDnHee+xwagr6+xRZ+T4+76efn3Rb8v4MHF65VeF8wTaicGiivVI9d/d9ax6j+fvDg0v1xiFG9J0EWNrPws0vY9Zwbgn7cWlZmDYsybHm3vIm4RnXKKa378MfGFvmf2z3PU1OLDcWentC4Qgo0HOcI8ejYlP5Nz3p9nJP1u7JOX1q3XvetW6871o64nleVWiemXR/+2NjCjd/Ewp8r9S3tbfX1LR2fCPfaMmbhi9Z76mSAgYEBnZ6eTiQv8+E3JziZw/eiGO/Uexdh2I/biJ4e+Id/WGzhB9mwAfbtgyOPhKeecpbtT37irN4NG+rOeImb8XG45BJnXC5blsz1aXT/tLvQKHgsiHmRUjizgA//scdg/52TrDhviIcegp/fNMGvsJPDeI2X11/EOaOLffgccwy89a0Lf4+OumO16cN/aecL7H7tKPouGmbNxvYrKyI7VHWg5s56T4IsbDZLJztMTTmDo2q0Vo2WZhZbYj2B0VHVN79Z9fjj3d+1umy1/LiNtqCvO+zDD1uPNahX9yTOSZI9sFZ6yJ2Wp53egU/CdQp2BuL0AvjwNtDAwk9d1BttJvjZYWRkaW+1WaP07iobHVVdudK5SILT6kZHa3eTwyoxNbV4wLOvz3Xv169XXbdO/3P123Xv0afoT9++tiVBb0S9uvs6J1lyqcUpymm5W2vVKYlz7uNcNhJ8G7TNAHl0J33wg7B2beMyexkMr56cJ55w0xgBnnsOHnnETaubnHTTDMOIOJdMcCljNRpZcOpdpUCL3A4/gXvXQJRLUa/uPs5JOy6SJNpW3CtY01ikVKtOSax4jn0Vbr0nQRa2Ilj4eRkwrhrH1fGlVsrZcd2qplSwH13LFVMdVK1l4Y+ONjXHghabbys1Tgu/1bIm2bay1OPwRVp1ipovZuFnlzSmhHa6LP0LX4A774R3vGNxwKd6x6plndXMuxoqoDqtb/t2d0JEFqbMiSx5UfJ8bx+l4IFuuw2WL4dPf9ottmlyDoJW8tatfi2repapD4u1VSswybaVdryfOEirTrHmW+9JkIXNLPzs5Fdrfc6yZQsz11o5Vs0l8dXR4FqDpj09Cwtjli/XR96xQX/ESt3BKXojI/rZ9Z2frLR8tL5opazBa93fH0s0ASODYBZ+dunU4uvUN9up1VdvBT44dW7lWD++dpxvz2yihzlmZpa5JfGrJmtPkxRxcwq3bj00vfEVBln7PjhwwCXp3w7vLfuzkrNqpda61q2Utdq2JiZcVIIvfhFuvz2j02iNRDDBzwDtCk2UFyh0OihU/b8DB5zY9/S4hYuq7uFR91iByftnfeMShNnKSsQDnMbk4gNX6etzsUpCsUQGgQsvhLGxhXw7dVP4HgwMBrD0GQIl6ssyqgPFs7P5XEmexwkNWcaL4IvIWcB1uJeY36yq19RJ9+vAA8D/UNWv+si7iETxzXYqdPVW4FfLMzTkQuhycSj4VFWtenromXdL5hXoKZXckvjBQbjvvpaDRg0POyvVh6/dl0VfLsP7Aj2P6uQhH8f24YdP+30BnZKFN5n5IjMPrnq+nlY3nMj/EPgloB94BDi5Trrv4F52/vutHLsIPvxOyNTMnuCMmqAvvuo0rjrKe3pcKIPqZ4R57lnztYfXc1UnD/kgy/P2G0U38JFPXPP7424/4eMnfb8Ssw9/LbBLVZ8GEJE7gHOBnaF0lwJ3Ar/uIc9C08hKT8ySqPowtm1beFNFMJhU9e+gaRnwx0cpXNZ87WGvVF+fPyval+vJ9zmrZX2DX4s8jp5J3L2GWsfPUnBGH4K/Eng28H0P8O5gAhFZCfwu8H6aCL6IbAQ2AqxatcpD8bJJVGGudQMHG1upBB/7mOeQ2h/5CHzjG/CGN8Dzzy+NBFgqoZV46fO9fZSGh10BMtGXjY+gVwoWn3MfD+CsPeCg/ntZfQpbHIuu4hbfWsfPlEutnunf6gb8d5zfvvr9o8AXQmn+H/Ceyt+3UXCXTlxdvHBIcBFPx69GeGwUFXL5ct01OqZjpRH9S0YWv4kocJgku9JZOE5mXG+eqVW3PNQ37jI2WnCXlBuSmF06e4DjAt+PBfaG0gwAd4gIwBHAOSIyq6pf85B/Zqln3cVlZVQtiWpIcFXnZrjqKre1msdj4+VDUQLXrMF1G159tXbiUunQi1y/MjnIZ4A5oBSaQZNGV7qT4/sMW5Clrrxv6lnf7VjkaQxkJhGq4YIL3Gewp5eZXlq9J0GrG84t9DRwAguDtr/aIP1tFMDCb2RJxGllVEN+B8Nw13qJUb1/3rt+RF9lmR6kpK+wXPeuH1ncbQhHUAsMvo6NLYzLhvNrZQCumRXUaL+vAT6fYQvyYPGmRTeem6zUiTgtfFWdFZFNwLdwM3G2qeoTIjJS2X9T1DzySCPrzqeVEbaSqtvwsLPq77mn9mtKax7o9NN5y6uvIWjl7T0zPL8Xjq46IFVd6IJjjnFmcCj42GWXubxKJTc+G8yrmR+zmWXdbL8vP6nPsAVpBf7KA3H0ftKe+piLHl29J0EWtm618Jv9X6u+vmZ5NNxf7QqccooLObxu3SHTdh50FtFXWO7eMtRCoaJY8FNTLsR8oxj7PnoIrdJu2IJusVCTxPf5y8L1yEIZVBtb+KmLeqMtz4Kv2r4AtdtgWhHB0BvzHKOjh14YvehF1729qqWSzvb1L32lnOeyh/+vmfspi/PRs7YeIG9UbQ4fMX7SeFFK3OsQOsUEPye022jbtfAfHVuYbTMf8sfPg7P0G7TWKD72VuoM7tWvPo8f/v8sWGCGw+f1SKJH7SO/KHm2SiPBt1g6GaJdP3QzH/HkJHz0tXEu1Fs48OphnHzxFMzPAlR89CE2bIArr6yZVyuzVzqZiTA0BL29zu8J8PDD9dM2On4r/ttc+FhzTLs+dJ/Xo5Pxkk5mdlXruHt3Z2VPPVxEvSdBFraiWfiqHp/+Y2P686OOX+S2mQ9Z9POgz3Cs7malPrthtOHh4uwyj4zUfiNhq7RqbWVhjnS30onFm3aPK0qPur/fzYQrldxnq26pJFxPmEunQExNuXe0Bl01IZGvbv9+5gZvA8RRixzl2O3cRGFxD44hRAzxU2g6FbI0H7ZRx8yq4w/tvhMi7odcI8E3l043MT4Ol1ziYuFWkMqnAgcpMcVvchiv8fL6izjnbzZS24GzlDinGEY9djuusLBbaHJyIeTz/Dxs2gRr1pirp106nRab5oKkdttduI7VqCFzc+25dmotzEqMek+CLGxm4bdIdbpDvQVSIvrjd6zT0/qnunbAMsrgW2/vwqnq6Ulmhkc3UgTXWL0eYiv3VVIuLMzC70Kq0Sp37oR//MfFgczARa9829vcNjrK4YODbMlKTO4OaTQo2KmlODgIN9zgLPu5OfeSrbzEi/eNr4B+5TJs2ZLfdtaIcDtrp5eQhUkDJvh5ZHwcLr544X2DQUTctJfrr1/yIu/MxPPogLhmN5TLLmLz9dd7idycW9KIRZRHOn3dJGQjaqYJft4ol5m/+JPI/Pwh//whAoHMuuouI76l+N0sTu3g6/xmwYqNCx+vm0w71IYJfs740cQkx87P0VP5XnXiSKkEN964xKpvRtrxR1olDuuom8WpXZKORZRHfLSXtHvZJvg547sM8SH66WEGgDmEv+05lxNvHGXNxvZaUp4s3Diso24Wp3bx+WattK3YuOiG9mKCnzNOGh7kzJsn+Z+z7vVKEwzzfQb53/thTZvHypuF69s66mZx6gRf5zdtKzYukmgvcfe4TfCzRpMrPjgIv/rxQT5508K+vlJn1kY3WCxRSUKc8uI2M5oTZ3tJosdtgp8lWrziw8Nw++1uwVBPj5th0uk0OrNw4yVPbjMjXZLocZvgZ4E2IzL5FOpu7X5nhby5zYz0SKLHbYKfNkETsFRyc+ih6RVPS6jNPdGY8PnJg9ssyWtq7ac+SfS4TfDTYnwc7rwTXve6BRMQ3Dz6VasyeUeYe6Ix9c5Plt1mSV5Taz/NiduQM8FPg/Fx+IM/WPje1+c+qxGZMnoXmHuiMfXOT/gmzpKVm+Q1tfaTPl4EX0TOAq7DvcT8ZlW9JrR/A3BF5et/Aher6iM+8s4ld965+Ps73wnr12dDARqQB/dEmrRyfrJm5SZ5Ta39pE9kwReREnADcAawB3hQRO5W1Z2BZP8GnKaqPxaRs4Fx4N1R806Ttqy0cOLzzoNvf3th/0UXtb1CNg2y7p5Im1bOT9as3CSvqbWf9BHVJS+6a+8AIoPAVar625XvVwKo6pY66d8IPK6qK5sde2BgQKenpyOVLw5attIqES3nb7kVZmfR/n5K91USV334552XC7E3/JA1C9/oPkRkh6oO1Nrnw6WzEng28H0Pja33i4C/85BvarRkpVXubH3tNUQVAQ4emGHPxCTHDw46kTehLxxm5Rpp4kPwlwRtpMb7sQFE5H04wf+tugcT2QhsBFi1alWkgsU1ONbUFzk+Dp/7HFTEXnExbw7Sz3cZYthfUYBsDQIazbG1D9mkCPeRD8HfAxwX+H4ssDecSER+DbgZOFtV99c7mKqO43z8DAwMdORvqr4b5FbnSfHedW5opYVm4CjCAfq4TT7GHX3DbBn225LMRWAY0SnKfeRD8B8EThKRE4DngA8D5wcTiMgq4C7go6r6rx7yrEv1wr322sILoOIYHKtppZXLzrIPIL98Ik9dPsGP9w+yZch/I8raIGDeKYKVZyylKPdRZMFX1VkR2QR8Czctc5uqPiEiI5X9NwF/DKwAbhQRgNl6gwpRqV644Fh0b2+MU8CqCrFiBVx2mXvSBLn8ctZsHGw7kmWr2FQ3f/i08uzBkS+Kch95mYevqtuB7aHfbgr8/XHg4z7yakb1wh044N4AKLJY/INEvimvuMJZ9KouLEL1XdgicOKJcPnlsQ/M2iCgP3xZeUVxD3QTRbmPum6lbfXCXXUV3HOPE/25uaU3b+Sb8oor4NprF77PzTnRL5XcAScmEms1NgjoB19WXre7B7q191KE+6jrBB/cRbvqKrj//vo3b6SbslyGz39+6e8f/CCsXZvandDOjditN20UfFl53ewesN7LUvJ0L3Wl4EPzmzfSTTk5udRPVCrB6GhqV7ydGzGYtrcXLrww0yF8EsWHldfN7oE89F7qCXAcwpy3B2DXCj40vnmDN+WKFe6z+nuV6vROCAni0BAcdtjCAO173wvXXJPqlW7nRgymnZuDsTH3QpWsN9Y80a3ugaz3XuoJcFzCnIcHYJCuFvxmVC9MvQYyNOR+Bzen/777Kv+TQROunRuxmrY6dVU1H43VSJ8MNv1F1BPguIQ56w/AMIUWfGjcQE6dKfNRnIn/fw8MMzk5uNBIKiZcuQyTW9Jv/K3eiNVu7dat8NBDsG2bq3seGquxlDT8x1nuvdQT4LiEOesPwCWoama3U089VeNmakp1+XLVUsl9Tk253x8dm9JX6dd50HnQV1mmj45NtfS/WaVWecfGVM88031mhakp1c2bs38+0yZv7S8p6rWforQrYFrraGrhLfyaT+hymTV3XoVy8FCgoGUyw5r9k8DCIzxv/rtweScmnO9+ZsbNaFqzJv3y520QLE3y1v6Sol4PJMs9k6QovOBDqCFUFefAASQQA05q9APz5r8LlxeyJxgmYq2Tt/ZnpI8JfpXwO2bn56GnBwYG4F3vqjlvsZ7/LqvzcsPlhQULPyuCYSLWOrnzHxupE/kFKHGS2AtQwqtm+/qc4HfgUyiX4X3vWxCsQzN7MkoWH05ZLJNh5IW4X4CSb2qtmo3wjtmJCRfHB9xnghEWOiKLfs0slskwuoFiC3657GIwhHs5OXnHrGEYRjv0pF2AVCiX4eKLnQV/zz0LES57elx4hAhiPzzsXDki7nPY9+utDMMwOqR4Fn6tN6T09MAHPuCs/Yi+hOqiLfNBG4aRNYon+OE3pIjAsmWLxD7qoKH5oA3DyCLFEfzgm6mq8/5qhIq0hT+GkTw2MysZiiH4IRX/4aVbefbh/aw4b4g1Gxe3Llv4Y9TCBCk+zMhKju4W/Opdunv3IRXXAzPc/uf72axX0n8/3BsKJ2ALf4wwJkjxYkZWcnSv4FfjGx886F5O0uuqOtvTz3fmhpibr924bPWiEcYEKV7MyEoOL4IvImcB1wEl4GZVvSa0Xyr7zwF+DvwvVf2Bj7xrUi7DZZctBLOfnYXf+R1Yu5Z/XjHEDy4bpNSgcdmgqxHEBClezMhKjsiCLyIl4AbgDGAP8KCI3K2qOwPJzgZOqmzvBv6y8umf8JtLqhx1FFx5JWtwbhxrXEarmCDFjxlZyeDDwl8L7FLVpwFE5A7gXCAo+OcCE5VYzQ+IyOEicrSqPu8h/8VMTCwV+2XLFq2AssZltIu1GaMb8LHSdiXwbOD7nspv7aYBQEQ2isi0iEzv27cveunWroX77qPMIFu2uA6AL8plvB/TMAwjLnxY+FLjt3AIzlbSuB9Vx4FxcNEy2y7N8LB7AW3V4bp1K2UGvc+ysJkbhmHkDR8W/h7guMD3Y4G9HaTxw+Cgi0n8p396KDZxrVkWUYnjmIZhGHHiw8J/EDhJRE4AngM+DJwfSnM3sKni33838NNY/PdVQg7XOGZZVI954ICLzrBiRfRjGoZhxElkC19VZ4FNwLeAJ4GvqOoTIjIiIiOVZNuBp4FdwBeBT0bNtx2qsyyuvtqf62VwELZudVP85+fdLFDz5RuGkWW8zMNX1e04UQ/+dlPgbwUu8ZFXp8Qxy2L/fif283UWcRmGYWSJYsbD90TVrVMq2YIcwzCyT/eGVkgAW5CTLBbAzDCiYYIfEVuQkww2DdYwomMuHSMX2DRYw4iOCb6RC2y8xDCiYy4dIxfYeIlhRMcE38gNNl5iGNEwl45HshRMLUtlMQwjG5iF74kszSLJUlkMw8gOZuF7IkuzSLJUFsMwsoMJfoWoLpAszSLJUlkMw8gO5tLBjwskS7NIslQWwzCygwk+tV0gnYhklmaRZKkshmFkA3PpYC4QwzCKgVn4mAvEMIxiYIJfoRMXiEVvNAwjT5jgd4jNdTcMI2+YD79DbK67YRh5wwS/Q2yg1zCMvBHJpSMibwL+GlgNPAN8SFV/HEpzHDABHAXMA+Oqel2UfLOADfQahpE3xL1fvMN/FrkWeFlVrxGRPwTeqKpXhNIcDRytqj8QkdcDO4D1qrqz2fEHBgZ0enq64/IZhmEUDRHZoaoDtfZFdemcC9xe+ft2YH04gao+r6o/qPz9H8CTwMqI+RqGYRhtElXw36Kqz4MTduDNjRKLyGrgncA/NUizUUSmRWR63759EYtnGIZhVGnqwxeRe3D+9zB/1E5GIvILwJ3AZar6s3rpVHUcGAfn0mknD8MwDKM+TQVfVT9Qb5+I/LuIHK2qz1d89S/WSdeHE/svqepdHZfWMAzD6JioLp27gQsqf18AfD2cQEQEuAV4UlX/PGJ+hmEYRodEFfxrgDNE5CngjMp3ROQYEdleSfObwEeB94vIw5XtnIj5GoZhGG0SaR6+qu4HTq/x+17gnMrf/wBIlHwMwzCM6NhKW8MwjIJggm8YhlEQTPANwzAKggm+YRhGQTDBNwzDKAgm+IZhGAXBBN8wDKMgmOAbhmEUBBN8wzCMgmCCbxiGURBM8A3DMAqCCb5hGEZBMME3DMMoCCb4hmEYBcEE3zAMoyCY4BuGYRQEE3zDMIyCYIJvGIZREEzwDcMwCkIkwReRN4nI34vIU5XPNzZIWxKRh0Tkb6PkaRiGYXRGVAv/D4F7VfUk4N7K93p8CngyYn6GYRhGh0QV/HOB2yt/3w6sr5VIRI4F/itwc8T8DMMwjA6JKvhvUdXnASqfb66TbiswCsw3O6CIbBSRaRGZ3rdvX8TiGYZhGFV6myUQkXuAo2rs+qNWMhCR3wFeVNUdIjLULL2qjgPjAAMDA9pKHmHKZZichKEhGBzs5AiGYRjdR1PBV9UP1NsnIv8uIker6vMicjTwYo1kvwn8NxE5BzgMeIOI/JWqfqTjUjegXIbTT4eZGejvh3vvNdE3DMOA6C6du4ELKn9fAHw9nEBVr1TVY1V1NfBh4DtxiT04y35mBubm3OfkZFw5GYZh5Iuogn8NcIaIPAWcUfmOiBwjItujFq4ThoacZV8quc+hoTRKYRiGkT1EtSM3eSIMDAzo9PR02/9nPnzDMIqKiOxQ1YFa+5r68PPI4KAJvWEYRhgLrWAYhlEQTPANwzAKggm+YRhGQTDBNwzDKAgm+IZhGAXBBN8wDKMgZHoevojsA37U4b8fAbzksTh5wOpcDKzOxaGTeh+vqkfW2pFpwY+CiEzXW3zQrVidi4HVuTj4rre5dAzDMAqCCb5hGEZB6GbBH0+7AClgdS4GVufi4LXeXevDNwzDMBbTzRa+YRiGEcAE3zAMoyDkWvBF5CwR+RcR2SUif1hjv4jIX1T2Pyoi70qjnL5pod4bKvV9VESmROSUNMrpk2Z1DqT7dRGZE5HfT7J8cdBKnUVkSEQeFpEnROS7SZfRNy207V8UkW+IyCOVOl+YRjl9IiLbRORFEXm8zn5/OqaqudyAEvBD4JeAfuAR4ORQmnOAvwMEeA/wT2mXO6F6/wbwxsrfZ+e93q3UOZDuO8B24PfTLncC1/lwYCewqvL9zWmXO4E6fxr4bOXvI4GXgf60yx6x3uuAdwGP19nvTcfybOGvBXap6tOqOgPcAZwbSnMuMKGOB4DDKy9bzzNN662qU6r648rXB4BjEy6jb1q51gCXAncCLyZZuJhopc7nA3ep6m4AVc17vVupswKvFxEBfgEn+LPJFtMvqvo9XD3q4U3H8iz4K4FnA9/3VH5rN03eaLdOF+GsgzzTtM4ishL4XeCmBMsVJ61c57cCbxSRSRHZISLDiZUuHlqp8/XA24G9wGPAp1R1PpnipYY3HcvzKw6lxm/hOaatpMkbLddJRN6HE/zfirVE8dNKnbcCV6jqnDP+ck8rde4FTgVOB5YDZRF5QFX/Ne7CxUQrdf5t4GHg/cCJwN+LyP2q+rOYy5Ym3nQsz4K/Bzgu8P1Y3FO/3TR5o6U6icivATcDZ6vq/oTKFhet1HkAuKMi9kcA54jIrKp+LZES+qfV9v2Sqr4CvCIi3wNOAfIq+K3U+ULgGnXO7V0i8m/ArwDfT6aIqeBNx/Ls0nkQOElEThCRfuDDwN2hNHcDw5VR7vcAP1XV55MuqGea1ltEVgF3AR/NsbUXpGmdVfUEVV2tqquBrwKfzLHYQ2vt++vAe0WkV0ReB7wbeDLhcvqklTrvxvVoEJG3AG8Dnk60lMnjTcdya+Gr6qyIbAK+hRvd36aqT4jISGX/TbjZGucAu4Cf46yDXNNivf8YWAHcWLF4ZzXHkQZbrHNX0UqdVfVJEfkm8CgwD9ysqjWn9uWBFq/z1cBtIvIYztVxharmOmyyiHwZGAKOEJE9wJ8AfeBfxyy0gmEYRkHIs0vHMAzDaAMTfMMwjIJggm8YhlEQTPANwzAKggm+YRhGQTDBNwzDKAgm+IZhGAXh/wMxSby+lsnn6AAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import os\n",
"# os.chdir(\"C:/Users/kevin/Documents/GitHub/DS_course/ipython/\")\n",
"os.chdir(\"../\")\n",
"import course_utils as bd\n",
"import imp\n",
"imp.reload(bd)\n",
"\n",
"#Generate Y, and X, where E[Y|X] is a 3rd order polynomial\n",
"betas = [0, 2, -2.5, 1]\n",
"n=200\n",
"sig=0.2\n",
"sp=20\n",
"\n",
"x_init = np.random.uniform(0, 1, n)\n",
"e_init = np.random.normal(0, sig, n)\n",
"\n",
"dat = bd.genY(x_init, e_init, betas)\n",
"dat = bd.makePolyFeat(dat, 6)\n",
"\n",
"#Plot the data vs. the real curve\n",
"plt.plot(dat['x'], dat['y'], 'b.')\n",
"plt.plot(dat['x'], dat[['x','x2','x3']].dot(np.array(betas[1:])), 'r.')\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEICAYAAABcVE8dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABdJklEQVR4nO3dd3xb1d348c/RlmV5752dOMNOYmaAAqUlbNqyV0LLKruMEtpfC0/HAw+UUihlFcp4SuFhjxBooRAoJYEEyE6cnXhb3pKtrfP7w9at7XjItjyinPfrpVcs3at7z3Xkr84943uElBJFURQl9unGuwCKoijK2FABX1EU5RChAr6iKMohQgV8RVGUQ4QK+IqiKIcIFfAVRVEOESrgKzFPCLFZCHH8eJdjvAgh9gohThrvcijjTwV8ZUx0BR23EMIlhGgWQrwrhMgfi3NLKWdLKVdG+7hCiKVCiGDXNXV/5ET7XIoSDSrgK2PpDCllPJAN1AF/HOfyRMMqKWV8r0f1eBdKUfqiAr4y5qSUHuBVoDj8mhDiNCHEN0KINiFEhRDi7u7vEUJcJoTYJ4RoFEL8onszhRDCKoR4ruvOYasQ4qdCiMpu7+2+791CiJeFEM8LIZxdzT1l3fZd0FUOpxDiFSHE/wkhfjPUaxRCTBFCNAkhFnQ9zxFCNISbloQQl3eV1SmE2C2EuLrbe48XQlR2XUe9EKJGCHG2EOJUIcT2ruP+rNv+dwshXu0qq1MI8bUQoqSfcumEEMuEELu6fpcvCyFShnp9ysFJBXxlzAkh4oDzgdXdXm4HLgOSgNOAHwshzu7avxh4FLiYzruDRCC323vvAoqAycB3gEsGKcKZwEtd53obeKTrPCbgDeBZIAV4Efje0K8QpJS7gDuAF7qu9xng2W5NS/XA6UACcDnwYPjLoUsWYKHzOn8J/LnruhYCxwK/FEJM7rb/WcArXeX+G/CmEMLYR9FuBM4GvgXkAM3An4ZzjcpBSEqpHuox6g9gL+ACWoAAUA3MHWD/PwAPdv38S+DFbtviAB9wUtfz3cDJ3bZfAVT2Ond437uBD7ttKwbcXT8fB1QBotv2z4Df9FPGpV3X0tLtsavXPm8DG4ENgHmA630TuKnr5+MBN6Dvem4HJHBEt/2/As7udk2ru23TATXAsX1c/1bg2932zQb8gGG8PyPqMfoPVcNXxtLZUsokwAxcD3wihMgCEEIcIYT4WAjhEEK0AtcAaV3vywEqwgeRUnYAjd2O22N7r5/7Utvt5w7AIoQwdB2nSnZFwgiPtVpKmdTtMaXX9j8Dc4A/Sim94ReFEKcIIVZ3Nc+0AKfyn+sFaJRSBrt+dnf9W9dtuxuI76ucUsoQUNl1Pb0VAm8IIVq6zrsVCAKZg1ynEgNUwFfGnJQyKKV8nc5Ac0zXy3+jszacL6VMBB4HRNe2GiAv/H4hhBVI7XbIHtuB4Y7+qQFyhRCi22vDHkkkhIin807laeDucFu5EMIMvAb8Dsjs+hJcwX+udzi0cgohdHT+PvrqPK4ATun1JWWRUlaN4NzKQUIFfGXMiU5nAcl01jChs9miSUrpEUIcDlzU7S2vAmcIIY7uamf/L3oGx5eBO4UQyUKIXDrvHoZjFZ1fQtcLIQxdZTx8mMcCeAj4Skp5BfAunV9iACY673IcQEAIcQrw3RGcB2ChEOL7XXcqNwNeevaRhD0O/FYIUQgghEjvuk7lEKACvjKW3hFCuIA24LfAEinl5q5t1wK/EkI46Wyzfzn8pq59bqCzo7UGcNLZ6RluIvkVnU0Ye4AP6fyC0JpPIiWl9AHfB35EZ3v8JcDyQY51VB/j8A/rCqKL6WyaArgFWCCEuFhK6aSz8/RlOjtNL6Lz7mYk3qKzI7wZuBT4vpTS38d+D3Wd6x9dv+vVwBEjPLdykBA9mysVZeLraippAaZJKff0sf3HwAVSym9F4VxfAI9LKZ8Z6bFGS9cQ1qlSysFGJymHOFXDVw4KQogzhBBxQggbnW3fG+kcfYIQIlsIsahrjPkM4FY6h1cO5zzfEkJkdTXpLAHmAe9H5yoUZXwZxrsAihKhs4D/pbPtfi2dNfjw7akJeAKYRGfN/yU6x+0Pxww6m1rigV3AOVLKmuEXW1EmDtWkoyiKcohQTTqKoiiHiAndpJOWliaLiorGuxiKoigHja+++qpBSpne17YJHfCLiopYu3bteBdDURTloCGE2NffNtWkoyiKcohQAV9RFOUQoQK+oijKIUIFfEVRlEOECviKoiiHiKgEfCHEYiFEuRBipxBiWR/bE4UQ7wgh1nctKXd5NM6rKIqiRG7EAV8IoadzibRT6Fw96MKuJem6uw7YIqUsoXM1nwe60twqiqIoYyQaNfzDgZ1Syt1d6WVfojPvSXcSsHctLBEPNNG5NJyiKIrSTUNDA/v372c00t5EI+Dn0nMZuEp6LjANnYtEz6JzBZ6NdK7dGerrYEKIq4QQa4UQax0ORxSKpyiKcvCor6+nqqqKnguvRUc0An5fper91XQysI7ONTZLgUeEEAl9HUxK+aSUskxKWZae3ufsYEVRlJjV0dGBzWYblWNHI+BX0nPdz77W0rwceF122knnykQzo3BuRVGUmCGlpKOjA6vVOirHj0bAXwNME0JM6uqIvYADl2vbD3wbQAiRSWfO8d1ROLeiKErM8Hg8hEIhampqRqUNf8TJ06SUASHE9cDfAT3wFynlZiHENV3bHwd+DTwrhNhIZxPQHVLKhpGeW1EUJZZ0dHQAEBcXNypt+FHJlimlXAGs6PXa491+rga+G41zKYqixCqn0wlAUlLSqBxfzbRVFEWZIFpbWwGoq6ubsMMyFUVRlCgIN+mYTKYJOyxTURRFGSEpJT6fD+hswx8NKuAriqJMAF6vV2vGUQFfURQlhjU2Nmo/+/1+QqE+kxGMiAr4iqIoE0Bzc7P2s8PhUG34iqIoscrlcmk/p6WlqYCvKIoSq7xerxbkbTabGpapKIoSizweD1JKLcjv27dvVM6jAr6iKMo4a2pq0n4WQpCQkKCadBRFUWJR94AvpZzQ2TIVRVGUEWhra+vxvLq6Wg3LVBRFiTUej0ebYRuWmpqKThf98KwCvqIoyjjq3pwTlpaWNirnUgFfURRlHHWfcNWdGpapKIoSQ0KhkFbDD4/K0el0bN++XQV8RVGUWOJ0OgkGg8B/avThzlo1LFNRFCWG9NecYzabVQ1fURQlloQzZPYekdPR0aGGZSqKosQKv9+vrWFrMERlefFBqYCvKIoyDlpaWgbcrpp0FEVRYkRDQ4P2s9/vP2C76rRVFEWJAVLKHgG/r9q8Xq+P+nlVwFcURRljLS0tBINBdDpdv4E9EAhE/bwq4CuKooyxqqoqoLMW319b/YRtwxdCLBZClAshdgohlvWzz/FCiHVCiM1CiE+icV5FOVRJKQkGg/h8PtxuNx6PB4/Hg9fr1SbyKBOTlFIbjjnQYuWj0aQz4rFAQgg98CfgO0AlsEYI8baUcku3fZKAR4HFUsr9QoiMkZ5XUWJdMBjE5XLR0dFBR0cHbrcbr9eL1+s9ILtib0IIjEYjFotFe8THxxMfH6/lWvd6vdqXRfgcgUCAQCBAKBRCSql1HOr1egwGAyaTibi4OBITE7HZbBiNxlHpXIxlDocDKSV2u10bltmXUCgU9aAfjcGfhwM7pZS7AYQQLwFnAVu67XMR8LqUcj+AlLI+CudVlJji8XhoaWmhtbWVtrY22tvbtW1CCKxWK2azGZvNhslkwmAwoNfre0zaCYVCWtD2+/14PB5aW1uprx+dPzkhBDabjeTkZNLT07Hb7eoLYBD79+8HIDk5ecCA7/V6MRqNUT13NAJ+LlDR7XklcESvfaYDRiHESsAOPCSlfD4K51aUg1YoFKK1tZWGhgaamppwu91A5yQcu91OWloadruduLg4rFbrgIFUSonb7aa1tRW/309bW5tWY4+WcC3fYDD0aFJyuVy4XC4qKirQ6XSkpKSQm5tLUlKSCv69hO/aTCZTjy/0vpjN5qifPxoBv6//0d69DQZgIfBtwAqsEkKsllJuP+BgQlwFXAVQUFAQheIpysQRCoVobm6mvr6exsZGAoEAOp2OpKQkcnJySE5Oxmaz9RsopZS0t7fT0NCA0+mkvb0dr9c7Kh18vQWDQYLBIF6vt8frRqMRk8lEIBDA6/XS0NBAQ0MDBoOB7OxsCgoKol5TPVjt3bsXgKysLKqrqwfcdzS+LKMR8CuB/G7P84DeV1IJNEgp24F2IcSnQAlwQMCXUj4JPAlQVlY2+p9iRRkDTqeT2tpa6uvr8fv9GAwGUlNTSUtLIyUlBZ1Oh9frxePxUF9fj8/nw+fz4fV6taAeDAZHFNiFEFoTkBACKSVSSkKhkNZmHwm9Xq8FeL/fj9/vJxAIaO/X6/UIIQgEAlRUVFBRUUFKSgrTp0/HYrEMu/wHOyklNTU1AGRmZmpNO30RQoxKuoVoHHENME0IMQmoAi6gs82+u7eAR4QQBsBEZ5PPg1E4t6JMWMFgkLq6OmpqanA6nQghSE5OJiEhAZ1Oh8fjoaamht27d+PxeEa9li6ljEoTTzAY1JqfLBYLCQkJeDwebV1WIUSPkSdCCJqamli9ejVJSUnMmjVrVJorJjqHw0EgEMBut1NbWzvgvlJKPB5P1L8gRxzwpZQBIcT1wN8BPfAXKeVmIcQ1Xdsfl1JuFUK8D2wAQsBTUspNIz23okwUgUBAG+nS1tZGS0sLHR0dPfaRUtLU1NTnknYHq/BwUOhsc46Pj9dGFOn1eqxWKx6PR/uiaWlpYdWqVaSnpzNz5sxRGXo4EUkp2bVrFwBFRUVs335A48YB/H7/xAv4AFLKFcCKXq893uv5/cD90TifooynUCiE0+mktbUVp9OJ0+nUgl5/dDodZrMZo9GIlFJrshmLtvexEh4yCv9pPnK5XAAkJSWh1+u18ecOhwOHw0FhYSFFRUUx37nb2NiI1+vFZDJhNpsP6Afpzr5tG8a2NnRlZVEvx9jk5FSUg1i4ozRcO29ra9OaLLoPjezejGGz2UhISMBoNOL1enG5XLjdbq0pJNaFv9TCwpkh7XY7iYmJ1NbWEggE2LdvH5WVlcyePZuUlJRxKu3o6l67z8/P12bZ9qfo2WeJ37GD0FVXRb0sKuArSh/CQyYdDgcNDQ1a8ApPOgqFQrhcLm3MO6BNbPL5fDidTq2DTvmP8B2RwWAgJSWF1tZWgsEgGzZsIC4ujnnz5sVcx25jYyNutxshBGlpaaxZs6bffeP27yf1iy/Yc/nlZE3QYZmKEhOklLS1tVFXV4fD4cDv96PT6UhMTCQpKQmPx4PT6ezRNq/T6TAajT3GoyuDCwQCB/RldHR0sHr1atLT05kxY8aYLQoymqSU7N69G+gcmdPU1DTgSla5r71GwGDgOYuFmzwebVZ0tBz8v1FFGSGv10ttbS21tbW43W4tyBuNRtrb27V1R/tKdBUKhQZsjx0TUmJoa0PU1uKrqCDU0IBsakK0taFzuxFeL8Ln6zFhJmQ2I81mpNWKTEpCl5GBISsLfX4+/owMpMk0bpcTbt/Pzc1lypQpByz/dzCprq7WKgg5OTls2bKl330NTidZ//gHb8bH8+Sbb3LHAw9EvTwq4CuHJCklLS0tVFVV0djYiJQSm81GUlISHR0dfS4u3VdSsnAuGb/fP6odsCIQQO7YgXfDBsT27VgqK7E3NpLmcpHh9TLQzb+36xGuVwrADAzUcOIwGKiPj6cpMxNPYSG6khJMRx5JcAzb2auqqqiqqsJisZCXl0dqamrUa7yjyefzsXv3bnQ6HTabjba2tgE797PffRe9x8N/eTycv3TpqHzRqYCvHFLCY+Orqqpob2/XRs/4/f5Bp7r3pa+VikasqYmOf/8bvv4a+65dZDU0UOh2073OXQdU6PVsjIujNT2djuRk/GlphDIyEKmp6NPS0KemYrDbMVutGAwGdDqdNtHK7/fjc7sJtLXhr68n5HCAw4G+rg5LfT0JLS1kOp1Mb2khs7wc/vEPAGqNRvZlZNAydy6mb38b3fz5yFEeWunxeNi5cyc7d+7EZDKRmppKSkoKycnJE7rZZ/fu3VolISMjQ+u47YsIBsl94w3W2O1UGgyccsopo1KmifvbUpQo8vv9VFZWUlVVRSAQ0IYBhkKhQYdUjibp89H+2WfIzz8nYds2JtXXU9DtS6QC2BUXx1dFRXQUFRGaPh3LvHmkT52K2WwmDogbpbL5gTUuF/XbtuFdswbL5s1kVFZSXFXFEVVV8P77tArBluxsmg4/HOv3vocuSulQDAZDn5PEfD4fNTU11NTUIIQgMTGRlJQUUlNTiYuLmzDDO1tbW6mtrdVmJe/duxcpZb9DMtP+9S8s9fX8Gjj3qqswm82jMkdBTORxwGVlZXLt2rXjXQzlINbR0cHOnTsnzGSnQGsrbe+/j3HVKnJ272aW06kF7Coh2JaQQF1hId45c4hbtIi0GTNG/IcfCAj8fh3BoCAYhFBIIKUgJaVz5FFbmwGvt/s5JDodpKZ2bne79QSDAr1eotNJQiEvtWu/IPjJSjI2bGCBw0Fe1zs322zsnT8f4/nnY5ozZ0TlFkJgsVgGHMpqNBq1uyyLxUJqaiqpqakkJSWNW9t/KBTi66+/pqOjg1AopH0JhTv3+3gDC6+5hraKCmbpdLz48svYbDaOPPLIYY1YEkJ8JaXscxC/quErMcXv99PU1ERtbS0tLS3jPrHJ29xM63vvYfn8cwp272a2240RCACbTSY+KCrCVVKC6VvfwjZjPi6XmY42Ay6XkZwpbej1IbZssbN2bQodHXo6Ogx0dOjxePQsW7aV+Pggr7+ey5tv5uLz6bSH3y94++1/YzaHeOyxKbz+el6Pcul0kn/+s3MdoieemMKKFdk9tttsAZYv/wyA++6bwcqVPZewSEs7jldeOQaAS++ewddr4jH424lrd5HwmZfpn23nN3GnsfOII3jB+jC1zkxMphBmcwizOUhOjptzzukcj75yZTputx6zOYjV2rk9JcVHUVHnjN22NoP23t4V+HCwt1qtmEwmampqqKqqQq/XazX/lJQUTGPYCb1nzx5cLpcW6G02Gy6XC4PB0GfAT//kE+w7dnAdcObFF2Oz2QAGHM0zXCrgKwet8ISo1tZWmpubtdTA48nv8dD4979j/OQTCrZv56j2dtwk8g1zed/4Ax7OnE1jxjzak2fxw6vryMnxsPrvmTx45/RetWx49tkvKSzsYMuWRJ55ZhImU5C4uM6HxRLE59MBQZKTfUyd6sJkCmEyhTAaOx/h4Hj00Q2kp3u1GrpOJ+l+03DyybUUF3fmwZGy82E0ygO2h0KCQEAQDAoslv90YJeUOrEnSPx+gd8fj7M+RHutn2CLn7M+/pjXWUOl/kg64lIImJLx+fRMntyuBfznniti715bj2svK2vi/vs3AHDVVWXU1VkQQmKxBLFYQhxzTAO33NKZnuBXvyomEBBYrUHi47NJTTUyb56bBQv243A4WLkynaQkM5mZ8WRlxZOZaSMzU5CaGpX/8h4cDgcVFRVdv0tJbm4uLS0tGAyGA1JtQGfb/aRnnmGvzcYrPh9/+8EPtG2jkWFUNekoBw0pJS6Xi+bmZm2hkPFezk9KyZ5PN7LxFQeePSFCHRnUUch+CrjBeCdpU/bzr9xruf+fN2jvMRpDJCf7uOuuzRQXOykvt/PRRxkkJvpJSPBjt/ux2wPMmuXEag3i9wuEAINh4v6t9kVKSf0nn2D6299YtHMnGVJSqdPx5fz5mK+/HltREQCtrQbcbgMejw6vV4/brcNqDTJjRuechhUrsmhtNeJ263G7O+9uJk928b3vdSblvfnmUtraOrd3dHTus3hxLb/9bSNJSSkUFOQeULZLLmngnns8GI2J5OfHY7MJbDa0x49/DFdeCS0tcP31na/FxXU+rFb47nfh8MOhrQ1WrJAI4aW9vZG2tnpMJklurof8fAvJyVmsXbsPgyGE0SgxGEIYDBK9XiJE58icGb/7HWcB1vPP55prrgE6m7OOOeaYYTXnDdSkowK+MqF5PB6amppobm6mublZ68jrncpgNAWDUFUVx/79nY/dO2D/Jj9nyoe4uuVJnMF5LORrAKzCRbq9icxJei79UT1z57bS2mpkx454UlO9pKb6sNsDBzRNHKzCaZYH429vp/G55yh6/32OdDrxAJ9kZ9N2xRWkn3hi1MpjMBiwWCy0tbnQ6TrvVqqr7Vit6eh0iTidIerrO0hObmTy5CY8Hh1//eskgkEzgYAZn8+E12vgggskS5YYqa4WHHsstLeD2w0dHZ2fh/vu87JkSStffeXh1FMP7Ki+/fZyli1L56WX9nLddQsO2H7XXZs58ahKOs5/gu+1/i8egthsVvR6EELy5JNbOPfckmF1QquArxw0pJQ4nU4aGhpobGzUhkqGhxUOtpbryM4NtbUWdu2KZ/duG5Mnt3PMMQ046uG884/X9suihuls5xoeZnLKSvbPKaN6wWVMOr4Ae0IwZoL5aGldtQrzn//MsXv2YAM+s9vZ+4MfkHPJJeiiODKlry8jo9FIeno6GRkZmEwmLQleeJZ070pE71xJgUAAjyek3XH5/TqcznQ8HkFjYwfx8anU17eyaFEiUEF9veDzz9MIBASBQGffSjAoOPbYBhZ98Syhxz7kJK4kY+ERFBRMBjo71S+7bB+nnbZwWGmkVcBXJrRwSoP6+nocDocW1K1W66jNZA0GweUykpjoR0q47bYStm2z09HR2a0lhOTEone53n0lJbV1rOICprEDs3EPtdMy8Z54IvGnnIKIG61BkcMXrhX2/tvW6/UYjUb0er22PGE4wIVfG0u+2lrcDz7IUWvXkhEK8Y3JxIbTTyfnqqswDjOPzFCuQwihLSWZmZmJyWTSFnUPP8KrfAWDQfR6vbbMo9Vq1RaEr6ysZM+ePWRlZeFwODAajYMO9TU7HBy2dCmfBgJclZ/PE088cUDzzRFHHDGsiWYq4CsTksvloq6ujvr6erxer7bKT/fVk6KlqcnIli2JbN6cwLZtdsrL7Uyb5uKhh9YBcM89U/HX7aOo4UOOdPyTc31fkUwHLiHYmJFB42GHYT37bPRTpkS1XMMhhDjgVj+8etVwhIcv9tdEptfrSU5OxuPxjEquIOl20/bII8z/xz8oCATYotfzxXe/S+6NN2KKYiK18MQzKWWfTYI6nY64uDgSEhJITk7WVvUKBAIEg0Ht3/AjFArR3t6Oy+VCr9dHvmqYlMz92c+IX7OGWcEgNz74IKWlpQfsdtxxxw1raKkK+MqEEO50ra6uprGxcdSaZ6SEioo4du60ceKJDgCWLZvLF1+kYjCEmDrVxYwZbRSkbmFe1e9J//prFnYbS77DamXvjBkEv/tdLN/+NozCkL5I2r51Oh0JCQlkZmZit9uxWq0RdeKFJ5OF0zGHH+ElFEdSk4+PjychIYFAIEBra2t0774CAZxPPcXsN99kstfLFr2ez086ifybbsI8xJpuOO9RJP08kfZDdN8fDryDilTGhx9S/NvfcptOx1fHHcddd93V536jMQ5fBXxlVASDQdrb23E6nbhcLlpbW3G73aM2Lr6+3szq1amsW5fEunVJNDd3Bum33vqMhIQAW7Yk0OF0Eb/9/7D/+59M272bBX4/BqBVCDbn5NByxBFYzz4bkZ8/8MlGkU6nIysri7y8POJGqbkoEAjQ3t6u1U7D/0cTJRbIQAD3c88x85VXmOz1slGv5/PFi5l0ww2YRnFpRCGENl4/EAgMaZ3fSBlbWjhsyRI2ezycHB/Pn//yFxITE/vc95hjjhlW6ggV8JVRJaXsmiDTpq0C1d7ePqoBpLXVwFdfJbNgQQtJSX7efDOHhx6aTlqal/nzmykpaWHOnCaoXEno/ffI3riRw1pbSaUzidi2+HgqZ89GnHwypuOOG/V8MH3pXrNMTk4mOzubtLS0cZkhGgqFeizPGGntPZxwzmAw9Jk+ekSCQTzPPkvxyy9T4PPxtcHAF6efzpRrrhnVwN/bSGv0Gikp/q//IvnTTymVksvvv5+yAVa1Ouyww7RJWEOhAr4SVVJKOjo6aGhooLm5GafTqTUTCCG0ZfyiOQkqFIIdO+JZvTqVL75IZds2O1IK7rhjK4sX19HaasDpNGIJbKP17bdI/PJL5tTUMLPrlr5er2dbYSHtixYRd+aZyLS0qJVtKLrniLFYLGRlZZGVlTUhF/3weDy0trZqK30N9P8phCApKUlLahbOA99X1tGhEsEg3j//mTmvv06O388qo5Gvf/ADZvzoR2OaPG2oTT+95b/4IlOefJI7gD3nnsu111474P6LFi0a1uQrFfCVYeveNBMeJjlWa7H6fDra2gykpfmorzdz/vlHIYRkxgwnRxzRxGGHNVGYW03jP1ZgXLmSybt3s8DrxQS4gS2pqdTPn4/59NPRzZvHWI6X7N4paDQatU4+vV5PRkYGmZmZJCYmTphkX4MJ38U1NTVRX19PW1vbgPsbjUYSExNJTEzE5/NRVVU18nkTXi+Bxx6jZPly0oNBPjKb2XLhhcy65JIJvxh6yqpVzP35z3ldr+en+fk89vjjg6Z7UAFfGXU+n09LU9Da2tpnymCTyaR1IEZ7nVa3W8/q1Sl8+mk6X3yRwsKFzfz615sB+PTTNIpnNeLb9BGhDz4ge/NmFra1kdT13m1xceyfPp3Qt7+N9TvfQY7RbX84xXL3IaQGg4FQKEQoFEKn05GamkpGRgYpKSkTPjhFIpxmev/+/X0OQRRCoNPptDs/vV5PXFwcbre7zyyYQyHcboIPPcT8Dz4gORTiPauVXZddRvF55415c1gk80Pi9u2j9Mc/Zpvfz2l2O/f96U9kZ2f3u394+OcRRxyhZtoq0RUeOeNwOGhqaupz2F14/dH4+Hgtb81ozHJ97LEpvPlmDj6fnqQkH4sWNXDCCQ6KzP/C8+67pHz9NXMdDvK6PrOVBgM7CgpoP/pobGecgcjIGOQM0WWxWAiFQtofe/cx4Hq9ntTUVNLS0khNTY2JIN+ftrY2tm3b1m/bvdFoxGKxEAgEolo5EE4n8ve/Z+Enn2CTknfj4ti3ZAnF55wzYVbJstTUMPcnP6HD4eBYs5nbHn6YqVOn9rt/YmIira2tgKrhK1HU0dFBbW0t9fX1Wg0tPO44FAppMxLNZjN1dXXR64jrEggI1q5N5tNP07n55h2YTCFeey2XqiorR84pJ7f8aRLWrGZ6RQUzumqETUKwOTOT5gULsJ5xBsYZM8asmUan0yGE6DGksXdnXlxcnJahMTExccIEnbHS2NhIeXn5gLXdcFOPyWQiGAzS3Nw84uG5upYW5O9+R9nnn2OVkuVxcey96CJmX3DBuH7Rxu3dy9xbb8XX0sJiITjv/vuZP39+v/t37yMwGo0cccQRapSOMnzBYBCHw0FNTY1Wi7DZbIRCIdxuN0IIkpOTsVgsNDU19XmrPpIZmVLCzp3xvP9+Fh99lEFLiwm73c/9//0l6dv/D8PKlRTu3Mkctxs90AFsSEqibvZs9N/9LrZFixBj/Ac8UEedyWQiOTmZ5ORkkpKSJmTH61iTUlJdXd1jtSeTydRnUA83daWkpODxeNi/f/+I+oZ0TU2IBx5g/qpVxEvJexYL5d/7HrOWLBlWioKRsJeXM/u222jt6GCxTsf3f/lLjj322AHf0/2zlp+fz+TJk1UuHWXofD4f1dXVVFVV4ff7sVqtWK1WXC4XPp8Pk8lEXFycNp28OyEEKSkpWCwWqqqqRlSO8vJ4rrmmDKMxxBEzt3Nc8HlOqX6aBS31xNGZI36T1creqVMJHn88SYsXo59AqQsSEhK0R2Ji4pgHkYNJMBhk165dVFd3ZrQ0Go3ExcVpFQ34z4zhcPNgfHy8Njx0JHQtLfDww8z75BOSQiFW6fWsOe448q+9luRhjs7qPrt5wOZMKcl55x2KHnmE6kCA79vtXHHPPRQXF/f7lu7NOAAFBQVMmjRp2B36KuAfonw+HxUVFdoIifB08YaGBoLBIBaLBSlln+Ot7XY7hYWF7N27d1jT6aWEb75JYvnyHFJTvVx37ipcb77Fv/+RxdLm55kRagBgu15PeV4erqOOIvHMM4kboDNrrHQfYZOamkphYSF2u/2gGVEzkXg8HrZs2aKN6klMTMRut1NdXd0jcIaTlEVz5q6uvZ3QU08x/b33yPF62Q18NHkycskSph57bNT/P00OB1PuuYfMb77hH8BdBQXcfN99ZGZmRnyM/Px8powwfceoB3whxGLgIUAPPCWlvLef/Q4DVgPnSylfHey4KuAPTyAQoKKigsrKSoLBIBkZGRiNRmprawkGg/3eYgshSEtLw+12DztnSlubgfffz2L5O9lUVNqw69u4SPc4j/vvAKAB+Co5mfqSEsynn076ggXjGkjDt9EGg4G4uDicTicAWVlZFBYWqmaaKGlubmbr1q34fD6EEOTk5GCxWKisrMTr9fb4krVareh0umEtKt8XEQwi33yT1P/7P+Y6HASAjywW9hx+OPaLLiJnxowRHd/Y1ET2yy+T/cYbSJ+P24Wg8dxzWbJ06ZCSn0Uj2MMoB3whhB7YDnwHqATWABdKKbf0sd8HgAf4iwr40SelpLa2lt27d+P3+0lPTycuLk5buDvc/t5Xu3Q0siXqW1v53Z2JvL31eA7n31zP45zNq2w2BNgxZQr+E08k+7TTsA5j9uBIGQwGgsGgFtzDtUmDwYDZbKa9vR0hBNnZ2RQUFKhAPwqklOzfv19b0Fuv1zNlyhT0ej2VlZU4nU4tgV54klf4c2mxWKKy2Lxu505Cf/4zxd98Q4bfjwf4LC6OqunTESeeSNZ3vhNRwjad203iunVY3n+fSZ99hiEU4lXgrzNm8L077mDSpEn9vre/BdotFgvFxcUkJCSM4ApHP+AfBdwtpTy56/mdAFLKe3rtdzPgBw4DlquAH11tbW1s374dl8ulJdyqqqqio6NDqz1Fe9GQYBDWLoeX/prGTwO3c27L+1RQxB4SaEncR92CBdjOPpv8uXPHvBYfnvEbzmEuhCAhIQGfz4fb7e4xftpoNJKTk0NOTo5qlx8Dfr+fbdu20djYCIDZbGbatGkYjUaqqqpwOBxIKbFYLASDQS34GwwG0tPTaWpqGnnTTyhE6N//Rrz6KlO3bSO/6463FdhjNuNITqYjOxtht6O32TCYzZibmrA2NpLY0MDk+npMUtIBvCAEnx5+OPPPO4/58+cP+Fnvq7I1a9YszGYzW7duxev1UlRURGFh4cRswxdCnAMsllJe0fX8UuAIKeX13fbJBf4GnAg8zQABXwhxFXAVQEFBwcJ9+/aNqHyxLhgMsmfPHiorKzGZTBQUFNDc3ExjY6P24RrplPDexNZ9LH+knVe3ncze0BQK2MdPuRZrzjrajjuO/O9/n9T09KidL1LhfO/hlLbhUUdms5nm5mY8Hk+PFLl2u53s7GwyMzNjepz8RNXa2sqWLVu04G2z2ZgyZQpxcXHU1dVRXV2N1+vVctB3D/Jmszmq7f1y/35c776L+ZtvSK6rI8flIrdX5SgE1AL7gfV2O/tmzCB49NEcefzxJCcnD/mc4dW5Fi5ciBACv9/Pjh07qK+vJzExkblz5068YZlCiHOBk3sF/MOllDd02+cV4AEp5WohxLOoGn5UtLS0sG3bNjweD9nZ2VitVvbu3Tsqk6Is+/Zheu010lZ+yknOL9nBdIpZzeHpf2Xm92H26YuJj4+P+nkHYzQatcRd4TznqampJCUl4Xa7qamp6dFUZTAYyMjIIDs7G7vdPublVXqSUlJVVcWuXbu0SklCQgKFhYUkJyfT3NxMbW0tDQ0N2nYhBGazOSpNPAPxtrfT0diIu7kZb3s7Ij0dW1ISCQkJI74TDN9tl5aWkpSUpL0e/n20trZSXFwc9WGZ0cg8VAl0zyebB1T32qcMeKmr8GnAqUKIgJTyzSic/5ATCoXYu3cv+/fvx2KxMHPmTKqqqqipqYnqeUxNTcS99Ra+d77in81ncDvL+RxJWeqfOPGEyZy+ZCrx8edE9ZyDMRgMWm4an8+H3+9Hr9eTk5NDcnIyXq+XmpoaHA6H9h6dTqetapScnHzITYiayIQQ5OXlkZ6ezq5du7Q8PRs3biQ+Pp6CggKKi4sJBALU19dTW1uL0+nE4/FoK0+FO9qjzWyzYbbZSC44cM3a4QgPlsjLy6OyshIAt9vdI+AD1NTU4Ha7CQaDUU8OF42jrQGmCSEmAVXABcBF3XeQUmo9GN1q+G9G4dyHHLfbzZYtW3A6nWRlZWGz2di2bVvUji/8fpI+/RTryy9j2O7mfu7kaZ4kiJ7dZ8Vx+uXzuKqf/N2jQafTaW3xfr9fWw0rMTERm82G0WjE7XbT2NjYY56AwWAgNTWVrKysQ3LW68HGbDZTXFxMbm4u27dvp729nY6ODrZs2YLFYiEvL4+srCxyc3Pp6Ohg8+bNWlK/3qLdhBkN6enpOBwOMjIyaG9v1xa3KS8vJxAIkN+1BsOWLVtob2/HarWOSr/XiAO+lDIghLge+DudwzL/IqXcLIS4pmv74yM9h9KpoaGBrVu3IoRg+vTpVFRUUFtbG5VjWysqSHj5ZXI++AC8Zq7nf3iRHyKF4FvH7+GKK5rIyRl4pmC0dO9cDgd6k8mkZRf0+/1ayt7uDAYDaWlpFBQUjNriIcroSkxMpKysjNraWnbt2qV1uu/cuVNbNzYnJ4eysjL27dvH3r17sVgsGI1GLfhPtGA/efJkqqqqsFgsJCcnU15eztSpU8nJyWHr1q3s2rVL63NyOBwYDAYWLFgwKv1KauLVQUBKyZ49e9i/fz/x8fHEx8dHJ9AHg6SuWoX9+ecp2rGDDgy8T4BPZs7lxYaPOPIoDxdfXEFmZvQXER+qcEbK8PBKr9dLMBhEp9ORkZFBTk6OmhwVYwKBAPv376eyspJQKITVatX6ahITE8nKykKn07F9+3ags9O3ra2NuLg4bbau0WiM6roMQ5WcnIxer6exsZHi4mK2b9+O2WzWOmqllJSXl2t/z0IIjjjiiBENCx7tNnxlFAUCAbZs2UJTUxPx8fG4XK4RLyStb28n4+23SX/pJVLa2thIMsfqH2Cj5SL++NhnfC8/jTMCmzEYRr8yEBcXh9lsxmq1EhcXh8Vi0dLD6nQ69Ho9HR0dNDc343A4tLHaKSkpZGZmxnwmykOZwWBg8uTJ5Obmsm/fPmpqapBSkpCQgNfrpby8XFt4Jbxal9FoxOfzkZ+fT0VFBX6/n/j4ePx+vzaqZyybfPR6PQ0NDUyaNEkbUNG7M7b753e0F8NRAX+CCgaD1NfXa7d7wIgDvamhgcy//Y3s5cuJ8/v5B1busv+Wb3w/weezcNLRdSQkZAP+UQn2Op2OpKQk0tPTSUxM7Led0uv10tTURHNzM01NTdrtbkpKCpMmTSI1NXVYaWOVg5PZbGb69Onk5+ezf/9+amtrkVKSkpKCwWCgpaVFmzkers3X1dUxa9YsbW5KeIhuc3OzFuytVmtU0zWHhSeJGY1GGhoaKCwspK2tjfb2dubNm6c1N3q9XjZt2qQ1RSUlJVFTU0NCQsKA+fJHQgX8CSQYDGorCnUfhjZSlpoaMp97jrwPPkB0zQh8derxfOh4m9ZWO0cd1cAVV2xi8uToTGXvLSMjg9zcXBISEvoN8OF1VFtaWnrcjofzyScnJ4/pcnbKxGO1WpkxYwaFhYVa/1UwGMRut5OVlUUoFKKhoQGPx4PP52Pr1q2YzWaEEAQCAZqbm9Hr9drzcLCP5pj+8JcKdH75FBQUEAwGaWxsZOrUqVpm0IqKCqqrq7WZ36WlpcTFxbFx40bKy8u1z360qTb8CcDtdlNVVUVtbW2PFAgjZa6tJfsvfyH/ww8JSMmzwIrSMznjujOYMmUqf/zjVL71LQclJa2DHapPBoNBG2oWroWHP082m42cnBwtj0+Y3+/H5XLhdDpxOp20tbVpf2x6vZ7ExEQt5bDNZlNt8kq/AoEAtbW1VFVV4Xa7tTkY8fHxVFRUaH9LZrMZt9sdcQXKZrNhMpm0wG0ymbBYLMTFxREKhbS7zu70ej3Jyck0NDRor+Xl5eH1enE4HKSnp5OUlERLS0uPIcMWi4XS0lKtGScQCLBu3To8Hg9HHnnkxJt4NZpiPeC7XC727NmjTTFPT0/HYDCMeDy9qamJ3GeeIXfFCoKhEE8Af19wLrXcz5YteTz//Bekpw990QkhBPHx8VqnWPhWtHsCsszMTDIzM7FYLNryh+3t7dqje03KYrFgt9tJSEggKSkJm82mhk8qQyalxOl0UldXR319vdas0z2lyKRJkzCbzezatUubvTtYpcpqtTJt2jQSExO1WnpdXR0tLS2Dlqn7jO7er4dCIYQQFBQUUFBQcEAflM/nw+PxDDunjgr4E0xHRwd79+6lvr4eg8GgjTDZs2fPiHKB691ucl98kdyXXkLn9/M0sGL+sQSzH+GDD+ag10suuWQf555bickU+Wxcu92O3W7X2tZ7p2uIj4/Xbp3DOfW714CEEMTFxWGz2bRRRvHx8YMu4qwoQxUKhXA6nVp6kb7G6Ucrp1RycjJxcXE9vmR6C+ds8vv92t9Neno6kydPHrXOWRXwJ4hQKMT+/fvZt2+fliLWaDSOfAnBUIis998n7/HHiXc6eQV4fuZMjr/iFu6771Lq6y1897u1XHnlbtLSIqvZhxedDi9IMdjnRKfTYbFYtEd4kZXwyBtVc1fGg9/vZ/369T0GPNhsNi01Q/jvzmq1YjKZ8Hq9eL3eAT/vfS13GZ4BHu4X0Ov1mEwmLTNrePnLpKSkUe+LUsMyJ4C2tjbKy8tpb2/X0hZXV1fj9/tH1E5t37KFwgceIG33bj4H7k1LY/YVy7jlu5155s88s5rS0hZmz24b0nGDweABtSOdTofNZiMhIUEbThl+GI1G1d6uTDhGo5GSkhK++eYbvF4vCQkJWiK9nJwcpk2bRlNTEw6Ho8eqU+F+p+419/BdbV93B+GEfQaDgXnz5o04xfFoUQF/lEkpqaysZNeuXRiNRvLy8qirq8PhcGC327VUAUNlaGuj6LHHyHv/faqBHxqN6C+7igTPLfz+94UUFq5j5kwnF1+8P6LjWa1WQqHQAaMV4uPjyczMJD09XeWIVw5K4aC/bt062tramDFjBk1NTVRUVFBRUUFaWhpFRUVYLBZt2U+/308wGNQqOOFJfQ6Hg6qqKpxOJwaDgaSkJFJSUmhpaaG+vp5Zs2ZN2GAPKuCPqkAgQHl5OQ6Hg+TkZC34x8fHk5KSQl1d3dAPKiUZ//wnhQ89hNnl4n7gk+OOY+Hxv+Lppw+jqiqOk06qIyNj8EyCBoMBu91Oa2trj/HI8fHxZGVlkZaWpoK8EhPMZjOlpaWsW7eOnTt3Mm/ePKZMmaKNjguPrrHb7VozpMlkIhAIaHmampubCQaD6PV6ioqKyM/P12bR1tfXk5ubS2pq6jhf6cBUG/4o8Xg8bNy4UWvCCeenz83NxeFwDGvCh9nhYMq995Lx9dd8Afy/9HQW33EH//rXBbz1Vi65uR385Cc7WLiwedBj9V51x2azkZ2dTXp6uloERIlZXq+XdevW4fV6mTVrFunp6doon8bGRlpaWvB4PD3udHU6HSaTieTkZNLS0npkXPX5fKxZswaTyTRq+W+GSrXhjzG328369evx+XykpKTgcDhITEwkJSVFW95tSKQk6733KHr4YYJeLzcKQeN553Hb5ZdjNpvZscPDRRft47LL9mE2Rzb6IBAIYLVaycnJITMzU42YUQ4JZrOZ+fPns2nTJjZv3qytLpWQkNCjKSactC2c5qMvoVCIrVu3EgwGKS4unhDBfjAq4EdZe3s769ev19bhbGpqIjc3F4/Hw549e4Z8PGNTE9PuuYeMtWtZCdyVl8eZN/+WTe8u5ssv6zn22AYuuKAi4uOZzWYyMzPJy8tTQV45JJlMJkpKSti+fTt79+6lvb2dadOm9fh7CNfq+xMIBNi8eTPNzc1Mnz4d2zis0zwcKuBHUUdHB+vWrdMmIXk8HqZNm6atLTtUqatWMeW//xvhcnGTEDRecAHHTfkpv/lNMS6XgblzI58hm5aWxpQpU7BarUMuh6LEGr1ez8yZM7HZbOzZs4fm5mYmT55Mdnb2oKPN3G43GzduxO12M336dHJycsao1COnAn6U+Hw+NmzYoC2zFwgEmDJlCnv37h1ymgSdz8ekRx8l/623WA/cmJrKaT+9hy0fnMkLL2YyfbqT3/9+PZMmDZ77xmq1UlZWdlDcbirKWArPdk1NTWXHjh1s376d6upqrS+rdw3f7XZrqRwA5s2bN6y1bMeTCvhREAwG2bhxI16vF4PBQCgUoqioiJ07d2rPI2WtqmL6L35B8p49/B744IQTuP2WW/jyy8l8/HE6S5fu4eKL90eUzTIzM5OZM2eq8fGKMgCbzUZJSQl1dXXs37+fHTt2sGPHDhITE7Vka36/n7a2zrksycnJTJs27aBcZEcF/BGSUmpLDppMJi3Y79q1S5u5F6n0Tz5h6j330OHz8QOjkck33cFphd8jPr6NE06oZ/p0J3l5kY3uKSwspKioSAV7RYmAEIKsrCwyMzNpb2+nvr6+R9plIQSTJk3S8kQdrFTAH6F9+/bR2NiIxWLB6/VSWFg45GAvgkGKnnySwpdfZjVwa24up1zzCM88cxLV1VZefHE1SUn+iIP9lClTtDUyFUWJXDhBYHx8/HgXZVSogD8CLS0t2pqaHo+H/Px89u7dO6T82sbWVmbcdRdp69fzKPDW8Scwf/bD/PrXs7DZgtx112aSkiJfom3mzJlkZWUN84oURYllKuAPk9/vZ+vWrdponJycHKqqqjAYDBEHe9vu3cy64w6MDQ1cIQSma67Ds/kX/OlPGRx2WBPLlm0lJSXyYF9cXExGRsZwL0lRlBinUhgOg5SSbdu2aQtpJycna4sa9F4YoT+pq1ZRcu21OBsbWWyzMfuBBzjvvB+Qne3hmmt2ce+9GyIK9uE2+tmzZ6tgryjKgFQNfxhqa2tpbGzUMuq53e5+82EfQEryXnmFyY89xjfAdXn5zDjxTcxmCTi55prdEZcjvIjD7NmzSU9PH/qFKIpySFEBf4h8Pp+W+dLv92Oz2Whvj3At2GCQqQ8/TN7bb/Mq8ND849FbX+T557Noba2iuPjAxRr6Ez6/CvaKokRKBfwh2rVrl9ZsM5Rgr/N4mPFf/0Xm6tXcD3x40s3s2fpbamqs/PjHOzn33MqIyxDuFC4uLlbBXlGUiKmAPwTNzc3U1dWh1+vR6XQRB3tjayvFt99O4o4d3AA0XHAvn755G3FxAR58cB3z5kWeIsFqteJ2u5k1a5Zqs1cUZUhUwI9QKBRix44dWrt5pOkSzHV1zP7JTzDV1HC+wcC8X/6Ss48+AqjinHMqSU2NbMlBIQQ2mw2Xy8XMmTPJzMwcwdUoinIoisooHSHEYiFEuRBipxBiWR/bLxZCbOh6fC6EKInGecdSZWUlHR0dQ8qLE7dvH/N+/GOoreVUSx57SjYxa9a30evh6qt3DxrswyNwdDodiYmJuFwupk2bpsbZK4oyLCMO+EIIPfAn4BSgGLhQCFHca7c9wLeklPOAXwNPjvS8YykQCGgLj0cqvrycudddh6ulhRNsR7ExfgsbNkxjxw57n/v3PrZer0dKidFo1JZQmzJlCrm5uSO6FkVRDl3RqOEfDuyUUu6WUvqAl4Czuu8gpfxcShlehmk1kBeF846Z/fv3EwwGI164JGHjRubefDP17e18K+EyNvtXAiYeeugbjjqq8YD9w4sjhxmNRi2ffmpqKg0NDdqSaoqiKMMVjYCfC3RfgaOy67X+/Ah4r7+NQoirhBBrhRBrw5OZhiIYDLJjxw6ampqG/N6+eDwe9u+PbCFwgKSvvmLOrbeyz+Ph5NQfU976LJMmuXnssa+YNavvYZfdg71er8fv92O320lPT6e2tpa8vDwKCwtHfC2KohzaohHw+2rn6LMqLIQ4gc6Af0d/B5NSPimlLJNSlg13yGFTUxPl5eURz3rtTygUYt26dRHvn/LFF8y+4w7K/X5+OHUqv3z4dC64YD9/+MM60tIG75zV6XQEg0FSUlJIT0+noqKCrKwspkyZorJeKooyYtEI+JVA97aGPKC6905CiHnAU8BZUsoD2zWiRAihZaoczpKCYaFQiPXr1+PxeCLaP2X1aop//nM+DSZxcvLz3PI/D5OTE8fVV++OaJ1ZIQShUIisrCzS0tLYvXs3aWlpzJgxQwV7RVGiIhoBfw0wTQgxSQhhAi4A3u6+gxCiAHgduFRKuT0K5+xX98VGqqqqaGlpGdYxNm/eTGtrZOPjU1atovj//T/eDk7mdNPXNLguoqpqaCNppJQUFRWRmprK9u3bSU5Opri4WAV7RVGiZsQBX0oZAK4H/g5sBV6WUm4WQlwjhLima7dfAqnAo0KIdUKItSM97wDlIRQKafmsN27cGHmeG/4T7BsbI7sJSVm9muJf/ILngkdwvv5LTJYsHnxw/ZDWmwWYMWMGCQkJbNmyhYSEBObMmYNOp3LbKYoSPSLSkSfjoaysTK5dO7TvBp/Px+effw5ARkYG9fX1GAwGZs2aRWpq6oDvDYVCbNq0KeIO3+Q1ayi+807+FFzMreI1snMC/M//bCQ3N7JmIOhsypkzZw4Gg4H169djtVopLS3VErMpiqIMhRDiKyllWV/bYm6mrRACnU5HKBTC7/eTlpZGQ0MDGzduJC0tjUmTJmGz2Q54n9frZf369XR0dER0nqSvv6b4Zz9jczDISzP0HJbUyp137iAxse+7CYPBcEAnsk6nY8GCBQCsW7cOs9nMvHnzVLBXFGVUxFzANxqNpKWlUV9fT3NzMyUlJej1eurq6mhsbKShoQG73U5mZiYpKSn4fD7q6uqoqamJ+ByJGzcy645lPBo4nTdn7+Tu3/0Yi2VLv/v3tQKWwWCgrKyMUCjEN998g06nY968eZjN5mFfu6IoykBiLuD7/X7q6+u15xs3bmTRokWEQiEcDgdJSUm43W527tw5rOPby8uZdtudnB/4M2+xhLvPXYvF4up3/77WtrVYLCxcuJBgMMiGDRsAKCkpwWq1DqtMiqIokYjJgN9dKBTiX//6l/Z8OKN2wmx79lD0k59zpu8lPuZ0LrmknOOO6z/Y6/V6fD5fj5m0iYmJzJs3j2AwyPr16/H7/ZSWlvbZzKQoihJNMRfwjUbjAakK4D+5aobbSW2tqiLnprs41f0GX7KIa6/dwLnnDty5GwwGe5QlMzOTmTNnEggEWL9+PV6vl5KSEuz2vvPrKIqiRFPMBXwpZZ9BPfyaTqfrd5/+mBwOZt10Ex86S1grFnLrLWs5/fTBO3e7B/uioiIKCwu1ZpyOjg7mzZtHYmJixOVQFEUZiZgL+IOlL+4+MSsShtZWJt/0U2hs5Im0DTxxz4dMnZoQ0XvDwX7mzJlkZWURCATYsGEDLpeL2bNnk5ycPKSyKIqijETMBXydTqcNyxwpfUcH5hsfYlHNB6RY7uLmP5xKbm5kwT6spKSE5ORkLdg7nU6Ki4tJS0sbcfkURVGGIuYCfltbW1SCvfD54CdPctb+v9KIiUt+ev6Qgr3BYGDhwoVYrVYCgQAbN26kra1NLTquKMq4ibmAn5AwtBp4n4JBAnf8lQu3P4mLILf//ANOOCF70LeF7yxMJhOHHXYYRqMRv9+vBXu16LiiKOMp5gJ+X8nGwuvQRkRK4u59nlPX/R4vHm684x+cdNLkiM4bCoWwWq2UlZVpQzI3bNhAe3s7c+bMUc04iqKMq5gL+CaTibi4OLxerxbkh7IObf6zzzL5w+c5jGym35jD4sXzInqflBK73U5paSl6vR6Px8OGDRvweDzMmTNn0Dw+iqIooy3mAn4oFIo4H05vDU9+TcaLa3gWmHWzlzPPiizYQ2dTUjiNg8vlYuPGjQQCAebOnatG4yiKMiHEXMAPBAIIITAYDASDwYg7cOte3M41L15OBkdz3NInufisswZ/Uxe73c68efPQ6/U0NzezadMm9Ho98+fP19I0K4qijLeYC/h6vR4p5ZBy4Fe/U8V1T56LlRZmfudpLrrs0ojfGw72BoOB6upqduzYgdVqZd68eVgsluFcgqIoyqiIyYAfHh0TieqVTdzw+8XE0cqcw5dx07KrB1xlym6343R2LkYeHx/PvHnz0Ol0bNu2jdraWm2lKpXiWFGUiSbmllRqa2uLONhbamp47bcSMx3MmHktt/3migFXmZo9e7aW+dJms1FSUoLf72fdunXU1tZSUFCg8tkrijJhxVwNP9I1YI0tLcy46SYeCrhYmn8Yt/5+2YCBevbs2dTV1eHz+bRVqWpra9mzZ4+2apUadqkoykQWcwE/Pj6+z2yZ3Tn2hvjbdXr+1OHjxswEbvjjLQPmok9ISKC9vZ2GhgaMRiMzZszQJlOlpqYyffp0tXCJoigTXswF/EAgoAX7vpYVbKoT3HH1FJp8iVxpm8mSh28ZNGNlWloau3fvRgiBzWZj3bp1GAwGZs6cSWZmZsR3FYqiKOMp5gJ+90lWvYN9W4uOn/8wh3pfKt8xnc7Zf/oxGRkZAx7ParWye/duoHNylcvloqioiLy8PAyGmPv1KYoSw2IuYvXXlOPu0HHX5Rns7sjnB7ozOPYP51NYWDjo8dxuN9CZJ6egoIDc3FzVKasoykEp5kbp9DfL1vTEG7S1GLmE81n4P99l1qxZER0vvILWwoULKSoqUsFeUZSDVswF/N7LBYZCkPHyK3z77T+wjGIKflnCwrKyiI6VnJyM3+8nPz9frTmrKMpBL+aadEwmE+np6TgcDqSEJ2+zkPTNXL4L7Lr1BhafcELEx/L5fJjN5oiafhRFUSa6qNTwhRCLhRDlQoidQohlfWwXQoiHu7ZvEEIsiMZ5++J0OnE4HAC8+mvJ/31zJK208OnVV7P49NMjPo7dbqe9vZ0pU6ag1+tHq7iKoihjZsQBXwihB/4EnAIUAxcKIYp77XYKMK3rcRXw2EjP259wG/tHj3h49OMTOIXnSLhwI2dfcMGQjuN2u0lISFALliiKEjOiUcM/HNgppdwtpfQBLwG9U02eBTwvO60GkoQQgy8hNQzBYJBvXmjjt699h6N5j8yzlnPBlVdE/H6z2YxOpyMQCDB16lQ1xl5RlJgRjYCfC1R0e17Z9dpQ94kKfUsLBc88TRmfMum7T3HZTdcOKWh7vV5CoRCZmZnRWS5RURRlgohGp21f0bT3YPhI9uncUYir6Gz2oaCgYMiFqfX7WRH3MTMOa2fpHT8bUrBPSUmhqakJIQSTJw++rKGiKMrBJBoBvxLI7/Y8D6gexj4ASCmfBJ4EKCsr6z8hTj8KCws583//l/j4+AEzX/ZzbgDy8/NVbhxFUWJONJp01gDThBCThBAm4ALg7V77vA1c1jVa50igVUpZE4VzHyAQCJCYmDjkkTU5OTm0tLQghBjWnYWiKMpEN+IavpQyIIS4Hvg7oAf+IqXcLIS4pmv748AK4FRgJ9ABXD7S8/bH5XIN630mkwkpJVlZWSpHjqIoMSkqkU1KuYLOoN79tce7/SyB66JxrgjKov08WJrksOzsbKqqqgCYMmXKqJVNURRlPMVcaoXube+RBHvonGTl9/ux2WwqV46iKDErJgO+0WgccEGT7jIzM9m/fz8AeXl5o1k0RVGUcRVzAd9gMDBnzhwtrfFgEhIS8Hg8AGqJQkVRYlrM9U4Gg0HKy8sj2jc9PZ2qqiqEECQkJKjmHEVRYlrM1fCBiMfQJyYm0tHRgZSS1NTUUS6VoijK+Iq5gK/X63E6nYPul5iYSHV1NSaTCeicZasoihLLYi7gd3R0HLCWbV9SUlLo6OjAbDZjMpnUAieKosS8mAv4fr9/0H0MBgP19fXExcXR3t5OamqqyoqpKErMi7mAb7FYBt0nKyuL9vZ20tLSCIVCqjlHUZRDQswFfJ1ON2gHrNPpxGw2I6VECEFycvIYlU5RFGX8xFzANxqNA9bY4+PjaW1tJS8vj6amJhISElTuHEVRDgkxF+mklOzZs6ff7TqdDoPBQGpqKrt27WLSpEljWDqlN7/fT2VlpTb5TVGUyFgsFvLy8oY0fyjmAr7X6x1wlE5bWxuFhYW0trYCqPH346yyshK73U5RUZHqOFeUCEkpaWxspLKyckiV1phr0rFYLCQmJva5zWAwoNPpyM3NpampSQ3HnAA8Ho8aJaUoQySEIDU1dch3xjEX8EOhUL8TrwKBAFlZWRiNRpqbm0lOTlaBZgJQ/weKMnTD+buJuYAfzovTn7y8PNrb2wkEAmp0jqIoh5SYC/hSSlpaWvrclpKSQlxcHM3NzQAkJSWNXcGUCUsIwa233qo9/93vfsfdd98dlWPr9XpKS0spKSlhwYIFfP7551E57nAsXbqUV199tc/XJ02aRGlpKaWlpRx99NH9HmPv3r387W9/056vXbuWG2+8MSrle/bZZ6mu7nOp64jfL4Tgn//8p/baG2+8gRBCu+7ly5czf/58SkpKKC4u5oknngDg7rvvJjc3V/sdlJaW9htHevvqq6+YO3cuU6dO5cYbb+xzHY4vv/xSO25JSQlvvPEG0DlEvPs509LSuPnmm7XrSU9P17Y99dRTw/7dhMVcp+1AbVrhfPctLS1YrdaIJmkpsc9sNvP6669z5513Rj1FttVqZd26dQD8/e9/58477+STTz6J6jmi4f777+ecc84ZdL9wwL/ooosAKCsro6ysLCplePbZZ5kzZw45OTnDPsbcuXN58cUX+fa3vw3ASy+9RElJCdA5Iuyqq67iyy+/JC8vD6/Xy969e7X3/uQnP+G2224b8jl//OMf8+STT3LkkUdy6qmn8v7773PKKaf02GfOnDmsXbsWg8FATU0NJSUlnHHGGdjtdu3zAbBw4UK+//3va8/PP/98HnnkkSGXqT8xF/D7G6JkNBpJTk4mFArR0tJCRkbGGJdMGczNN9/c48MfDaWlpfzhD38YcB+DwcBVV13Fgw8+yG9/+9se2/bt28cPf/hDHA4H6enpPPPMMxQUFLB06VISEhJYu3YttbW13HfffYMGzLa2Nq0Z0eVycdZZZ9Hc3Izf7+c3v/kNZ511FgC//vWveeGFF8jPzyctLY2FCxdy2223sWbNGn70ox9hs9k45phjeO+999i0aRPBYJBly5axcuVKvF4v1113HVdffTVSSm644QY++ugjJk2aFPEKcGGffPIJN910E9B5F/Tpp5+ybNkytm7dSmlpKUuWLGH+/Pn87ne/Y/ny5dx9993s2bOHmpoatm/fzu9//3tWr17Ne++9R25uLu+88w5Go5Ff/epXvPPOO7jdbo4++mieeOIJXnvtNdauXcvFF1+M1Wpl1apVbNmyhVtuuQWXy0VaWhrPPvss2dnZA5b52GOP5V//+hd+vx+v18vOnTspLS0FOmvTgUBAG5lnNpuZMWPGkH4nvdXU1NDW1sZRRx0FwGWXXcabb755QMCPi4vTfvZ4PH22v+/YsYP6+nqOPfbYEZVpIDHXpKPX6/v8Zebk5CCEwOVyEQwGVXOO0sN1113HCy+8oA3XDbv++uu57LLL2LBhAxdffHGP5ouamho+++wzli9fzrJly/o8rtvtprS0lJkzZ3LFFVfwi1/8AugcTfbGG2/w9ddf8/HHH3PrrbcipWTt2rW89tprfPPNN7z++uusXbtWO9bll1/O448/zqpVq9Dr9drrTz/9NImJiaxZs4Y1a9bw5z//mT179vDGG29QXl7Oxo0b+fOf/zxgc9Ltt9+uNR1cfPHFQGfT1p/+9CfWrVvHv/71L6xWK/feey/HHnss69at4yc/+ckBx9m1axfvvvsub731FpdccgknnHACGzduxGq18u6772q/0zVr1rBp0ybcbjfLly/nnHPOoaysjBdeeIF169ZhMBi44YYbePXVV/nqq6/44Q9/yM9//vPB/hsRQnDSSSfx97//nbfeeoszzzxT25aSksKZZ55JYWEhF154IS+88AKhUEjb/uCDD2q/gxNOOAGA8vLyHk0uvZt8qqqqeqyUl5eXp62P3dsXX3zB7NmzmTt3Lo8//vgBEz5ffPFFzj///B7x67XXXmPevHmcc845VFRUDHr9g4m5Gr6Uss+aTH5+PoDWLqcC/sQzWE18NCUkJHDZZZfx8MMP91gec9WqVbz++usAXHrppfz0pz/Vtp199tnodDqKi4upq6vr87jdm3RWrVrFZZddxqZNm5BS8rOf/YxPP/0UnU5HVVUVdXV1fPbZZ5x11llaGc444wyg83PrdDq19vWLLrqI5cuXA/CPf/yDDRs2aO3Ura2t7Nixg08//ZQLL7wQvV5PTk4OJ554Yr/X31eTzqJFi7jlllu4+OKL+f73vx/REqCnnHIKRqORuXPnEgwGWbx4MdDZ1BJuPvn444+577776OjooKmpidmzZ2vXGVZeXs6mTZv4zne+A3QubDRY7T7sggsu4OGHH6a1tZUHHniA//7v/9a2PfXUU2zcuJEPP/yQ3/3ud3zwwQc8++yzQN9NOjNmzBjwrrOvWNPf6JkjjjiCzZs3s3XrVpYsWcIpp5zSo1n5pZde4n//93+152eccQYXXnghZrOZxx9/nCVLlvDRRx9F8ivoV8wF/L4mXcXFxWnfps3NzcTFxUW8SIpy6Lj55ptZsGABl19+eb/7dP9j7v4ZiqS55KijjqKhoQGHw8GKFStwOBx89dVXGI1GioqK8Hg8/R5noONLKfnjH//IySef3OP1FStWjGjI67JlyzjttNNYsWIFRx55JB9++OGg7wn/TnQ6HUajUTu/TqcjEAjg8Xi49tprWbt2Lfn5+dx999199rtJKZk9ezarVq0acrkPP/xwNm3ahNVqZfr06Qdsnzt3LnPnzuXSSy9l0qRJWsDvS3l5Oeeff36f21auXEleXh6VlZXaa5WVlYP2QcyaNQubzcamTZu0/o/169cTCARYuHChtl/3SaFXXnkld9xxx4DHjURMNun0lpmZCXSO0W9tbVW1e6VPKSkpnHfeeTz99NPaa0cffTQvvfQSAC+88ALHHHPMsI+/bds2gsEgqamptLa2kpGRgdFo5OOPP2bfvn0AHHPMMbzzzjt4PB5cLpfWDJKcnIzdbmf16tUAWpkATj75ZB577DEtNfj27dtpb2/nuOOO46WXXiIYDFJTU8PHH388pPLu2rWLuXPncscdd1BWVsa2bduw2+0RLTDUn+7rR7tcrh6jhrofe8aMGTgcDi3g+/1+Nm/eDMAjjzwyaEfmPffc06NmD539JitXrtSer1u3jsLCwgGPE67h9/VISkoiOztb+3+RUvL8889rfTHd7dmzR6uM7tu3j/LycoqKirTtL774IhdeeGGP99TU1Gg/v/3228yaNWvAskYi5mr4fdWE0tPTgc5Om1AopMbfK/269dZbewSThx9+mB/+8Ifcf//9WqftUITb8KHzs/ncc8+h1+u5+OKLOeOMMygrK9Pa+AEOO+wwzjzzTEpKSigsLKSsrEybOf70009z5ZVXYrPZOP7447XXr7jiCvbu3cuCBQuQUpKens6bb77J9773PT766CPmzp3L9OnT+da3vtVvOW+//XZ+85vfaM+//PJL/vCHP/Dxxx+j1+spLi7mlFNO0XJRlZSUsHTpUubPnz+k30dSUhJXXnklc+fOpaioiMMOO0zbtnTpUq655hqt0/bVV1/lxhtvpLW1lUAgwM0338zs2bPZtm0bixYtGvA8vTtNofP3f99993H11VdjtVqx2Ww9avcPPvggf/3rX7Xnb775Zo+g3J/HHnuMpUuX4na7OeWUU7Rzv/3226xdu5Zf/epXfPbZZ9x7770YjUZ0Oh2PPvpojxFhL7/8MitWrOhx3Icffpi3334bg8FASkrKgHcikRJD7bkfS2VlZbJ7p1UkfD5fj84pvV7PMcccgxCCvXv3snfvXhYtWqQWLJ8gtm7dGpWaSyxxuVzEx8fT0dHBcccdx5NPPsmCBQu01wHuvfdeampqeOihh8a5tGPv9NNP5/XXX9eWJz2U9fX3I4T4SkrZ51jZmK/hJyUlae2ILS0txMfHq2CvTGhXXXUVW7ZswePxsGTJEhYsWADAu+++yz333EMgEKCwsDAqNb6DUbizWhm6EQV8IUQK8H9AEbAXOE9K2dxrn3zgeSALCAFPSilHrVpiNBrR6/UEg0EArfkmFArR1tYWcU+/ooyX7jNZuzv//PP77UBUlEiMtNN2GfBPKeU04J9dz3sLALdKKWcBRwLXCSGKR3jefkkptWAPaO2cLpeLUCjUbyZNRVGUWDfSgH8W8FzXz88BZ/feQUpZI6X8uutnJ7AVyB3heQclhECv12ttnuEJNSrgK4pyqBppwM+UUtZAZ2AHBsxXIIQoAuYDXwywz1VCiLVCiLUOh2PIBQqFQqSmpmIwGEhISNDa71tbW7FYLGr8vaIoh6xBA74Q4kMhxKY+HgcONh34OPHAa8DNUsq2/vaTUj4ppSyTUpaFh1MOhdFoZObMmfj9fm28vZSStra2AdMmK4qixLpBA76U8iQp5Zw+Hm8BdUKIbICuf+v7OoYQwkhnsH9BSvl6NC+gL72bbzweDz6fTzXnKH1SKYwjp1IY90xhvH//fk444QTmz5/PvHnzDhhLP+GEc88M5wHcDyzr+nkZcF8f+wg6R+n8YajHX7hwoRyOnTt3ypUrV8pAICCllLKmpkZ+/PHH0ul0Dut4yujZsmXLeBdB2mw27ef3339fHnfcceNWliVLlshXXnkl4tf78vHHH8vTTjst2kWTUkr5rW99S65Zs2bY73/mmWfk3Llz5Y9+9CPttfPOO0+WlJTIV155Rfp8PpmdnS0rKiqklFJ6PB65bds2KaWUd911l7z//vuHdd7DDjtMfv755zIUCsnFixfLFStWHLBPe3u79Pv9Ukopq6urZXp6uva8uwULFshPPvlESinllVdeKR999FEppZSbN2+WhYWFwyrfcPX19wOslf3E1JGOw78XeFkI8SNgP3AugBAiB3hKSnkqsAi4FNgohFjX9b6fSSlH7auwtbUVu92upVlobW1Fr9er9WsnuptvhiinR6a0FIaQlE2lMFYpjCHyFMZCCNraOluoW1tbR5TLfyyMqNNWStkopfy2lHJa179NXa9XdwV7pJSfSSmFlHKelLK06zFqwT68pm33fDnh9nu1dqrSF5XCWKUwDhtqCuO7776bv/71r+Tl5XHqqafyxz/+cdDrH08xN9NWp9Nx9NFHazUlv99Pe3s7w+kAVsbYOKVHVimMVQrjsKGmMH7xxRdZunQpt956K6tWreLSSy9l06ZN6HQTMy9lzAV86LnqVfh2S3XYKpFQKYxVCmOIPIXx008/zfvvvw90fnY8Hg8NDQ0TdkW9ifk1FEXhETtqSKYSCZXCWKUwhshTGBcUFGgjjrZu3YrH45nQrQkxWcPvrq2tjfj4+D7z5CsKqBTGvakUxpGnMH7ggQe48sorefDBBxFCaMNOJ6qYS4/cXSgU4rPPPiM7O5tp06ZFsWRKtBys6ZFVCuOBqRTGY+OQT4/cXUdHB6FQCLvdPt5FUWKMSmE8MJXCeGKK6YAfbmtU7fdKtKkUxsrBKKY7bdva2tDr9drwOUVRlENZTAd8p9OJ3W6f0J0oiqIoYyVmA34wGKS9vV015yiKonSJ2YDvcrmQUqoOW0VRlC4xG/BVh60ymMbGRi3nSlZWVo/Uuz6fb8D3trS08Oijj2rPV65cyemnnz7aRVaUEYnZUTpOpxOTyaRWuFL6lZqaquVgufvuu4mPj++RpyUQCByQPCssHPCvvfbasSiqokRFzAb8trY21ZxzEDr++ANfO+88uPZa6OiAU089cPvSpZ2PhgbolVuMbrP0I7J06VJSUlL45ptvWLBgAXa7vccXwZw5c1i+fDnLli1j165dlJaW8p3vfIfTTjsNl8vFOeecw6ZNm1i4cCF//etf1YABZUKJyYDv9/txu91kZWWNd1GUg9D27dv58MMP0ev13H333X3uc++997Jp0ybtDmHlypV88803bN68mZycHBYtWsS///1vjjnmmLEruKIMIiYDvsvlAlA1/IPQQDXyuLiBt6elDb1G35dzzz13WLmXDj/8cC09cWlpKXv37lUBX5lQYrLTNpwSWQV8ZTi6r4xmMBh6LMDRV5rgsO79RXq9Xsu8qCgTRUwGfKfTidVq7ZEXX1GGo6ioiK+//hqAr7/+mj179gCMOAWxooyHmAz4qsNWiZYf/OAHNDU1UVpaymOPPaYt1pGamsqiRYuYM2cOt99++ziXUlEiE3PpkUOhENu3byc5OZnMzMxRKpkSLQdremRFmQgO+fTIOp1OW6hCURRF+Y+YbNJRFEVRDqQCvjLuJnKzoqJMVMP5u1EBXxlXFouFxsZGFfQVZQiklDQ2NmKxWIb0vphrw1cOLnl5eVRWVuJwOMa7KIpyULFYLNpEv0ipgK+MK6PRyKRJk8a7GIpySFBNOoqiKIcIFfAVRVEOESrgK4qiHCIm9ExbIYQD2DfMt6cBDVEszsHgULxmODSv+1C8Zjg0r3uo11wopUzva8OEDvgjIYRY29/04lh1KF4zHJrXfSheMxya1x3Na1ZNOoqiKIcIFfAVRVEOEbEc8J8c7wKMg0PxmuHQvO5D8Zrh0LzuqF1zzLbhK4qiKD3Fcg1fURRF6UYFfEVRlEPEQR3whRCLhRDlQoidQohlfWwXQoiHu7ZvEEIsGI9yRlsE131x1/VuEEJ8LoQoGY9yRtNg19xtv8OEEEEhxDljWb7REsl1CyGOF0KsE0JsFkJ8MtZljLYIPt+JQoh3hBDru6758vEoZzQJIf4ihKgXQmzqZ3t0YpmU8qB8AHpgFzAZMAHrgeJe+5wKvAcI4Ejgi/Eu9xhd99FActfPpxzs1x3JNXfb7yNgBXDOeJd7jP6vk4AtQEHX84zxLvcYXPPPgP/p+jkdaAJM4132EV73ccACYFM/26MSyw7mGv7hwE4p5W4ppQ94CTir1z5nAc/LTquBJCFE9lgXNMoGvW4p5edSyuaup6uBoeVQnXgi+b8GuAF4Dagfy8KNokiu+yLgdSnlfgAp5cF+7ZFcswTsQggBxNMZ8ANjW8zoklJ+Sud19CcqsexgDvi5QEW355Vdrw11n4PNUK/pR3TWDA5mg16zECIX+B7w+BiWa7RF8n89HUgWQqwUQnwlhLhszEo3OiK55keAWUA1sBG4SUoZGpvijZuoxLKDOR++6OO13mNMI9nnYBPxNQkhTqAz4B8zqiUafZFc8x+AO6SUwc6KX0yI5LoNwELg24AVWCWEWC2l3D7ahRslkVzzycA64ERgCvCBEOJfUsq2US7beIpKLDuYA34lkN/teR6d3/hD3edgE9E1CSHmAU8Bp0gpG8eobKMlkmsuA17qCvZpwKlCiICU8s0xKeHoiPQz3iClbAfahRCfAiXAwRrwI7nmy4F7ZWfj9k4hxB5gJvDl2BRxXEQllh3MTTprgGlCiElCCBNwAfB2r33eBi7r6uE+EmiVUtaMdUGjbNDrFkIUAK8Dlx7ENb3uBr1mKeUkKWWRlLIIeBW49iAP9hDZZ/wt4FghhEEIEQccAWwd43JGUyTXvJ/OOxqEEJnADGD3mJZy7EUllh20NXwpZUAIcT3wdzp79v8ipdwshLima/vjdI7WOBXYCXTQWTM4qEV43b8EUoFHu2q8AXkQZxiM8JpjTiTXLaXcKoR4H9gAhICnpJR9Du07GET4f/1r4FkhxEY6mzrukFIe1CmThRAvAscDaUKISuAuwAjRjWUqtYKiKMoh4mBu0lEURVGGQAV8RVGUQ4QK+IqiKIcIFfAVRVEOESrgK4qiHCJUwFcURTlEqICvKIpyiPj/ZuBydDrD06MAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#Now let's do bagged prediction using linear regression\n",
"from sklearn import linear_model\n",
"\n",
"n_boots = 100\n",
"\n",
"boot_est = dict()\n",
"x_grid = np.arange(0, 1, 0.01)\n",
"d_grid = pd.DataFrame(x_grid, columns=['x'])\n",
"d_grid = bd.makePolyFeat(d_grid, 6)\n",
"\n",
"\n",
"fig = plt.figure()\n",
"\n",
"#Generate and plot each bootstrap model\n",
"for i in range(n_boots):\n",
" D_b = dat.iloc[np.random.randint(0, dat.shape[0], size=dat.shape[0])]\n",
" regr = linear_model.LinearRegression(fit_intercept=True)\n",
" regr.fit(D_b.drop('y', 1), D_b['y'])\n",
" boot_est[i] = regr.predict(d_grid)\n",
" plt.plot(x_grid, boot_est[i], color='0.75')\n",
"\n",
"#Now aggregate the bootstrapped models for a single prediction\n",
"bag_est = pd.DataFrame(boot_est).mean(axis=1)\n",
"\n",
"#For comparison, we'll also see what a single fit on the original data looks like\n",
"regr = linear_model.LinearRegression(fit_intercept=True)\n",
"regr.fit(dat.drop('y', 1), dat['y'])\n",
"non_bag = regr.predict(d_grid)\n",
"\n",
"#Truth\n",
"truth = d_grid[['x','x2','x3']].dot(np.array(betas[1:]))\n",
"\n",
"#Now Get MSE of estimates\n",
"mse_bag = round(np.sqrt(((truth-bag_est)**2).sum()), 3)\n",
"mse_non = round(np.sqrt(((truth-non_bag)**2).sum()), 3)\n",
"\n",
"\n",
"plt.plot(x_grid, non_bag,'k-', label='Non Bagged Estimate, MSE={}'.format(mse_non))\n",
"plt.plot(x_grid, bag_est,'r-', label='Bagged Estimate, MSE={}'.format(mse_bag))\n",
"plt.plot(d_grid['x'], d_grid[['x','x2','x3']].dot(np.array(betas[1:])), 'b--', label='Truth')\n",
"\n",
"\n",
"plt.title('Bagging Example')\n",
"plt.legend(loc=4)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
We can see in the above example that variance of the bootstrapped predictions is very high, especially around the extreme values of $X$, where there is less data. The goal of Bagging is to reduce the variance of the prediction (and thus improve accuracy). The lever that controls the amount of variance reduction is the number of bootstrap samples used.\n",
"
\n",
"When to use Bagging
\n",
"According to the original Bagging paper by Leo Breiman, bagging will produce a much better reduction in training error when the underlying classifier or regression model is unstable, and sensitive to minor variations of the data $D$. The above example was meant to show Bagging at work. We used a linear regression on polynomial features, and even though the underlying model was misspecified (we used degree 6 where the truth is degree 3), the bagged estimate isn't too far from the estimate without bagging. One must be careful when choosing to use bagging, because it can actually hurt performance if the underlying model is fairly stable (if you rerun the above regression example with different data sets or different bootstrapping iterations, you might find that the bagged MSE is worse).\n",
"
\n",
"\n",
"## Random Forests\n",
"### The Basic Algorithm\n",
"
The Random Forest algorithm is probably the most well known and utilized implementation of the Bagging technique. A RF is an ensemble of Decision Trees, where both bagging and random feature selection are used to reduce the variance of the forest. The basic algorithm goes as follows:
\n",
"Assume we have a data matrix $D=[X,Y]$ with $N$ records and $M$ features.
\n",
"\n",
"Train \n",
"For each $b$ of $B$ iterations:\n",
"
\n",
"
Draw a bootstrap sample $D^b$ of size $N$ from $D$.
\n",
"
Sample $p$ features from $X$, where $p<\n",
"
Grow a Decision Tree $T_b(X)$ on this data
\n",
"
\n",
" \n",
"Score \n",
"Take the average of all of the tree predictions, i.e. \n",
"$RF(x)=\\frac{1}{B}\\sum\\limits_{b=1}^B \\: T_b(x)$\n",
"\n",
"\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
We'll start by building a forest using SKlearn.ensemble. Check out here for a brief overview of Ensemble methods in Python. We'll also compare this to a single tree. This test is close to an off-the-shelf test, in the sense that we are not doing any intelligent hyper-parameter optimization. However, we do use a few well reasoned starting parameters."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from sklearn.ensemble import RandomForestClassifier\n",
"from sklearn.tree import DecisionTreeClassifier\n",
"import pandas as pd\n",
"import course_utils as bd\n",
"imp.reload(bd)\n",
"\n",
"\n",
"# f = 'C:/Users/kevin/Documents/GitHub/DS_course/datasets/Cell2Cell_data.csv'\n",
"f = 'data/Cell2Cell_data.csv'\n",
"dat=pd.read_csv(f, header=0, sep=',')\n",
"train, test = bd.trainTest(dat, 0.8)\n",
"lab = 'churndep'\n",
"\n",
"\n",
"#We'll build a RF and compare to a DT\n",
"clf_def = DecisionTreeClassifier(criterion='entropy', min_samples_leaf = 20)\n",
"clf_def = clf_def.fit(train.drop(lab, 1), train[lab])\n",
"dt_pred = clf_def.predict_proba(test.drop(lab,1))\n",
"\n",
"rf_def = RandomForestClassifier(criterion='entropy', n_estimators=100)\n",
"rf_def = rf_def.fit(train.drop(lab, 1), train[lab])\n",
"rf_pred = rf_def.predict_proba(test.drop(lab,1))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
We can see that with little optimization involved, the RF is already a bit better. Let's see how much better we can do with a grid search on both the DT and the RF. We have already covered DT optimization in a different module, so here we'll focus on the input parameters of the RandomForestClassifier object that let us tune the forest.
\n",
"\n",
"The following two are forest specific:\n",
"
\n",
"
n_estimators - the number of trees (and bootstrapped samples) to be used
\n",
"
max_features - the number of features that will be randomly sampled for each tree.
\n",
"
\n",
"The default in RandomForestClassifier is max_features=sqrt(total_features), which is generally a good suggestion. The default for n_estimators is 10, which is probably too low. The other design parameters are specific to the individual decision trees, which are covered in the decision tree module.\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Out-of-Bag Error\n",
"
Usually we would want to use cross-validation for performance tuning. For Random Forests, this becomes problematic, as cross-validation with RF's can be painfully slow. That's because each cross-validation step requires building k*n_estimators trees. One advantage of a random forest is that it performs what is called an \"out-of-bag\" error calculation. Remember that each tree is built from a bootstrap sample of the data, so that for each tree, some portion of the data is not used for that tree. The RF method computes an out-of-bag prediction for each record $[x_i, y_i]$ by averaging the prediction $f^b(x_i,y_i)$ on record $i$ for the bootstrap iterations in which record $i$ was not chosen in the bootstrap. The out-of-bag prediction can then be used to compute out-of-sample error for model selection and validation. This method should be equivalent to $N$-fold cross-validation. \n",
"\n",
"
\n",
"We'll start by building a tree of oob predictions enabled and then compare the performance to a true hold out set.\n",
"\n",
"
The training AUC is 1.00. This means that the RF perfectly predicts the training set!
\n",
"
The OOB AUC is exactly equal to the test AUC up to 2 digits of precision.
\n",
"
\n",
"Now we'll use OOB error to do model selection."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"n_est = [50, 100, 200, 500, 1000]\n",
"m_feat = [1, 3, 6, 11]\n",
"\n",
"aucs_oob = {}\n",
"aucs_test = {}\n",
"\n",
"for m in m_feat:\n",
" aucs_oob[m] = []\n",
" aucs_test[m] = []\n",
" for n in n_est:\n",
" rf_oob = RandomForestClassifier(criterion='entropy', n_estimators=n, max_features=m, oob_score=True)\n",
" rf_oob = rf_oob.fit(train.drop(lab, 1), train[lab])\n",
" aucs_oob[m].append(roc_auc_score(train[lab], rf_oob.oob_decision_function_[:,1]))\n",
" aucs_test[m].append(roc_auc_score(test[lab], rf_oob.predict_proba(test.drop(lab,1))[:,1]))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABaJUlEQVR4nO3dd3RU1drH8e+TThJIowVC6B0EJAKCBUQFK4L32lBRuoq9YXvVa7l2r71AANtVrkqwF0BRUUGKlCSAhJ6ElkJ6n+f945yEIQZSyGRS9metWTNz6p6ZZH5z9j5nb1FVDMMwDKOqPNxdAMMwDKNhMcFhGIZhVIsJDsMwDKNaTHAYhmEY1WKCwzAMw6gWExyGYRhGtZjgMOodEdklIme7uxxNlYg8IiLvu7scFRGRiSLyvbvL0dSZ4GggROQ6EdkkIrkisl9E3hCR4HLL9BGRz0UkQ0SyRORHERnuNL+TiKiIZNu3AyLyuoh4V7JvEZEdIhJfwby/fcnbZV3h9NzH/jLaJiI59jrzRKRTTd+PqhCRkfbrXVRu+gB7+nIX7PMRESlyeo+zReSeWthmvfgid3pPXys3fYWIXHec9dT+7Kv8vjj9vXqVTlPVD1T13BN+IRXvb7mITHXFthsbExwNgIjcCTwN3A0EAcOAjsASEfGxl+kK/ApsAjoD7YAY4HsRObXcJoNVNRDoD5wK3FRJEc4AWgNdROSUGryET4CLgavs8g8A1gKja7Ct6joEDBeRMKdpk4C/XLjPhaoa6HR7xoX7cocc4NoaBP+Axvq+2D+umsz3aZN5oQ2ViLQAHgVuVtVvVbVIVXcBl2GFx9X2oo8Av6vqA6qapqpZqvoy8B5W6PyNqh4ElgB9KinGJOAz4Gv7cXXKfzZwDjBOVVerarGqZqjqa6oafZxVTxGReBFJF5H5IuJnby9WRC5y2r63iKSIyMBjbKcQWAxcYS/vifXefVCunC+JyF4RyRSRtSJyutO8r0XkeafnC0VkXnXeB3u9ySKy2X5N34lIx8r2LyJjgfuBy+1f6RuOse3ZIrLdPtKMF5HxTvOus48InrP3vVNEznOa31lEfrLXXQK0rOSlHAYWAA9X9z04RtmHiMga+7UfEJEX7Fk/l+7Pfu2nVnA0qyJyo300myUij4lIVxH53d7e/5x+XIWIyJcicsh+H74UkQh73hPA6cCr9r5etacPF5HVYh3Fr5ajj+CXi8gTIvIrkIv1w+o6sY7Os+z3eWJtvEf1jqqaWz2+AWOBYsCrgnnvAB/aj/cD11ewzCigBPAHOgFaui2so5INwOTj7N8fyATOBy4FUgAfp/m7gLPLrXMdsMJ+/BTwUzVf8y4gFugAhGIdST1uz7sH6xd96bLjgE3H2M5IIBEYDqyyp50PfAdMBZY7LXs1EAZ4AXfa76efPa8tcBA4C5gI7ACaH2OfjwDvVzD9EiAB6G3v40Hgtyruv8Jtltv+P+3P0wO4HOuoINzp8ygCpgGewA1AMiD2/N+BFwBfrKPLrGPtz+k9bWv/XfS0p68ArjtO+RTodox5vwPX2I8DgWH24044/b2W/9ty2u7nQAugL1AALAO6YB3dxgOT7GXDsP6G/YHmwMfAYqdtLQemOj0PBdKBa+zP5Ur7eZjT8nvs/XrZ+3N+T8KBvu78/nDVzRxx1H8tgRRVLa5g3j6O/DpsaT+vaBkPIMRpWoqIHAaSsL5gPjnO/idg/TN+D3yJ9Q9yQTXKH3aMclXmVVXdq6ppwBNY/7QA7wPn20diYP1Tv3e8Danqb0CoiPQErgXerWCZ91U1Va0jouexvkR72vP2AzOxgvol4FpVzTrOLi8TkcNOt3bADODfqrrZ/iyfBAaWHnUcb/9Voaofq2qyqjpUdSGwDRjitMhuVZ2jqiX26wgH2ohIJHAK8JCqFqjqz8AXVdjffuBN4F9VLSOwrtz7MsaeXgR0E5GWqpqtqiursU2Ap1U1U1XjsH5wfK+qO1Q1A/gGGGSXOVVVP1XVXPvzewI48zjbvQDYpqrv2Z/Lh8AW4CKnZRaoapz9mRYDDqCfiDRT1X12mRodExz1XwrQ0rmB0Em4Pb90ufBjLOPA+qVUqqWqBmP98voV+PY4+58E/M/+xykAFnF0dVUxUL5x3RvrywAg9Rjlqsxep8e7sX5No6rJdpkvFevkgPMoV+10DO8Bs7COwGLKzxSRO+1qpAw7VIM4usrmS6xf61tVdUX59cv5n6oGO92SsaoVXyr90gTSAAHaV3H/xyUi14rIeqft9yu3/v7SB6qaaz8MxHpf01U1x2nZ3VXc7dPAGBEZUK4scXKkAfx0p1knl3tfvrOnTwF6AFvs6qALq7j/UgecHudV8DzQLpe/iLwlIrtFJBOrKixYrOrLirTj7+/FbuzPzFb2d2q/h5dj/cjYJyJfiUivar6WBsEER/33O9Yv/gnOE0UkAOtLc5k9aSlWdUV5l2G1feSWn6GqeVh11aeKyN++pOz637OAq8U6k2s/8A+sX/yly+/BqlJw1pkj/3BLgSGldcnV0MHpcSRW1Uqpd7Cqdv6J9dqSqrC994Abga/Lvxf2l9u9WO9ViB2qGVhf7KWeADYD4SJyJdW3F5hR7ouzmar+VoX9H7cLa/uoZQ5WMIbZ68eWK/+x7ANC7L+nUpFVeUGqmgr8B3is3PS+eqQB/JcqbGebql6JdQLG08Andnlqu+vuO7GO4oaqagusajk49vtcGvjOIrGO1EsdtY6qfqeq52D9WNqC9bk0OiY46jn7cPtR4BURGStWY3AnrPrZRI5U0zyKdfbQEyISKiLNReRmrKqZeyvatoj4YlX17Mc6MijvGqyzj3oCA+1bD3u/pV+eC4HbRKSXWKKAycBHdvmXYjXAx4jIYBHxsss2U0QmH+el3yQiESISitU4vNBp3mLgZOBWKqh2qoiq7sSqlniggtnNsY6cDgFeIvJ/WHXmAIjIGcD1WO/ltVifRfsKtnM8bwL3iUhfe5tBIlIa9MfdP9Yv6E5y7LN2Sr9kD9nbvh7riKNSqrobWAM8KtZp06dxdFVMZV7AakPqXY11jiIiV4tIK1V1YDW8g9UudwjraLlLTbddTnOsI5DD9t9V+cb9A+X29TXQQ0Susv9uL8c6keTLY7yONiJysR16BUC2/ToaHRMcDYBapy3eDzyH1fi2CusX7Gi7+ghV3QachnWq6y6sX5KXAmNU9ddymzwsItlY/yinAherakW/7iYBr6vqfucb1pdgaXXVHGA+Vr14BtYX+QOq6lz99Q+sf8KF9jKxQBTW0cix/BerXWWHfXvc6f3IAz7FOrJZVOHaFVDVFXa1UXnfYdWF/4V1pJSPXQVht6W8C8xS1SS7mioamC8iVflFX7rvGKxf0x/Z1SSxWEeMx92/7WP7PlVE1lWw7Xjgeayj0wNYp1mX/8yP5ypgKFb12cNUMYztfWcCz2A1JFdmgxx9Hcd/7OljgTj7b/Il4ApVzbePDJ8AfrWr4IZV/SVV6D9AM6xq3ZX8vYr2JeAfYp1x9bJ9RHUh1pFKKtaJGReqagoV87CXTcZ6L8/EOsptdKTi7wvDqN/sX+U9VPXqShc2DKNWVdTgahj1ml3NMAWrKs0wjDpmqqqMBkVEpmFV43xjnzpqGEYdM1VVhmEYRrWYIw7DMAyjWppEG0fLli21U6dO7i6GYRhGg7J27doUVW1VfnqTCI5OnTqxZs0adxfDMAyjQRGRCnsRMFVVhmEYRrWY4DAMwzCqxQSHYRiGUS0mOAzDMIxqMcFhGIZhVIsJDsMwDKNaTHAYhmEY1dIkruMwDMNoTBzqIK84j5yiHHKKcsgtyi17nFN89POLu15MZIsqjc1VZSY4DMMwXExVKSgpOPIlX5xTpS9952WPml78twE9KyQIA1oNMMFhGIZRF4pKisq+zP/2BW9/eZf/Mv9bGDitW6JVGwzQz9MPf29/ArwDCPAOwN/LnzC/MCKbR1rPS+d5BRy1XNk8ryPPm3k1oxrjjVWZCQ7DMBqFEkdJ2Zd3RV/c5b/0K/qSd16uyFFUpf16eXj97Ys80CeQNgFt8PfyP/pLvfRL/hhf+v5e/nh51P+v5fpfQsMwjArkFeexZPcSYrbFEJcaR15xXpXWE+Rvv9wDvAMI8Qs55i/38kcAzl/4Pp4+Ln6l9Y8JDsMwGgxVJTYllkUJi/hm5zfkFOUQ2TySS7tfSnOf5n/70v/br3ovf5dV3zQlJjgMw6j30vLT+HL7l8QkxJBwOIFmXs04p+M5jO82nsFtBpsgqGMmOAzDqJdKHCX8mvwrixMW8+PeHyl2FHNSy5N4+NSHGdtpLIE+ge4uYpNlgsMwjHplb+ZeYhJi+Gz7ZxzMPUiIbwhX9bqK8d3G0y2km7uLZ2CCwzCMeiCvOI+lu5cSkxDD6v2r8RAPRrQbwX1D7uPMiDPx9vR2dxENJyY4DMNwC1UlLjWORdushu7somw6NO/ALYNu4eKuF9MmoI27i9hwqULKX7BtCQy4AgJa1urmTXAYhlGn0vPT+XKH1dC9LX0bfp5+VkN39/FEtYkyDd01lZ8JO3+GhCWQsAwy9lrTQzpC74tqdVcuDQ4RGQu8BHgCc1X1qQqWGQn8B/AGUlT1zOOtKyKPANOAQ/Ym7lfVr135OgzDODEljhJ+3/c7i7YtKmvo7t+yPw8Ne4jzOp9Hc5/m7i5iw6MKB+KOBMWe38FRDD7NocuZcMZd0HU0BHeo9V27LDhExBN4DTgHSARWi8jnqhrvtEww8DowVlX3iEjrKq77oqo+56qyG4ZRO/Zm7WVxwmI+S/iMA7kHCPYN5oqeVzC++3h6hPRwd/EanrzDsONHSFhqhUXWPmt6m/5w6izofg50GAoubhNy5RHHECBBVXcAiMhHwDgg3mmZq4BFqroHQFUPVmNdwzDqofzifJbuWUrMthj+2P8HHuLB8HbDueeUexjVYZRp6K4OhwP2b7CCYttSSFwNWgJ+QdBllBUUXUdDi/A6LZYrg6M9sNfpeSIwtNwyPQBvEVkONAdeUtV3q7DuLBG5FlgD3Kmq6bVcdsMwqkFViU+NJyYhhq93fE1WURYRgRHcPOhmLu56MW0D2rq7iA1Hbhps/8Fq2N6+DHLsWvnwgXD6HdDtbGgfBZ7ua6J25Z4rauHSCvY/GBgNNAN+F5GVlaz7BvCY/fwx4Hlg8t92LjIdmA4QGVm7XQobhmE5nH+Yr3Z+xaJti/gr/S98PX05p+M5TOg+gcFtBuMhZqy4SjlKIPlP+6hiCSStBRSahUK30VZQdB0Nga3cXdIyrgyORMC5VSYCSK5gmRRVzQFyRORnYMDx1lXVA6UTRWQO8GVFO1fVt4G3AaKiosoHlmEYNVTiKGHlvpXEJMTww54fKHIU0TesLw8Ne4ixncfSwqeFu4tY/2UfdDqq+AHy0gCBiCgYORu6nQPtBoKHp7tLWiFXBsdqoLuIdAaSgCuw2jScfQa8KiJegA9WddSLwJZjrSsi4apqtwgxHoh14WswDMOWmJVoNXRv/4z9OfsJ8g3i8p6Xc0m3S+gZ2tPdxavfSoohaY0VFAlLYd96a3pAK+gxxj6qOAv8Q91azKpyWXCoarGIzAK+wzqldp6qxonITHv+m6q6WUS+BTYCDqzTbmMBKlrX3vQzIjIQq6pqFzDDVa/BMJq6/OJ8lu1ZRsy2GFbtX4UgDG83nLui7mJUh1FNskvxKsvcZ5/9tNQ6Eyo/A8QTOgyBsx6ywqLtSeDR8KrzRLXx1+JERUXpmjVr3F0Mw2gQVJXNaZtZtG0RX+/8mqzCLNoHtueSbpdwSbdLTEP3sRQXwt5VR8LigF0Z0jzcColuZ0OXkdAs2J2lrBYRWauqUeWnmyvHDcMAjjR0x2yLYWv6Vnw8fDi749lM6D6BU9qeYhq6K3J4r9NRxU9QmAUe3hA5DM5+1DpdtnUfaGRXw5vgMIwmzKEOViZbDd3L9iyjyFFEn7A+PDD0Ac7rfB5BvkHuLmL9UlwAu387EhaHtljTgzpA/39YQdH5DPBt3FfCm+AwjCYoKTup7IrufTn7CPIN4rKelzG+23jT0F1e2s4jQbHzZyjKBU8f6DgCTr7WqoJq2aPRHVUcjwkOw2giCkoKWLZ7GTEJMazatwqAU9udyh1RdzCqwyh8PX3dXMJ6oigPdq04cl1F2nZrekhnGHS1FRSdTgOfAPeW041McBhGI7c51Wro/mrnV2QVZtEuoB03DLyBcV3H0S6wnbuL536qkJpwJCh2/wrF+eDVDDqfDkNnWGER1tXdJa03THAYRiOUUZDBVzu+IiYhhi1pW/Dx8GF0x9GM7zaeoeFDTUN3QTbs+uXIdRWHd1vTw7pD1GQrKDoOB+9m7i1nPWWCwzAaCYc6WLVvFTHbrIbuQkchvUN7c//Q+zm/8/lNu6FbFQ5uPtJWsed3KCkE7wCrC/IRt1hhEdLJ3SVtEExwGEYDl5ydzGcJn7E4YTHJOcm08GnBpT0uZXy38fQO6+3u4rlPfoZ1imzpeBWZSdb01n3s6qdzrNNmvUzbTnWZ4DCMBqigpIAf9/zIom2LWLlvJQDDwodx2+DbOCvyrKbZ0K0K+zcdCYq9q6yBjXxbWBfenXmv1WlgUIS7S9rgmeAwjAZkS9oWYrbF8OWOL8kszCQ8IJyZA2Yyrts42ge2d3fx6l5JsVX1tPlz6z7b7gO17Ukw/BbruoqIU1w+sFFTY4LDMOq5jIIMvt75NTHbYticthlvD2/OjjybS7pfwrDwYU2zoTttJ/z5Pqz/wBoFzy/Y6iSw29nWUUVz0y2KK5ngMIx6yKEO/tj/B4u2LWLZbquhu2dIT2YPmc2FXS5smg3dRfmw5UtY9451IZ54WO0U5z9n9TBrjipQVRw5uZSkp1GSnk5xWhrNTjoJr9Da7XXXBIdh1CP7svexeLt1RXdSdhLNfZozofsExncfT5+wPu4unnsciIN178LGhZCXDsGRMOpBGHgVBDXu6jktKaHk8OGyEChJP0xJetqRx2lp1vOyx+loYeFR2+jw9lsEnnFGrZbLBIdhuFlhSSE/7P2BmG0x/J78O4oyNHwotwy6hbMiz8LPy8/dRax7BVkQ+6kVGElrrS4+el1odfHR+cwG2RU5gCMvzw6BdErS0/8eAofteXYIlGRkWI3+FfAIDMQzNBSvkBC827TBr3dvPEOC8QoNxTMktOyxT9fav3DRBIdhuMnWtK3EJFgN3RkFGbQNaMuMATMY13UcEc2b4Jk/qpC42qqKio2Bohxo1RvG/BtOuhwCwtxdwqOow4EjM9P6oj9sfdkXp6dTUvrFXy4EitPT0by8ijfm6YlnSAheIcF4hoTi27On9cUfEopnaKhTIITgGRKKV0gw4uO+sVBMcBhGHcoszOSbHd+wKGER8anxeHt4c1bkWUzoNoGh4UPxrKdDhbpUTips/Mg6uji0xboor/+lcPIkaD+4zjoP1MJCq8onPe3oEEh3flx6lJBOyeHDUFJS4bbE3x+v4GDrSz8sFN9uXe2jgBA8Q0OcQsB67NG8OdKAjqJMcBiGiznUwer9q4lJiGHp7qUUlBTQI6QHs4fM5oLOFxDsF+zuItY9hwN2LrfCYvOX4CiyTpu9+BXoO/6EuyW3GolzrF/7ziFw2G4rsEOgOP3IY0d2dsUbE8EzKMj+5R+CT6dONBt0sv2lH2JND7YDIcR67uHXuKsXTXAYhovsz9nP4oTFLE5YbDV0ezfnkm6XWA3doX2QJtQNd5mMRFj/X1j3HmTsgWYhMGQaDLoG2lS98V9VKU5OJi8+nvz4eIp277FCwLmRuKiownXFx6csBLxCQvDpEHkkBEqrgkofh4biGRSEeDbBI8HjMMFhGLWosKSQH/f+SMy2GH5L/g1FGdJ2CLMGzeLsyLObZkN3SRFs/cY6uti+DNQBXUbBOY9YDd6VdPmhqhTt3Ut+XBz58fHkx1lhUXL4sLWApyfe7dvjFRqKd3g4fn37WFVBwdYX/1EhEByCR4B/0wztWmSCwzBqQXZhNh//9THvxb/HobxDtPFvw7STpnFJt0vo0LyDu4vnHinbrLDY8CHkHILm7eD0u2DQxGN2JqglJRTu3l0WDvlxceRv3owjK8tawNsb3+7daH7O2fj16YNfnz749uzZ6KuG6hsTHIZxAlLzUvlg8wd8tPUjsgqzGBo+lEeHP8rwdsObZkN3YS7Ef2YFxp7fwMMLeoy1Grq7jQan90SLiynYseNISMTHk795M5qbC1hVSr69etHigvPx69vXConu3fFw49lEhsUEh2HUQFJ2EgtiFxCTEENhSSFndzybyf0m069lP3cXzT2S11thseljKMiE0K5w9qMw4Epo3gYtLKRgy1by4+PJs6ucCrZsRQsKAJBmzfDr3ZvgCROsI4m+ffHt0hnxNleD10cmOAyjGralb2Ne7Dy+2fkNIsJFXS7iun7X0SWoi7uLVvfyDltBse5d2L8RvPygzyU4+l1JQX6odQTx3WtWSPz1V1ljtUdgIH69exNy5ZX49bWqm3w6dTIN0A2ICQ7DqIL1B9cTvSma5YnLaebVjKt6X8W1fa6lbUAT60xP1Rpadd17EL8YR34B+R69yPedSH6GP/nvJFCQMKvs+gaPoCCa9e1D6KRry44kvDt0aFDXLBh/Z4LDMI5BVVmRtILo2GjWHlhLkG8QNw64kSt7Xdn0rr3IOkDJygXkL/0v+XsOkZ8RQH5OJIWHcsCRAfyIZ2gofn37EjhqpN1w3Rfv9u3MGUyNkAkOwyin2FHMkt1LiN4Uzdb0rbTxb8M9p9zDpd0vxd/b393FqxMlhw+THxdL/i9fkL/6F/L3HKIwq/TrIgiv1q3w69uPFn36WNVNffvi1bq1CYkmwgSHYdgKSgr4LOEz5sfOJzE7kc5BnXlsxGNc0PkCvBtxl93FaWnWaa+lZzfFbqQoeX/ZfO9Axa9LR4KGjsLvlNPx690br1at3Fhiw91McBhNXnZhNgu3LuS9+PdIzU+lX1g/7oq6i1GRoxrdIElFBw6SH3/0hXTF+51CIsQbv8BMggcU49evP37nTsIr6lIz1oVxFBMcRpOVkpfCB5s/YOGWhWQVZTEsfBhP93+aIW2HNPgqF1WleN++slNfS28lh1KsBUTw6dwZ//498BvaEr+Ctfj5p+HZqgOcfKM91kUT7KHXqBITHEaTk5iVyIK4BSxOWFx2DcaUflPo27Kvu4tWI2VdcpQeRdhh4dwlh2/XrgSOOM1qtO7RGb+SzXjEfwRJH1ljXURdYI91MbLBjnVh1B0THEaT8Vf6X8yLnce3O79FRLi468Vc1/c6Ogd1dnfRqkwdDgp37T7Sb5N9q7RLDl9fSFxjjXWx7F57rIte9XasC6N+M8FhNHp/HvyTuZvm8nPizzTzasbE3hO5ts+1tAlo4+6iHVdZlxylAREXT8HmzTgq6pKj9Grr8l1y5KTCn/PssS42W2Nd9JtgdQESEVVnY10YjYsJDqNRUlV+SfqF6E3RrDu4jmDfYG4ceCNX9bqKIN8gdxfvb7SwkILt24/uAXbrVjQ/H7C75OjVi6CyLjn64NulS8VdcjgcsPMnKyy2fAklhdA+Ci562QqNExzrwjBMcBiNSrGjmO93fU90bDR/pf9F24C2zB4ym/HdxtebazAcBQUU/PXXUe0RR3XJERCAX58+hFx+edk1ElXqkiMjyRrr4s934bA91kXUFDj5GmjTMNtvjPrJpcEhImOBlwBPYK6qPlXBMiOB/wDeQIqqnnm8dUUkFFgIdAJ2AZeparorX4dR/xWUFLB422Lmx80nKTuJLkFdeHzE45zf+Xy3XoPhyM0l3+7cr7Sb8ILt26G4GLC65PDr0/tIlxx9+uAdGVn1LjlKiuCvb60uQBKWWGNddD4TRj9sjXXhbbobN2qfy4JDRDyB14BzgERgtYh8rqrxTssEA68DY1V1j4i0rsK6s4FlqvqUiMy2n9/rqtdh1G9ZhVks3LqQ9+PfJzU/lf4t+3P3KXczqoN7rsEo3LuX7B9+KDsNtnDHTqvqCI50yTFypN253wl0yZGSYB1ZrP8Qcg5C83A4/U4YOBFCG05jv9EwufKIYwiQoKo7AETkI2AcEO+0zFXAIlXdA6CqB6uw7jhgpL3cO8ByTHA0OSl5Kbwf/z4Lty4kuyib4e2GM6XfFE5pe4rbrsHI/fNP9k6fgSMrC6/WrfHr25cWY8aW9QDr1abNiZWtMBc2f261Xez+FcQTep5nnUbbdTR4mppno2648i+tPbDX6XkiMLTcMj0AbxFZDjQHXlLVdytZt42q7gNQ1X2lRynlich0YDpAZGTkib0So97Ym7WXd+LeIWZbDEWOIs7peA5T+k+hT1jVx6t2hZzff2fvTbPwbtWKDh//D59OnWpv4/s2WGGx8WMoyIDQLnD2IzDgKmhev88MMxonVwZHRT+ttIL9DwZGA82A30VkZRXXPS5VfRt4GyAqKqpa6xr1z9a0rdY1GLu+xUM8GNd1HNf3u56OLTq6u2hk/fAjSbfdhk+nTkRGz62dfpzyDkPsJ1Zg7Ntgj3Uxzjq66DjCnEZruJUrgyMRcB5sOQJIrmCZFFXNAXJE5GdgQCXrHhCRcPtoIxw4iNForTuwjrmb5vJL0i/4e/lzTe9ruKbPNfXmGoyMr74i+d7Z+PXpQ+Tbb+EZHFzzjanCnt+tsIhbDMV50KY/nP8c9P+HdZaUYdQDrgyO1UB3EekMJAFXYLVpOPsMeFVEvAAfrOqoF4Etx1n3c2AS8JR9/5kLX4PhBqrKz4k/Ex0bzZ8H/yTEN4RZA2dxRa8r6tU1GOkff8z+/3sY/6goIt54A8/AgJptKPugfRrte5CaAL4tYOCV1tFF+EBzdGHUOy4LDlUtFpFZwHdYp9TOU9U4EZlpz39TVTeLyLfARsCBddptLEBF69qbfgr4n4hMAfYA/3TVazDqVrGjmG93fcu82HlsS99GeEA4s4fMZkL3CTTzaubu4h0ldcECDj71NAFnnE7ESy/h0aya5XOUQMIyqwuQv74FRzFEDofT77KqpHzqxzUnhlERUW381f9RUVG6Zs0adxfDOIb84nwWJyxmQdwCkrKT6BrUlcn9J3Ne5/Pw9qhf3XmrKimvv07KK6/S/Nxzaf/cs4hzFx+VSd8Nf74P6z+AzCTwb2n1RDvoGmjVw3UFN4waEJG1qhpVfro5f89wm8zCTP639X+8F/8eaflpnNTqJO495V7O7HBmvRwHQ1U5+OxzpM2bR9AllxD++GOIVxX+hYoLYMtXVtvFjuXWtG5nw9inoMdY8KpG8BhGPWCCw6hzKXkpvBf/Hv/b+j+yi7IZ0W4EU/pPIapNVL0dB0MdDvY/+i8OL1xIyFVX0ebBByq/uvvgZuuK7g0fQl4aBEXCyPtg0EQz1oXRoJngMOrM3sy9zI+bz2cJn1GsxZzb8Vwm95tM77De7i7acWlxMcn33U/mF18QNm0are64/dgBV5ANcTHW0UXiH+DhDb0vNGNdGI2KCQ7D5bambSV6UzTf7f4OT/FkXLdxXN/3eiJb1P8LMx2FhSTdcQfZS5fR6vbbaTljesULFhfAmnnw87OQmwote8KYJ+2xLlrWbaENw8VMcBguoaqsPbCW6NhoViStwN/Ln0l9JnF1n6tp7V/hxf71jiM3l8RZN5Pz22+0eeABQq+5uoKFSmDjQvjxScjYC11GwpmzIXKYOY3WaLRMcBi1yqEO6xqMTdGsP7SeUL9Qbh50M5f3vLxeXYNRmZKsLPbOvIG8P/8k/IknCL50wtELqMLWr2HZv+DQFmg3CC5+BbqOck+BDaMOmeAwakWRo4hvd1rXYCQcTqBdQDvuG3If47uPr3fXYFSmOD2dvVOnkb91K+1feJ4WY8cevcCuX2HpI1YbRlg3+Oc71rUX5gjDqAcKikvYdiCbuOQM4pIzmXJaZzqG1fDi1GMwwWGckLziPBYnLOaduHdIyk6iW3A3njztScZ2HlvvrsGoiqKDB9kzeTJFexPp8NqrBJ555pGZ+zZaRxgJS6B5O2tEvYETTa+0httk5RexeV9WWUjEJWey7UAWxQ7r+rwAH09G925jgsOoHzILM/loy0d8sPkD0vLTGNBqALOHzOaMiDPq5TUYVVGYmMSeyZMpSUmhw9tvEzB0iDUjdbvVhhH7CfgFwzmPwZBp4N2wjqSMhu1QVkFZQMQnZxKXnMGu1Nyy+S0DfejTLoiRPVvRt10L+rYLomOoPx4etX8kbILDqLZ1B9Yxa9kssoqyOK39aUzpN4XBbQbX22swqqJgx072TJ6MIzeXyPnzaDZgAGTth5+esboF8fSxBkoafgs0C3Z3cY1GTFVJTM8rC4nYJOv+YFZB2TIdQpvRNzyIS0+OoG97KyRaN/ets/9BExxGtaw7sI6ZS2fSNqAt886YR6/QXu4u0gnL37KFPVOmAtDxvXfxi2xjVUmtfANKCmHwdXDG3dC8rXsLajQ6xSUOth/KcapqyiA+OZPMfGtoYU8PoVurQE7r1pI+9lFEn3YtCGrm3mpgExxGlTmHRvS50bTyr4VxJ9wsb/169kyfgUdAAJFvv47vwW8h5gXIPwz9/wmj7rcGTjKME5RXWMKW/ZllbRHxyRls2Z9FQbE1tLCftwe92rbgogHt6NsuiL7tWtCzbXP8vD3dXPK/M8FhVEljDI2clavYe+ONeIWF0fHOC/D+fDxkJUO3c2D0/0H4Se4uotFAZeQWHXUUEZecyfZD2dht1rTw86JvuyCuGdaxrKqpS8sAvDwbRvugCQ6jUusOrOOGpTfQxr9NowmNrOXLSbr1NnxaB9NhdAreK/8PIobApXOg02nuLp7RQKgqBzKtRuvYpCMhkXQ4r2yZti386NuuBef1a0sf+0giIqRZg24TNMFhHFdpaLT2b828MfMaRWhkfvMNSXfdhV9LDzpErcereS+45EPoeZ65FsM4JodD2ZWaU1bVVNoekZpTCFh/Op3DAhgUGczVwzraZza1ICzQ180lr33HDQ4R6Qa0UdVfy00/HUhW1e2uLJzhXn8e/LPRhcbhuS+w7/k5NGtZQIfz/fEc8wacdBl41L96ZMN9Cosd/HUgq+y017jkTDbvyySnsAQAb0+hR5vmjO7duqw9old4CwJ9m8Zv8cpe5X+A+yuYnmfPu6iWy2PUE38e/JOZS2Y2ntBI2Ubak7M48PUeAto7iPi/2/AYMQO8avfXYH5RCXvTcvHz9sTXywPf0nsvjwZdNdGYZRcUs3lfJnFJThfRHcyiqOTIRXS9w1vwj8ERVki0b0H31s3x8WoY7RGuUFlwdFLVjeUnquoaEenkmiIZ7taoQiMjCV3+b1I/+IxDGwNpPjCSdnP+i0fzsFrfVdLhPK58eyV70nIrnO/r5VEWKBXd+3l74Ovlia9973ci9yawKpSSXXBUg3V8cia7UnMoHQg1LMCHPu1acEaPLmVVTZ3CAlxyEV1DVllw+B1nnrlsthFqNKGRmwYrXkBXvs2h9c1IjQ8k6IIxhD/9XNVG7aumvWm5XDlnJRl5RTw1oT+eHkJBsYP8ohIKih0UFJWQX3pf5KCg+Oj73MJi0nKcp1vLFhQ7KCxxnFDZKgss53vf4wRaVe+t4HNvYB25iM467bX0SGJ/Zn7ZMhEhzejbrgXjB7Uvu9K6TYu6u4iuIavsP2i1iExT1TnOE0VkCrDWdcUy3GH9wfUNPzQKc2Dl6/Dry2h+Fgd2DyY9PpngKy6n7f/9X+Wj9tXAnlQrNLILivnv1GH0j6jdXoBLHEqhUwjlF5WQX1xCQdHR0yq6L6hkvqsDqyyEvD3w8/KsVmAd897e1pFtenA4zz79Ncm+RmJfJhl5RQB4CHRrHcipXcPo266FdSFdeBBB/g2vL7X6orLguA2IEZGJHAmKKMAHGO/Cchl1bP3B9cxYMoPW/q2JHtMAT7ktLrS6BvnpGcg5iHY/n31/BJOxcjmhUybT+q67XPJLcldKDlfNWUluUQkfTB1Kv/a133W8p4fQzMeTZj5124BfUWBVdl+dwErPdZ5uhVdB0YkFlq+XB73aNuf8/uH0s6+P6FVPL6JryI4bHKp6ABguIqOAfvbkr1T1B5eXzKgz5UOjoQy0BIDDAbGfwo+PQ/ou6DgCx4QFJP9nIVlLltDq1lsImznTJaGx41A2V81ZRUFxCf+dOow+7VrU+j7cyV2B5XBolYOq9D7A15M+4UF0bdVwLqJryCo7HTfUfrjBvilw2MVlMupQgw0NVdi2xOpT6sAmaNMfJn6Co/0IEm+9jZxffqHNfbMJnTTJJbtPOJjNVXNWUuJQPpw+jF5tG1douJOHmwLLqLrKqqrWYoWFON0HisgGYKqq7nJt8QxXWn9wPTOXzqSVf6uGFRp7VsGyR2H3rxDSCS6Nhr4TKMnNJXH6DHLXriX88ccI/sc/XLL7bQeyuHLOKsAKjR5tmrtkP4ZRX1VWVdW5oukiMgF4Exhb0Xyj/isNjZbNWjJvzLyGERoH4uGHx6whWwPbwAXPw6BrwcvHGrVv+gzyN2+m/fPP0eL8811ShK37s5g4dyUiwofThtGttQkNo+mp0XmJqrpIRB6s7cIYdaPBhUb6blj+b9jwEfg2h7MegmE3gI81qlnxoUPsmTyFwt27iXj5ZZqf5Zpxvzfvy2Ti3FV4ewr/nTaMrq0CXbIfw6jvahQcIhIImBaoBsg5NKLPrefVU9mH4JfnYU00iAcMvxlOux38Q8sWKUpKYvfkyRQfSqHDW28ScOqpLilKXHIGV89dha+XJx9OH0bnlrU7FKdhNCSVNY7fUcHkEOBi4FWXlMhwmfKh0SagjbuLVLH8TPj9Nfj9VSjKhUFXw5mzIaj9UYsV7trF7usn48jOJjJ6Lv6DBrmkOLFJGUycu4oAHys0anv8ZsNoaCo74ihfgavAfuBqVd3kmiIZrtAgQqO4AFZHwy/PQW4q9BkHox6EVj3+tmj+1r/YM2UKOBx0fPcd/Hr3dkmRNuw9zDXRq2ju581H04fRIdTfJfsxjIakssbxRyuaLiJ+IvJPVf3YNcUyatOGQxuYuXQmYX5h9TM0HCVW+8Xyf0PGXugy0hpIqf3gChfP27iRPdOm4+HnR+S77+DbxTUj9P25J51r5/1BsL83H04bRkSICQ3DgGq0cYiIJ3AucCUwBvgFMMFRz204tIEZS2YQ5hfGvDHz6ldoqFpnSC37FxzaAu0GwcWvQNdjN27n/PEHiTNvwDM0lMgF8/GJiHBJ0dbuTmPSvNWEBfrw32nDaB9sumYzjFKVBoeInAFcBVwA/AGMADqrasVdgBr1Rr0OjV0rYOkjkLgawrrDZe9C74uPO5BS9s8/k3jzLXhHRBA5LxrvNq55Pat3pXHdvD9o3cKP/04bSniQCQ3DcFZZ43gisAd4A7hbVbNEZKcJjfpvw6ENzFxiV0+NqUfVU/s2WEcYCUuheTu46GUYOBE8j/8bJvPb70i6+258u3cjcu5cvEJDj7t8Ta3akcr1C1bTNsiPD6cNo02L43UQbRhNU2Wn1H4KtAcuBy4SkQCsBvIqEZGxIrJVRBJEZHYF80eKSIaIrLdv/+c071YRiRWROBG5zWn6IyKS5LSOa670asBKQyPUL5ToMdG0DWjr7iJB6nb4ZDK8dQYkroFzHoNb1sHgSZWGxuGYxSTdcQfN+vWj44IFLguN37ancN381bQLbsZH001oGMaxVNY4fqv9pT0Kq23jWSBIRC4DvlbV7GOta7eJvAacAyRiddH+uarGl1v0F1W9sNy6/YBpwBCgEPhWRL5S1W32Ii+q6nNVfZFNSb0Ljaz9Vo+1694BTx84/U4Yfgs0C67S6mkffMCBxx7H/9RhdHjtNTz8XdNAvWJbClPfXU1kqD8fTB1Gq+aNb5xow6gtlbZxqKoCPwA/iIg3VjcjVwKvAy2Ps+oQIEFVdwCIyEfAOKB8cFSkN7CytEpMRH7C6sb9mSqs22RtPLSRmUtmEuIX4v7QyDsMv74EK98ARxEMvg7OuBuaV71MKW/P4dALLxA4ejTtX3geD1/XfJn//Nchpr27hs4tA/hg6lDCAk1oGMbxVOvKcVUtEpFwVb1KRCprMWwP7HV6nggMrWC5U+1OE5OBu1Q1DogFnhCRMKzxzc8H1jitM0tErrWn3amq6eU3KiLTgekAkZGRVXuBDdjGQxuZsWQGIX4hzBszz32hUZQHq96CFS9C/mHo/08YdT+EVv2UWVXl0Iv/IfXtt2lx4YW0+/eTiLdrBt35cetBZry3lm6tAnl/6lBCA3xcsh/DaExq0m3ITABVzatkuYpOjynfPrIO6KiqA4BXgMX2tjcDTwNLgG+xunQvttd5A+gKDAT2Ac9XtHNVfVtVo1Q1qlWrBjYoUTXVi9AoKYa1C+DlQbD0YYg4BWb8ApfOrV5oOBwceOJJUt9+m+DLLqPd00+5LDSWbT7AjHfX0qNNIP+dZkLDMKqqJn1VVXVEnESgg9PzCKyjijKqmun0+GsReV1EWqpqiqpGA9EAIvKkvb3SwaWwp88BvqzBa2g03B4aDgds/gx+eBxSEyBiiBUWnU6r9qa0pIR9Dz5ERkwModdfT+t77nbZ+M/fx+3npv+uo3d4C96bPNQMI2oY1VCT4LioisutBrqLSGcgCbgC63qQMiLSFjigqioiQ7COgFLtea1V9aCIRAITgFPt6eGqus/exHisaq0madOhTe4Nje0/wNJHYd96aNUbrvgQep533GsxjkULC0m6516yvv2WlrNm0fKmG10WGt/G7mfWf9fRr30Q70weQlAzExqGUR1VuQDwTCBdVTfaZ1OdISLbgddVteBY66lqsYjMAr4DPIF5qhonIqVVXW8C/wBuEJFirLaMK+zGeIBP7TaOIuAmp3aMZ0RkIFa11y5gRrVfdSOw6dAmpi+Z7p7QSNsJX94GO5ZDUCRc8iacdBl41GzENkd+Pom33krOTz/T+t57Cbv+utos7VG+2riPWz76kwERQSyYPIQWfiY0DKO65Mj3dAUzRV4DTgJ8gb+AQKw2h+GAp6pOrItCnqioqChds2ZN5Qs2EKWhEewbzPyx8+s2NDb+D768w+rmfNR9EDUZvGp+FlJJdg6JN95I7urVtH3kEUIuv6wWC3u0LzYkc9vC9ZwcGcz864cQ6FujUQUMo8kQkbWqGlV+emX/OaNUtY+I+GFVN7VW1RIReQvY6IqCGsdXWj1V56GRnwlf3wUbF0LkqTDhbQg+sbPVSg4fZs/0GeTHxdHumWcIuujCyleqoc/WJ3H7wvVEdQpl/nWnEGBCwzBqrLL/nnwAVc0Xkd2qWmI/VxEpcnnpjKOUhkaQb1DdhsbeP+DTqZCRCKMegNPuqPRq78oUp6SwZ8pUCnfsIOLll2g+enQtFfbvPl2byN2fbGBo5zCir4vC38eEhmGciMr+g1rbgzmJ02Ps5437HNd6JjYltu5Dw1ECv7xgdXce1B6u/wYiK7oUp3qK9u1jz/WTKTpwgIg33yBwxIhaKGzF/rdmL/d+upERXVsy59oomvnUrB3GMIwjKguOORwZzMn5McBcl5TI+JvYlFimfz+dIN+gumsIP7wXFk2HPb9Bv3/AhS+AX9AJb7Zw9252X389jswsa9S+k0+uhcJW7KM/9nBfzCZO62aFhp+3CQ3DqA01GsjJqDvlQyM8MLwOdrrIOmvKUQLj34KTLq/RKbbl5f9lj9pXVEzkOwto1rfviZf1GD5YtZsHYmI5s0cr3rpmsAkNw6hFlV45LiLnicjPIpIiIodE5CfTI23dKA2NFr4t6iY0CrLhs5vgk+utMTJm/gIDrqiV0MjbFMuea65FEDq+/55LQ+Pd33fxQEwsZ/VqzdvXmtAwjNpW2Xgc07Cuk7iHI31FRQFPiUiEqr7t4vI1Wc6hMX/MfNeHRtI6qwE8bQecfheMnA2etXONQ+6aNeydMRPP4GBr1L4OHSpfqYbm/7qTR7+I5+zebXht4iB8vUxoGEZtq6yN43bgNFVNc5r2g4icB6wATHC4QFxKXN2FhsMBv70MPzwGgW3gui9r1F3IsWT/soLEm2/GOzycyPnz8G7ruvaZub/s4PGvNjOmbxteufJkfLxq0hWbYRiVqSw4pFxoAKCqqa7qDqKpi0uJY9r30+omNDKTIWYm7PwJ+oyDC/8D/rU3SFLm99+TdOdd+HbtSmT0XLzCwmpt2+W99dN2/v3NFs7v35aXrhiEt6cJDcNwlcqCI1NEBqjqBueJIjIAyHJdsZom59BweZvGlq/gs1lQnA8XvwKDrqmVtoxSGZ99RvL9D9CsXz86vP0WnkEnfkbWsby+PIFnvt3KRQPa8eJlA/AyoWEYLlVZcNwJfC4i84G1WP1DnQJMAq52cdmalLiUOKYtORIa7QLbuWZHhbnw/QOwZh6ED4BLo6Fl91rdRfqHH7L/0X/hP2wYHV57FY+AgFrdvrNXlm3j+SV/MW5gO57/pwkNw6gLlZ2Ou8LutfZOrLAQYCcwTFX310H5moSy0PBxcWjs3wSfTIGUrdbwrWc9BF61OwZFanQ0B599jsCRI2n/0n9cNmofwH+W/sV/lm5jwqD2PPvPAXh6mOpTw6gLlZ1V5YUVGpOBPVjBcSbgEJEHVNV0O3KC4lLrIDQcDlj1pjXAUrNQuCYGup5Vq7tQVQ69/DKpb7xJi/PPo93TT7tsACZV5cUlf/HyDwn8c3AET116kgkNw6hDlVVVPYt1tXhnVc0CEJEWwHP27VbXFq9xi0u12zRcGRrZB2HxDZCwFHqcB+NehYDjDRVffarKgX//m/R33yPo0gmE/+tfiKdrToNVVZ79biuvL9/OFad04Mnx/fEwoWEYdaqy4LgQ6OE0RgaqmikiNwBbMMFRY3USGn99D5/dCAVZcMHzEDWlVhvAwR617+GHyfjkU0KuvYY2s2cjHq5pZ1BVnvpmC2/9vIOrhkby+Lh+JjQMww0qCw7VCgbssLtWP/ZAHsZxOYdG9Jjo2g+NonxY+gisegNa94VJX0Dr3rW7D0CLiki+914yv/6GsBtm0uqWW1w2ap+q8sRXm5m7YifXntqRRy/u67J9GYZxfJUFR7yIXKuq7zpPFJGrsY44jGqKT423Lu6zQ6N9YPva3cHBzdYV4AdiYehMOPtR8Par3X0AjoICkm69jezly2l9152ETZ1a6/sopao8+kU8C37bxXXDO/HwRX1MaBiGG1UWHDcBi0RkMkefjtsMa7xvoxriU+OZ9v00Ar0Daz80VGFNNHz3APgEwlUfQ49za2/7Thx5eey94UZyV66k7cP/R8iVV7pkP2CFxsOfx/Hu77uZclpnHrygtwkNw3Czyk7HTQKGishZQF+ss6q+UdVldVG4xsQ5NOaNnVe7oZGTCp/Pgq1fQ7ezYdzr0LxN7W3fiaOwkMRZN5P7xx+EP/Vvgi+5xCX7AXA4lIc+i+WDVXuYcUYXZp/Xy4SGYdQDVRoKTVV/AH5wcVkaLZeGxvYfrW5D8tJgzL+t6ilXNU4XF5N8513k/Por4U887vLQuD9mEx+t3suNI7ty95ieJjQMo54wY2i6mMtCo7jQ6pjwt5ehZU+4+hNo2792tl0BdTjY9+BDZC1ZQpv7ZhN86aUu21eJQ7n30418sjaRm8/qxh3n9DChYRj1iAkOF3JZaKQkwKeTYd8GiJoM5z4BPv61s+0KqCoHnvw3GYsX03LWLEInTXLZvkocyt0fb2DRn0ncdnZ3bju7h8v2ZRhGzZjgcBGXhIYq/PkefHMvePnC5R9A7wtPfLuVOPTyy6S//z6hkybR8qYbXbaf4hIHd368gc/WJ3PnOT24eXTt9qFlGEbtMMHhAptTNzPt+2kEeAfUXmjkpcMXt0L8Z9D5DGtI1xYu6tPKSWp0NKlvvEnQPy6l9ex7XVZlVFzi4LaF6/ly4z7uHtOTm0Z1c8l+DMM4cSY4atnm1M1M/X6qFRpjaik0dq2ARdMh+4B1XcbwW1zWAO4s/aOFHHz2OZqfN5bwRx91WWgUlTi49aM/+XrTfu47rxczzuzqkv0YhlE7THDUovKhEdE84sQ2WFIEy5+CX56H0M4wZQm0P7l2CluJjC+/Yv+jjxJw5hm0f/ppl/U9VVjs4OYP1/Fd3AEevKA3U0/v4pL9GIZRe0xw1JLNqZuZtmRa7YVG2k7rCvCkNTDoahj7NPgG1k5hK5H1w48k33sv/lFRRLz0EuJTu12vlyooLuGmD/5k6eYDPHxRH64f0dkl+zEMo3aZ4KgFpaHh7+VfO6GxYSF8dSeIB/xjPvSbUDsFrYKclStJuu02/Pr0IeKN1/Hwq/3uSsAKjRveX8cPWw7yr3F9ufbUTi7Zj2EYtc8ExwnakralLDSix0SfWGjkZ1iBseljiBwOE96G4A61V9hK5K1fz94bb8KnY6Q13Guga45w8otKmPn+WpZvPcQT4/sxcWhHl+zHMAzXMMFxArakbWHq91PLQqND8xP4kt/7B3w6BTKSYNSDcPod4OGadoWK5G/dyp7pM/Bq2ZIO0dF4hYS4Zj9FJUx7dw0rElJ4akJ/rhgS6ZL9GIbhOiY4aqg0NJp5NTux0HCUWI3fy5+CoPYw+VvoMKR2C1uJwl272DNlKh7NmhE5bx7erVu7ZD95hSVMfXc1v21P5elLT+KyqLo7mjIMo/aY4KgB59CYN2ZezUPj8B7rNNs9v0P/y+CC58AvqHYLW4mi5GR2T54MJSVEvvsOPhG13M27LbewmCkL1rBqZyrP/3MAE04+wXYgwzDcxgRHNW1N21o7oRG7CL64DdQB49+GAZfXajmrojglhT2Tp+DIyqbjOwvw7eKaU2FzCoq5fsFq1uxK48XLBzJuoGvCyTCMuuHSq8hEZKyIbBWRBBGZXcH8kSKSISLr7dv/Oc27VURiRSRORG5zmh4qIktEZJt975rK+ApsTdvKlO+nnFhoFGTD4pvgk+uhZXeY+YtbQqMkI4M9U6ZSdOAAHd56E78+fVyyn+yCYibN+4O1u9N56YpBJjQMoxFwWXCIiCfwGnAe0Ae4UkQq+nb6RVUH2rd/2ev2A6YBQ4ABwIUiUtpx0Wxgmap2B5bZz13uqNA4t4ahkbQW3jod1n8AZ9xttWeE1v21C46cHPZOn0Hhjh1EvPIK/ie75qLCrPwiro1exfq9h3nlykFcNMD1XaQYhuF6rjziGAIkqOoOVS0EPgLGVXHd3sBKVc1V1WLgJ46MODgOeMd+/A5wSe0VuWJ/C40W1QwNhwNWvAjR51rdoV/3FZz1IHh6u6bAxytKQQF7Z80ib9Mm2j3/HIGnjXDJfjLyirgm+g82Jmbw6lUnc37/cJfsxzCMuufK4GgP7HV6nmhPK+9UEdkgIt+ISF97WixwhoiEiYg/cD5Q+m3dRlX3Adj3rjkFyFbapuHn6Vez0MhMhvfGwdJHoNcFcMMK6OSaL+vKaFERSXfcSe7vKwl/8glanOuaoWUzcou4JnoVcckZvD7xZMb2a+uS/RiG4R6ubByvqEc8Lfd8HdBRVbNF5HxgMdBdVTeLyNPAEiAb2AAUV2vnItOB6QCRkTW7VqA0NHw9fZk/Zn71Q2Pzl9aQrsUFcPGrVtchbhqQSB0Oku9/gOxly2jz0IMuG73vcG4hV0ev4q/92bx59WBG93bNELaGYbiPK484EjlylAAQASQ7L6CqmaqabT/+GvAWkZb282hVPVlVzwDSgG32agdEJBzAvj9Y0c5V9W1VjVLVqFatWtXoBXy45cOahUZhrnXG1MKJENwRZvwCJ1/jvtBQZf+//kXmF1/Q6rbbCJ040SX7Scsp5Ko5q/jrQDZvXWtCwzAaK1cecawGuotIZyAJuAK4ynkBEWkLHFBVFZEhWEGWas9rraoHRSQSmACcaq/2OTAJeMq+/8xVL+CBYQ+QmpdK24BqVLXs22hdAZ7yF4y41boK3Ms1nQRW1aEXXuDwRwsJmzqFsBnTXbKP1OwCJs5dxc6UHOZcG8WZPWoW1oZh1H8uCw5VLRaRWcB3gCcwT1XjRGSmPf9N4B/ADSJSDOQBV6hqaXXWpyISBhQBN6lquj39KeB/IjIF2AP801WvwdvDu+qh4XDAqjestoxmoXDNYug6ylVFq7KUt94mdc5cgq+4nFZ33umSMTUOZRUwce5K9qTlEj3pFE7r3rLW92EYRv0hR76nG6+oqChds2aN63aQdQAW3wDbl0HP8632jIAw1+2vitI++IADjz1Oi4suot3TTyEuGPzpYFY+V81ZRVJ6HtHXRTG8qwkNw2gsRGStqkaVn26uHD9Rf30Hi2+Ewmy44AWImuy2tgxnhxcv5sBjjxM4ejTtnnzCJaFxIDOfK+esZH9GPvOvP4VhXdwfloZhuJ4Jjpoqyocl/wd/vAVt+sGl0dC6l7tLBUDmkiXsu/8B/E8dRvsXnke8a/96kf0ZVmgczMznnclDOKVTaK3vwzCM+skER00ciLdG5zsYB8NuhNEPg7drBjyqruxffyX5jjtp1r8/HV59FQ9f31rdvqryw5aDPPpFPGk5hbw7ZSiDO9ZZry+GYdQDJjiqQxVWz4XvHwTf5jDxE+h+jrtLVSZ33ToSZ92MT5cudHj7LTwCAmpt26rK9/EHeHnZNuKSM4kM9ee9KUMYFGlCwzCaGhMcVZWTAp/Ngr++gW5nwyVvQKBLL1qvlvz4ePbOmIl3mzZERs/FM6h2umd3OJTv4vbz8g8JbN6XSacwf5775wDGDWyHt6dL+8g0DKOeMsFRFdt/gJiZkJcOY5+CITPABY3NNVWwY4c1EFNgIJHzovFqeeJnNpU4lG9i9/HKsgS2HsiiS6sAXrx8ABed1A4vExiG0aSZ4Die4kL44V/w2yvQqhdc/Sm07e/uUh2lMDGJPddPBg8PIudF493uxHqgLXEoX25M5pUfEkg4mE231oG8dMVALjypHZ4e7j9bzDAM9zPBcTyf3wwbP4KoKXDu4+Dj7+4SHaXo4EH2TJ6MIy+Pju+9i2/nmnfRXlzi4As7MHYcyqFnm+a8etUgzu8XjocJDMMwnJjgOJ7Tboc+F1u92tYzxenp7J0yleKUFDrOi8avZ8+abafEQcyfSbz2YwK7UnPp1bY5b0w8mTF925rAMAyjQiY4jqd1r3pzbYazkmx7IKbdu+nw9ls0Gziw2tsoKnGwaF0ir/24nT1pufRt14K3rhnMOb3bmMAwDOO4THA0MI78fBJvuIH8+HgiXnmZgGHDqrV+YbGDT9Ym8tqPCSQdzuOkiCAeviiKs3q1dkk/VoZhND4mOBoQLSwk6dbbyF2zhnbPPEPzs86q8roFxSX8b00ib/yYQHJGPgM7BPP4+H6M7NHKBIZhGNVigqOB0JISkmfPJvunn2j76KMEXXRhldbLLyph4eq9vLF8O/sz8xncMYSnLj2J07u3NIFhGEaNmOBoAFSV/Y88QubX39D67rsIufyyStfJKyzhv3/s4a2ftnMwq4AhnUJ5/rIBDO8aZgLDMIwTYoKjnlNVDj79DIc//oSwmTMImzLluMvnFhbzwco9vPXzDlKyCzi1SxgvXTGIU7uanmsNw6gdJjjquZTXXydtwQJCrr6aVrfeeszlcgqKeW/lbub8vIPUnEJO69aSW0afzJDOptdawzBqlwmOeizt3XdJeeVVgsaPp83991VYxZSVX8S7v+9m7i87SM8t4owerbh1dDcGdzSBYRiGa5jgqKcOf/opB578N83PPZfwx/71t4GYMvOLeOfXXcxdsZOMvCJG9WzFLaO7m95qDcNwORMc9VDmt9+y76H/I2DECNo99yzideRjysgtYt6vO5n3606y8os5u3drbhndnZMigt1XYMMwmhQTHPVM9s8/k3T3PTQbNIiIV17Gw8cHgPScQub9upMFv+4iq6CYMX3bcPNZ3enXvna6TzcMw6gqExz1SO7q1STefAt+3bvT4c038PD3Jy2nkLm/7OCd33aRU1jC+f3bMmtUd/q0a+Hu4hqG0USZ4Kgn8jbFsnfmDXi3b0+HuXNIFx/mfLOZ937fTV5RCRf0D+fms7rTs21zdxe1XigqKiIxMZH8/Hx3F8UwGjw/Pz8iIiLw9vau0vImOOqBgm3b2DttGp7Bwfi/8gZP/baf91ftprDYwUUD2jFrVDe6tzGB4SwxMZHmzZvTqVMnc0GjYZwAVSU1NZXExEQ6V3FoBhMcbla4dy97Jk/B4eXF51ffx5z5cRSVOLhkYHtuOqsbXVsFuruI9VJ+fr4JDcOoBSJCWFgYhw4dqvI6JjjcqOjAAXZeex15OXncOeJGdm0rYMKg9tw0qhudWga4u3j1ngkNw6gd1f1fMsHhJnt3JpM06Vq801J54PQbOGXkYOaP7EZkWP0aZdAwDKM8Exx1bG9aLnO+3cQpLz9Ih8wDLJ/6IHOuu4gOoSYwDMNoGDwqX8SoDXtSc7n3k42Mfep7TnrtUbpk7iP4uRe47Y7LTWgYx1VQUMDZZ5/NwIEDWbhwYbXXX7x4MfHx8VVePjU1lVGjRhEYGMisWbOqvT+j7j8zgI0bN3LqqafSt29f+vfv79IzDs0Rh4vtTMnhtR8TiPkzCV8t4dUtCwlP3037F56nxdhz3F28RuHRL+KIT86s1W32adeChy/qW6vbrKk///yToqIi1q9fX6P1Fy9ezIUXXkifPn2qtLyfnx+PPfYYsbGxxMbG1miflXn6j6fZkralVrfZK7QX9w65t1a3WVN1/ZkVFxdz9dVX89577zFgwABSU1OrfGptTZgjDhfZfiib2xeuZ/Tzy/liQzLXDenA5xnfEf7XesIf+xctxo51dxGNE7Rr1y569erF1KlT6devHxMnTmTp0qWMGDGC7t2788cff/DHH38wfPhwBg0axPDhw9m6dSsAL7zwApMnTwZg06ZN9OvXj9zc3L/t4+DBg1x99dWsX7+egQMHsn37dtauXcuZZ57J4MGDGTNmDPv27QNgzpw5nHLKKQwYMIBLL72U3NxcfvvtNz7//HPuvvvusvUrExAQwGmnnYafn18tvlv1Q2P9zL7//ntOOukkBgwYAEBYWBienp619bb9nao2+tvgwYO1rvy1P1Nv/u867TT7S+314Df6+Jdxuv9wjibNvk/je/bS1AUL6qwsjVl8fLy7i6A7d+5UT09P3bhxo5aUlOjJJ5+s119/vTocDl28eLGOGzdOMzIytKioSFVVlyxZohMmTFBV1ZKSEj399NN10aJFOnjwYF2xYsUx9/Pjjz/qBRdcoKqqhYWFeuqpp+rBgwdVVfWjjz7S66+/XlVVU1JSytZ54IEH9OWXX1ZV1UmTJunHH39cNu+ZZ57RAQMG/O128803H7Xf+fPn60033XSib1O90lg/sxdffFGvvvpqPffcc3XQoEH69NNPV/u9qeh/ClijFXynmqqqWrJlfyav/JDA15v20czbk+lndGHa6V0IC/DhwJP/JiMmhpazZhE6aZK7i2rUos6dO9O/f38A+vbty+jRoxER+vfvz65du8jIyGDSpEls27YNEaGoqAgADw8PFixYwEknncSMGTMYMWJElfa3detWYmNjOeccq5qzpKSE8PBwAGJjY3nwwQc5fPgw2dnZjBkzpsJt3H333dx9990n+tIbrMb4mRUXF7NixQpWr16Nv78/o0ePZvDgwYwePbrK70t1mOA4QfHJmbzywza+id1PoK8XN47sypTTuhAaYHVOePCll0h/7z1Cr7uOljfd6ObSGrXN19e37LGHh0fZcw8PD4qLi3nooYcYNWoUMTEx7Nq1i5EjR5Ytv23bNgIDA0lOTq7y/lSVvn378vvvv/9t3nXXXcfixYsZMGAACxYsYPny5RVu49lnn+WDDz742/QzzjiDl19+ucplaaga42cWERHBmWeeScuWLQE4//zzWbduncuCw7Rx1FBsUgbT3l3D+S//woptKdxyVjdW3DuKu8f0KguN1Oh5pL7xJsH//Aet773HXLDWBGVkZNC+fXsAFixYcNT0W2+9lZ9//pnU1FQ++eSTKm2vZ8+eHDp0qOxLqKioiLi4OACysrIIDw+nqKjoqC+Z5s2bk5WVVfb87rvvZv369X+7NYXQqIqG+JmNGTOGjRs3kpubS3FxMT/99FOVG9ZrwqXBISJjRWSriCSIyOwK5o8UkQwRWW/f/s9p3u0iEicisSLyoYj42dMfEZEkp3XOd+VrKG/D3sNMWbCaC19Zwaodqdx2dndW3HsWd5zbk2B/n7Ll0hf+j4PPPkuL88+j7SOPmNBoou655x7uu+8+RowYQUlJSdn022+/nRtvvJEePXoQHR3N7NmzOXjwYKXb8/Hx4ZNPPuHee+9lwIABDBw4kN9++w2Axx57jKFDh3LOOefQq1evsnWuuOIKnn32WQYNGlSlhlaATp06cccdd7BgwQIiIiKqfWpoQ9YQP7OQkBDuuOMOTjnlFAYOHMjJJ5/MBRdcUINXXzVitX+4YMMinsBfwDlAIrAauFJV452WGQncpaoXllu3PbAC6KOqeSLyP+BrVV0gIo8A2ar6XFXLEhUVpWvWrDmh17NuTzovL9vG8q2HCGrmzdTTOjNpRCda+P39lLeML78i+e67CTjjdDq88gri41PBFo0TsXnzZnr37u3uYhhGo1HR/5SIrFXVqPLLurKNYwiQoKo77AJ8BIwDqvrTxQtoJiJFgD9Q9UrFWrRmVxovLdvGL9tSCPH35u4xPbn21I40ryAwALJ++JHk2bPxj4oi4qWXTGgYhtHouDI42gN7nZ4nAkMrWO5UEdmAFQx3qWqcqiaJyHPAHiAP+F5Vv3daZ5aIXAusAe5U1fTyGxWR6cB0gMjIyBq9gEc+j2PBb7sIC/Bh9nm9uGZYRwJ8j/2W5axcSdJtt+HXuzcRb7yBRyM8D95wnfnz5/PSSy8dNW3EiBG89tprbiqRUZmm+pm5sqrqn8AYVZ1qP78GGKKqNzst0wJwqGq23Vbxkqp2F5EQ4FPgcuAw8DHwiaq+LyJtgBRAgceAcFWdfLyy1LSqavnWg2w7kM3EYZH4+xw/Y/M2bGD39ZPxad+eyHffwSskpNr7M6rOVFUZRu2qTlWVKxvHE4EOTs8jKFfdpKqZqpptP/4a8BaRlsDZwE5VPaSqRcAiYLi93AFVLVFVBzAHq0rMJUb2bM20M7pUGhr5W7eyZ/oMvFq2pEP0XBMahmE0aq4MjtVAdxHpLCI+wBXA584LiEhbsU83EpEhdnlSsaqohomIvz1/NLDZXi7caRPjAdd0plNFhbt2sWfKVDz8/IicNw/v1q3dWRzDMAyXc1kbh6oWi8gs4DvAE5inqnEiMtOe/ybwD+AGESnGasu4wr7MfZWIfAKsA4qBP4G37U0/IyIDsaqqdgEzXPUaKlOUnMzuyZPB4SDy3XfwiWjvrqIYhmHUGZdeOW5XP31dbtqbTo9fBV49xroPAw9XMP2aWi5mjRSnpFhDvmZl0/GdBfh26eLuIhmGYdQJc+V4DZRkZLBn6jSKDhygw1tv4ufCKzQNo67Hdvjjjz8YOHAgAwcOZMCAAcTExFR7n01dfRpD5YEHHqBDhw4EBgZWuxzHYvqqqiZHTg57Z8ykcPt2It54A/+TT3Z3kYxvZsP+TbW7zbb94bynanebNVTXYzv069ePNWvW4OXlxb59+xgwYAAXXXQRXl6193Wx/8knKdhcu+Nx+PbuRdv776/VbdZUfRpD5aKLLmLWrFl07969RmWpiDniqAZHQQF7Z80ib9Mm2r3wPIGnVa13TKNxaqxjO/j7+5eFRH5+fqPqLqexfmbHG0Nl2LBhZb3x1pqK+lpvbLfaGI/DUVioe268SeN79tL0mJgT3p5xYsx4HK4dj2PlypXap08fDQgI0EWLFtXG21UvNObPTPX4Y6gEBAQc970x43HUMnU4SH7gAbKXLaPNQw8SfMkl7i6SUU80xrEdAIYOHUpcXBybN29m0qRJnHfeeY1mRMDG+pnVJRMclVBV9j/2GJmff0Gr224jdOJEdxfJqEca49gOznr37k1AQACxsbFERf3tAuIGqbF/ZnXBtHFU4tALL3L4w48ImzaVsBnT3V0co4FpiGM77Ny5k+LiYgB2797N1q1b6dSpU43fg4amIX5mdc0Ex3GkvD2H1DlzCL7yClrdcUejaiQ06kZDHNthxYoVZdseP348r7/+etnIck1BQ/zM4NhjqNxzzz1ERESQm5tLREQEjzzySDXejYq5rJPD+qSmnRxmfvst2T//QvjjjyEeJmPrE9PJoWHUrvoyHkeD12LsWFqMHevuYhiGYdQrJjgMo55oqmM7NGRN9TMzVVVGg7R582Z69epl2p0MoxaoKlu2bKkX43EYhsv4+fmRmppKU/jhYxiupKqkpqZW6zodU1VlNEgREREkJiZy6NAhdxfFMBo8Pz8/IiIiqry8CQ6jQfL29qZz587uLoZhNEmmqsowDMOoFhMchmEYRrWY4DAMwzCqpUmcjisih4Dd7i7HcbQEUtxdiHrKvDcVM+/LsZn3pmI1eV86qmqr8hObRHDUdyKypqJzpQ3z3hyLeV+Ozbw3FavN98VUVRmGYRjVYoLDMAzDqBYTHPXD2+4uQD1m3puKmffl2Mx7U7Fae19MG4dhGIZRLeaIwzAMw6gWExyGYRhGtZjgcDMRCRaRT0Rki4hsFpFT3V0mdxORniKy3umWKSK3ubtc9YWI3C4icSISKyIfikjVuzVtxETkVvs9iWvqfy8iMk9EDopIrNO0UBFZIiLb7PuQmm7fBIf7vQR8q6q9gAHAZjeXx+1UdauqDlTVgcBgIBeIcW+p6gcRaQ/cAkSpaj/AE7jCvaVyPxHpB0wDhmD9H10oIt3dWyq3WgCUH750NrBMVbsDy+znNWKCw41EpAVwBhANoKqFqnrYrYWqf0YD21W1Pl/5X9e8gGYi4gX4A8luLk990BtYqaq5qloM/ASMd3OZ3EZVfwbSyk0eB7xjP34HuKSm2zfB4V5dgEPAfBH5U0TmikiAuwtVz1wBfOjuQtQXqpoEPAfsAfYBGar6vXtLVS/EAmeISJiI+APnAx3cXKb6po2q7gOw71vXdEMmONzLCzgZeENVBwE5nMDhY2MjIj7AxcDH7i5LfWHXS48DOgPtgAARudq9pXI/Vd0MPA0sAb4FNgDFbi1UI2aCw70SgURVXWU//wQrSAzLecA6VT3g7oLUI2cDO1X1kKoWAYuA4W4uU72gqtGqerKqnoFVTbPN3WWqZw6ISDiAfX+wphsyweFGqrof2CsiPe1Jo4F4NxapvrkSU01V3h5gmIj4i4hg/c00+RMqAESktX0fCUzA/O2U9zkwyX48CfisphsyV467mYgMBOYCPsAO4HpVTXdroeoBu556L9BFVTPcXZ76REQeBS7Hqor5E5iqqgXuLZX7icgvQBhQBNyhqsvcXCS3EZEPgZFYXakfAB4GFgP/AyKxfoD8U1XLN6BXbfsmOAzDMIzqMFVVhmEYRrWY4DAMwzCqxQSHYRiGUS0mOAzDMIxqMcFhGIZhVIsJDqPBE5HsWtjGOSKyVkQ22fdnOc0TEfnB7lsMEVERed5p/l0i8sgxtrvL3mZpT78vH6cMI0VkuNPzmSJy7Ym+Nntb99fGduxtPef8/hhNj5e7C2AY9UQKcJGqJts9rX4HtLfnnQ9sUNVM+3kBMEFE/q2qKVXY9qgqLjcSyAZ+A1DVN6vzAipxP/BkVRe2Ly4UVXVUMPsVYA7wQy2VzWhgzBGH0SiJyEARWSkiG0UkpnTsARE5xZ72u4g8Wzpegar+qaqlvczGAX4i4ms/n8jRV9kWY43ffPsJlO8WEYm3y/KRiHQCZgK320cmp4vIIyJyl738chF5UUR+tsdtOUVEFtljKzzutN3F9hFTnIhMt6c9hdWb7noR+cCedoc9dkVs6dgVItLJ3vbrwDqgg4gssJfZJCK32+/VbiBMRNrW9PUbDZyqmpu5NegbkF3BtI3AmfbjfwH/sR/HAsPtx08BsRWs+w9gqdPz3UBz5/0BLYBdQBBwF/DIMcq2C9gErLdvt9vTkwFf+3Gwff8IcJfTumXPgeXA0/bjW+31wwFfrD7Pwux5ofZ9M/u1hpV/j7DGONkEBACBWEE5COgEOIBhTsstcVov2OnxHOBSd3/25uaemzniMBodEQnC+pL7yZ70DlaX28FYAfCbPf2/FazbF6uX1RlOk0NVNct5ObWqrd7FGlSpMqPUHphKVV+0p20EPrB7tq1qL66f2/ebgDhV3adWVyM7ONKF+C0isgFYaU+raDCj04AYVc1R1WysjhJPt+ftVtWV9uMdQBcReUVExgKZTts4iNU7r9EEmeAwmhI57kyRCKyRBq9V1e1Os4pFpKL/lf8AU7B+uSMink6N4P+qpCwXAK9h/apfaw/KVJnS/qgcTo9Ln3uJyEis3nNPVdUBWP1YVTSs7PHeh5zSB2r1mTYA62jnJqw+1Ur5AXlVKLPRCJngMBodtTpFTBeR0l/R1wA/2V+EWSIyzJ5eNuSqfTTyFXCfqv5abpNbsQbdKr+fNKxO46bYz0ucjiz+71jls0Oog6r+CNwDBGNVGWUBzav5cp0FAemqmisivYBhTvOKRMTbfvwzcIndw24A1kh5v1RQzpaAh6p+CjzE0V3+98CqCjOaIHNWldEY+ItIotPzF7C6jX7T7mV3B3C9PW8KMEdEcrB+SZf2vDsL6AY8JCIP2dPOVdWDWIEyEkioYN/P2+sez48iUmI/3miX4X27Sk2AF1X1sIh8AXwiIuOAmyt/2X/zLTBTRDZihd1Kp3lvAxtFZJ2qThSRBcAf9ry5qvqn3UDvrD3W6JSlPzDvA7ADqBuwpgZlNBoB0zuu0aSISKBdr4+IzAbCVfXWStYJB95V1XPqooz1nYiMB05W1YcqXdholMwRh9HUXCAi92H97e8GrqtsBVXdJyJzRKSFHrmWoynzwjrSMpooc8RhGIZhVItpHDcMwzCqxQSHYRiGUS0mOAzDMIxqMcFhGIZhVIsJDsMwDKNa/h+DGqqjwHG15wAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"#We'll plot in this block\n",
"\n",
"x = np.log2(np.array(n_est))\n",
"for m in m_feat:\n",
" plt.plot(x, aucs_oob[m], label='max_feat={}'.format(m))\n",
" \n",
"plt.title('OOB AUC by Max Feat and N-Estimators')\n",
"plt.xlabel('Log2(N-Estimators)')\n",
"plt.ylabel('OOB-AUC')\n",
"plt.legend(loc=4, ncol=2, prop={'size':10})\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
We can see in the plot above that the forest gets better with more trees, but we see that this effect tapers off asymptotically. We also see somewhat of a Golidlocks phenomenon with respect to the optimal number of features. Having too few features is suboptimal (i.e., max_feat=1), probably due to high bias in the trees. Conversely, having too many features (max_feat=11) is also bad, and this is because the individual trees are too correlated with each other, which prevents the Bagging approach from reducing the overall variance (see the analysis below to understand the math behind this statement). The right number of features is in the middle. We see that choosing $3$ or $6$ perform relatively the same. $3$ happens to be closer to the good heuristic default, which is max_feat=sqrt(total feats). \n",
"
\n",
"Next, for the purpose of educating ourselves, we'll also look at the Test set AUC's on each of the design options above.\n",
"
We can see similar patterns in the true holdout evaluation as we did in the OOB evaluation. We also tend to see a few more things. For one thing, having more trees actually hurts the performance for some levels of max_feat. We also see that the best choice for max_feat was indeed $3$.\n",
"\n",
"\n",
"
Much like with Decision Trees, the Random Forest Classifier has a built in mechanism for evaluating feature importance. Quoting the sklearn documentation:
\n",
"\n",
"Features used at the top of the tree contribute to the final prediction of a larger fraction of the input samples. The expected fraction of the samples they contribute to can thus be used as an estimate of the relative importance of the features.\n",
"
\n",
"The above computation is made for each feature in each tree and then averaged over all trees in the forest. The Random Forest Classifier returns an attribute with an importance score for each feature, and these scores sum to $1$ across all features.\n",
"
\n",
"We'll train a RF on the best options found above and then look at the feature importances.\n",
"\n",
"
We can show with a little algebra why Random Forests lead to generally better estimation performance.
Let's say that the variance of a single tree of size $N$ and $p$ features is: $Var(T_b(X))=\\sigma^2$. Similarly the variance of the Random Forest is the variance of the sum of such trees, i.e., $Var(RF(X))=Var(\\frac{1}{B}\\sum\\limits_{b=1}^B \\: T_b(X)).$ Because the invidual trees are trained with overlapping records and features, the individual tree estimates are correlated. Thus, the variance does not factor completely cleanly into a sum of individual tree variances. We'll solve this by using the following relation:
\n",
"This last quantity shows the two factors that reduce the variance of the Random Forest. 1). The term $\\rho\\sigma^2$ shrinks to $0$ as the correlation between trees reduces to $0$, which is achieved by randomly downsampling the features. 2). The right-hand term shinks asymptotically to $0$ as $B$ increases.\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}