{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem Sheet 7 - Subset selection" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this and the following exercises, we want to explore different methods of linear model selection.\n", "\n", "In the lecture you've learned about problems that might occur in datasets with many predictors (high $p$)and a low number of samples (low $n$).\n", "Ways out could be:\n", "* Subset selection - try to find a suitable subset of predictors\n", "* Skrinkage/Regularization - increase weights of *important* predictors, decrease weights of *unimportant* ones\n", "* Dimension reduction - Build linear combinations $v_i, i=1,\\ldots,M$ of predictors and fit a model using these vectors instead of predictors with $M < p$\n", "\n", "We always have to keep in mind that it's in general not wise to select the model with the minimal training error, due to the danger of overfitting.\n", "Our goal is to find a model that performs well on a test set.\n", "This refers to a subset of samples that are completely held out from training (and also cross-validation)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 1 - Best subset selection" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to implement the **best subset selection** algorithm from the lecture:\n", "1. Let $\\mathcal{M}_0$ denote the *null model*, which contains no predictors.\n", "2. For $k = 1, 2, \\ldots, p$:\n", " 1. Fit all $p \\choose k$ models that contain exactly $k$ predictors.\n", " 2. Pick the best among these and call it $\\mathcal{M}_k$, while the best is the one with highest $R^2$ score.\n", "3. Select a single best model from among $\\mathcal{M}_0, \\ldots \\mathcal{M}_p$ using one of the following methods:\n", " * cross-validated prediction error\n", " * $C_p$ (or equivalently AIC - Akaike information criterion)\n", " * BIC - Bayesian information criterion\n", " * adjusted $R^2$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1.1 - Understanding steps 1 and 2\n", " \n", "The algorithm `bestSubsetComputation` belows contains the implementation of step 1 and step 2.\n", "\n", "**Task**: Understand the following code and add comments as you wish." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "from itertools import combinations\n", "\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.metrics import r2_score\n", "\n", "def bestSubsetComputation(X, y, scoring_func = r2_score):\n", " \"\"\" Input: X - predictor array of size (n,p)\n", " y - array of size (n,)\n", " scoring_func - function that takes two arguments y_true\n", " and y_pred, and returns a score\n", " \n", " \"\"\"\n", " \n", " # Get the number of samples and the number of predictors \n", " n, p = X.shape\n", " \n", " # Prepare lists that keep the best scores and models:\n", " # best_score[i] keeps the best score in a model i predictors\n", " # best_model[i] keeps the best model using i predictors\n", " \n", " best_score = []\n", " best_model = []\n", "\n", " ### First step in best subset selection algorithm\n", " \n", " # The model containing no predictors simply predicts the sample mean\n", " ybar = y.mean()\n", " yhat = ybar * np.ones_like(y)\n", "\n", " best_score.append(scoring_func(y, yhat))\n", " best_model.append( () )\n", " \n", " ### Second step in best subset selection algorithm\n", "\n", " # Loop over k - number of predictors in our model\n", " for k in range(1,p+1):\n", "\n", " best_model.append( () )\n", " best_score.append(0.)\n", "\n", " for l in combinations(range(p),k):\n", "\n", " lr = LinearRegression()\n", " lrfit = lr.fit(X[:,l],y)\n", " yhat = lrfit.predict(X[:,l])\n", "\n", " this_score = scoring_func(y,yhat)\n", " \n", " if this_score > best_score[k]:\n", " best_score[k] = this_score\n", " best_model[k] = l\n", " \n", " return (best_model, best_score)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**: Load the diabetes data set form sklearn. If you forgot how to do this have a look at problem sheet 6.\n", "Store the predictors in an array `X`, and the targets in an array `y`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from sklearn.datasets import load_diabetes\n", "data = load_diabetes()\n", "\n", "#print(data.DESCR)\n", "X = data.data\n", "y = data.target\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**: Use the function `train_test_split` from `sklearn.model_selection` to split the data into a training and test set. Use as `test_size=0.2` and `random_state=1`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import train_test_split\n", "X_train, X_test, y_train, y_test = train_test_split(\n", " X, y, test_size=0.2, random_state=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**:\n", "Apply the `bestSubsetComputation` function from above to the training set." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "best_model, best_score = bestSubsetComputation(X_train, y_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**: Now, you can execute and interpret the output of the following cell." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Score of model with 0 predictors has score 0.0000\n", "\tSelected predictors ()\n", "\n", "Score of model with 1 predictors has score 0.3632\n", "\tSelected predictors (2,)\n", "\n", "Score of model with 2 predictors has score 0.4741\n", "\tSelected predictors (2, 8)\n", "\n", "Score of model with 3 predictors has score 0.4918\n", "\tSelected predictors (2, 3, 8)\n", "\n", "Score of model with 4 predictors has score 0.5053\n", "\tSelected predictors (2, 3, 6, 8)\n", "\n", "Score of model with 5 predictors has score 0.5272\n", "\tSelected predictors (1, 2, 3, 6, 8)\n", "\n", "Score of model with 6 predictors has score 0.5302\n", "\tSelected predictors (1, 2, 3, 4, 6, 8)\n", "\n", "Score of model with 7 predictors has score 0.5320\n", "\tSelected predictors (1, 2, 3, 4, 5, 7, 8)\n", "\n", "Score of model with 8 predictors has score 0.5329\n", "\tSelected predictors (1, 2, 3, 4, 5, 7, 8, 9)\n", "\n", "Score of model with 9 predictors has score 0.5332\n", "\tSelected predictors (0, 1, 2, 3, 4, 5, 7, 8, 9)\n", "\n", "Score of model with 10 predictors has score 0.5332\n", "\tSelected predictors (0, 1, 2, 3, 4, 5, 6, 7, 8, 9)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAHjCAYAAACaZwbkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3X2c1XWd9/HXh2FAMVQMzBsoxEVbwRJDlCtMYxGkGy27WSzLrXZtK3Ntq6291pKsru3m2q1217a87P7K3K6yjbo8B5BUvHRDwRtU0ELUwDtIEYXQcZjv9cfvTAwwwAzM7/zO+Z3X8/GYx+/8zvnNOe8ZRpk33+/5fiOlhCRJkiSpXAYVHUCSJEmSNPAse5IkSZJUQpY9SZIkSSohy54kSZIklZBlT5IkSZJKyLInSZIkSSVk2ZMkSZKkErLsSZIkSVIJWfYkSZIkqYQGFx2gv0aOHJnGjh1bdAxJkiRJKsSyZct+n1Iatafrmq7sjR07lqVLlxYdQ5IkSZIKEREP9+U6p3FKkiRJUglZ9iRJkiSphCx7kiRJklRClj1JkiRJKiHLniRJkiSVkGVPkiRJkkrIsidJkiRJJWTZkyRJkqQSsuxJkiRJUglZ9iRJkiSphCx7kiRJklRClj1JkiRJKiHLniRJkiSVkGVPkiRJkkrIsidJkiRJJWTZkyRJUvOaO7foBM3F71f/NPn3K1JKRWfol8mTJ6elS5cWHUOSJCkfc+c2/S+YdRUBTfb7bL91f30pbfvoed6fxw46CDZurE/evNXjdUaMaMifr4hYllKavKfrBtcjjCRJamGWl0xK8MIL0Nm5++NnPgNvfjN0dW372Lq199t7Om+FawGOPXZgC9GuHsv7+Xt7LA8HHZTfc6uhOLInSZLy1ZeRl66uPZegvTnm8Zx7+1rdxaRZRUBbGwwalH30vL2n84G+9v774b77ds44cSIcf3yWtfujO/uOt/f2sYF4jnq/9oIFcN11O3+/Zs6EWbN6//PeV90Z8pbH61SrMH/+zvdfemnD/MOVI3uSJOWlmUequrrg+eezj+eeyz66b+/puLfXAowbt/tSVEQRGjQI2tth8OD+HdvbYf/99+5zd3X85S/h5z/fOeOcOXDeeftWkAaiaPUsEI2mFaZx7quPf3zbbb9fe3bxxdtuN/n3y5E9SZL6a2/+8u/qgo6O+hWsXR1feGFgvgf77QdDh+76+Mgj8OCDO3/eiSfCKacMXEkaPHjvP2dQg65T1+S/XNad36/+8fvVPw36/XJkT5KkgbJ1azZlbMkSuPPO7L53vKN/RaujY2Cy7K5g7bcfDBuWLSiwpzLW89ifa4cOhSFD+jfK06C/LKkkLr206ATNxe9X/zT598uRPUmSdvToo1mxu/XW7Lh0KTz7bO/XHn44HH30wBeq3j6nvyWrUVj2+qeZpwlLqgtH9iRJ6otNm7Iy113sbr0V1q7NHmtvh1e+Et71Ljj55Oxj/PjsfUyWl75r8n8ZrzuLnqQBYtmTJLWOzk5YsSIrdd3F7t57ty0OcvTRcOqp24rdCSdko2raN5YXSSqEZU+SVE4pZSN0PYvd0qXwhz9kjx9yCEyZAueckxW7k06CkSP79tyOVEmSmoBlT5JUDs88A7fdtm065pIl8Pjj2WNDhsCkSfCXf5kVvJNPzkbx9vb9b45USZKagGVPktR8XngB7r57+/fZrVy57X10xxwDZ5yxrdi94hXZAieSJLUQy54kqbGlBA8/vP10zGXLtm3WPWpUVujmzNk2HXPEiGIzS5LUACx7kqTGsmFDNh2zu9jdeiusW5c9tt9+2abcH/hAVuymTIGxY5tzOwJJknJm2ZMkFaejA+66a/v32f3mN9ljEfDyl8PrXret2B1/fLYdgiRJ2iPLniSpPlKCBx7YvtjdcUdW+AAOOywrdeefnx0nT4aDDio2syRJTcyyJ0nKx5NPbl/sbr0Vnnoqe2zYsKzMXXTRtj3tRo92OqYkSQPIsidJ2nfPPQd33rmt1C1Zko3iAQwaBBMmwJvfvG065oQJMNi/giRJypN/00qS+qerC3772+2L3V13ZdshQDZCN2UKXHBBdnzVq2D48GIzS5LUgix7ktTq5s7d/Sbh69ZtPx3zttvg6aezx170omyrg49+dNuedkccUY/UkiRpDyJ1b0DbJCZPnpyWLl1adAxJKo+IbZuR/+EP2aIpPUftHnooe6ytLVsNs3sq5sknZ6tltrUVFl2SpFYUEctSSpP3dJ0je5LUirq64PHHYfXq7PwDH8iK3fLlsHVrdt/LXpaVugsvzIrdiSdmC6tIkqSmYNmTpLJ69ll48MGs0HUfuz9++9ttpQ7gG9/IjqeeCh//eDY187DDisktSZIGRK5lLyLOBL4GtAFXppS+sMPjfwF8GXikdte/pZSuzDOTJJVGZyesXdt7mXvwQVi/fvvrDzwQjj46WwnzDW+AceOyjzPPzIrfoEHFfB2SJCkXuZW9iGgDLgfOANYCt0XEvJTSih0u/Y+U0oV55ZCkppUSbNiwfYHrWeh+97us8HUbPDibejluXLbNQXeZGzcOjjoKRozY9T52Fj1Jkkonz5G9KcCqlNJqgIi4Gjgb2LHsSVLrev55ePjhXY/Obdy4/fWjRmXlbcoUmDNn+zI3evTe7V136aUD87VIkqSGkmfZOxJY0+N8LXByL9e9JSJeA/wG+EhKac2OF0TEBcAFAC996UtziCpJOUkJnnhi12Vu7dptK2ECDB26rbxNm7btdvcxj/3qdrftgiRJalp5lr3e5grtuM/DL4AfpZSej4i/Br4HTN/pk1K6ArgCsq0XBjqoJO2TzZuz4tZbmVu9GrZs2f76I4/MittrX7v9yNy4cdmiKE6plCRJAyDPsrcWGNPjfDTwaM8LUkpP9jj9X8AXc8wjSXtn61Z45JFdl7knntj++he9KCtu48fDrFnbity4cTB2LOy3XyFfhiRJai15lr3bgPERcRTZaptzgHf0vCAiDk8pPVY7PQtYmWMeSa1i7tz+T018+umdp1p2337oIXjhhW3XDhoEL31pVt7e+Mbty9y4cfDiF+96IRRJkqQ6iZTymxUZEa8Dvkq29cK3U0qfj4jLgKUppXkR8Y9kJa8TeAr4QErpvt095+TJk9PSpUtzyyypBCK2fx8cQEdHtnrlrkbnNmzY/vpDDtl5imX3x5gx0N5ev69HkiSph4hYllKavMfr8ix7ebDsSdql55+H++6DE06Az31u+zK3Zg10dW27dsiQbEplb2XuqKPgoIMK+zIkSZJ2p69lL9dN1SUpF52d8MADcM892z6uvx6e7PE24EsuyY6jR8Npp+1c5o44AtraiskvSZJUB5Y9SY0rpWxE7u67ty92K1dmo3iQTdk8+mg49VSYODH7mDMHNm2CAw4oNr8kSVKBLHuSGsO6ddsXuu6PZ5/dds3o0VmZmzFjW7H70z+FYcO2f645cyx6kiSp5Vn2JNXXxo1w7707l7r167dd8+IXZ0Xu/PO3lboJE+Dgg/v2Gpdemk92SZKkJmLZk5SPLVuy6ZY7lro1a7Zdc8ABWZE766xtpW7iRHjJS/Zt64L+brsgSZJUQpY9SfvmhRdg1aqdS92qVdtWvxwyJJtu+ZrXbF/qXvrSbM86SZIkDTjLnqS+6eqChx/eudTdd1+2hx1kxe1P/mTbIikTJ8Lxx2f3DfZ/N5IkSfXkb1+StpcSPPHEzitg3nsvbN687bqXvjQrc2eeuW2k7uUvh/33Ly67JEmS/siyJ7WyDRt6Xyyl5351o0Zlo3Pve9+2UnfccW46LkmS1OAse1Ir+MMfYMWKnUvdI49su2b48KzInXPO9u+rO/TQ4nJLkiRpr1n2pDLp6IDf/GbnUrd6dTY9E2Do0Gxkbvr07UvdmDH7tgKmJEmSGoplT2oGc+duv51AVxc8+ODOpe7++7PVMQHa2mD8eDjxRHj3u7eVunHjXCxFkiSpBfgbn9ToOjvhM5/Jpll2l7oVK7Kpmd3Gjs2K3BvesK3UHXss7LdfYbElSZJULMue1Mg6O+Etb8luf+xjcNhhWZG74ILtF0sZPrzYnJIkSWo4lj2pUV16KVx22fb3Pf44vP/920/plCRJknph2ZMaUUrw7LPZ7c98Jit+3QusSJIkSX0wqOgAknrxP/4HfOUr8OEPw6c+VXQaSZIkNSHLntRo/v3f4ZJL4Lzz4KtfzbZDuPTSolNJkiSpyVj2pEZy9dXwoQ9lq2p++9swqPafqO/RkyRJUj9Z9qRGUa3Cu94F06bBj38M7e1FJ5IkSVITs+xJjeCWW7ItFiZMgF/8Avbfv+hEkiRJanKWPalod98Nr389HHEEzJ8PBx1UdCJJkiSVgGVPKtLq1TBzJgwbBgsXwkteUnQiSZIklYT77ElFeewxOOMM6OiAxYth7NiiE0mSJKlELHtSETZsgFmz4IknYNGi7L16kiRJ0gCy7En19oc/ZFsr3Hcf/N//CyefXHQiSZIklZBlT6qnjg5461vhv/4r217hjDOKTiRJkqSSsuxJ9dLVBX/xF1CpwBVXZKVPkiRJyomrcUr1kBJ8+MPwox/BP/4j/NVfFZ1IkiRJJWfZk+rh0kvh61+Hj30MPvGJotNIkiSpBVj2pLx97Wvw2c/Ce98LX/oSRBSdSJIkSS3Asifl6Qc/gIsvhje/Gb75TYueJEmS6sayJ+XlF7+A97wHpk+Hq66Cwa6HJEmSpPqx7El5WLwY3v52mDQJ/vM/Yb/9ik4kSZKkFmPZkwbaHXfAG98IY8dm2ywMH150IkmSJLUgy540kH7zG5g1Cw46CBYsgJEji04kSZKkFmXZkwbK2rUwc2Z2e+FCGDOm2DySJElqaa4YIQ2EJ5/MRvSeegpuuAGOPbboRJIkSWpxlj1pXz37LLzudfDAA1CtwoknFp1IkiRJsuxJ++T557M99JYtg5/+FE4/vehEkiRJEmDZk/be1q3wznfCokXw3e/C2WcXnUiSJEn6IxdokfZGSvDXf52N5v3zP8P55xedSJIkSdqOZU/aG3//93DllfAP/wAf+UjRaSRJkqSdWPak/vryl+GLX8xG9j772aLTSJIkSb2y7En98a1vwd/9Hfz5n8O//RtEFJ1IkiRJ6pVlT+qra66BCy7I9tP7/vehra3oRJIkSdIuWfakvli0CM49F04+OVuUZciQohNJkiRJu2XZk/bkttvgTW+CY46BX/4SDjig6ESSJEnSHln2pN1ZuRJmz4ZRo2DBAjjkkKITSZIkSX1i2ZN25eGH4YwzYPBgWLgQDj+86ESSJElSnw0uOoDUkNaty4repk2weDEcfXTRiSRJkqR+sexJO9q4Ec48E9auzUb0XvGKohNJkiRJ/WbZk3rasgXOOgvuvhvmzYNXv7roRJIkSdJesexJ3To7Yc4cuOkm+OEPs4VZJEmSpCZl2ZMAurrgfe/LRvMuvzzbU0+SJElqYq7GKaUEH/0ofP/7cNll8MEPFp1IkiRJ2meWPenzn4evfhUuugguuaToNJIkSdKAsOyptf37v8OnPgXnnQdf+QpEFJ1IkiRJGhCWPbWuq6+GD30I3vhG+Pa3YZD/OUiSJKk8/O1WralahXe9C049Ff7jP6C9vehEkiRJ0oCy7Kn13HILnHMOHH98tvrm/vsXnUiSJEkacJY9tZbly+H1r4fRo7PRvYMOKjqRJEmSlAvLnlrH6tUwaxYccAAsWACHHlp0IkmSJCk3bqqu1vDYY3DGGdDRATfdBGPHFp1IkiRJypVlT+W3YUM2ovfEE/CrX8FxxxWdSJIkScpdrtM4I+LMiLg/IlZFxCd3c91bIyJFxOQ886gFbd4Mb3gD3H8//Od/wpQpRSeSJEmS6iK3shcRbcDlwGzgOODciNhpSCUihgMXAUvyyqIW1dEBb30r/PrXcNVVMGNG0YkkSZKkuslzZG8KsCqltDql1AFcDZzdy3WfBb4EPJdjFrWarVvh/POzFTe/8Q14y1uKTiRJkiTVVZ5l70hgTY/ztbX7/igiJgFjUkq/3N0TRcQFEbE0IpauX79+4JOqXFKCiy6Cq6+GL3wB/uqvik4kSZIk1V2eZS96uS/98cGIQcBXgI/u6YlSSleklCanlCaPGjVqACOqlC69FL7+dfj4x+ETnyg6jSRJklSIPMveWmBMj/PRwKM9zocDE4EbIuIh4BRgnou0aJ987Wvw2c/C+94HX/xi0WkkSZKkwuRZ9m4DxkfEURExBJgDzOt+MKW0MaU0MqU0NqU0Fvg1cFZKaWmOmVRm3/8+XHwxnHNO9j696G1wWZIkSWoNuZW9lFIncCEwH1gJ/DildG9EXBYRZ+X1umpR8+bBe98Lf/Zn2cqbg91CUpIkSa0t19+IU0rXAtfucN+nd3Ht6XlmUYndeCO8/e1w4onws5/B0KFFJ5IkSZIKl+um6lLu7rgDzjoLxo2Da6+F4cOLTiRJkiQ1BMuemtdvfgOzZsHBB8OCBTByZNGJJEmSpIZh2VNzWrsWzjgju71wIYweXWweSZIkqcG4ioWaz+9/DzNnwoYNcMMNcMwxRSeSJEmSGo5lT83l2Wfhda+D1ath/vxsURZJkiRJO7HsqXk8/zy8+c1w++1wzTVw2mlFJ5IkSZIalmVPzWHrVnjnO2HRIvje97IVOCVJkiTtkgu0qPGlBO9/P/z0p/CVr8C73110IkmSJKnhWfbU+D75SfjWt+CSS+Dii4tOI0mSJDUFy54a25e+lH188INw2WVFp5EkSZKahmVPjevKK+ETn4A5c+Bf/xUiik4kSZIkNQ3LnhrTNddk79M788xsQZZB/qhKkiRJ/eFv0Go8110H554Lp5wCP/kJDBlSdCJJkiSp6Vj21FhuvRXe9CY49lj45S/hgAOKTiRJkiQ1JcueGseKFTB7NrzkJTB/PowYUXQiSZIkqWlZ9tQYHn4YZs7MpmwuXAiHH150IkmSJKmpDS46gMS6dXDGGbB5M9x4I4wbV3QiSZIkqelZ9lSsjRuzFTfXrs1G9F7xiqITSZIkSaVg2VNxtmyBs86Cu++GX/wCXv3qohNJkiRJpWHZUzE+9Sm46y646Sa46qpsdE+SJEnSgLHsqf66uuBzn8tuf/3rMGdOsXkkSZKkEnI1TtXfv/xLdvzsZ+EDHyg2iyRJklRSlj3Vz9y5EAEf+Uh2/qlPZedz5xaZSpIkSSqlSCkVnaFfJk+enJYuXVp0DO2tJ5+EQw/NpnI22c+eJEmS1AgiYllKafKernNkT/W1cGFW9CRJkiTlygVaVF/VKhxyCHzoQ0UnkSRJkkrNkT3VT1dXVvZmzoTLLis6jSRJklRqlj3Vz113wRNPuKeeJEmSVAeWPdVPtZodZ80qNockSZLUAix7qp9qFSZNgsMOKzqJJEmSVHqWPdXHxo1w881O4ZQkSZLqxLKn+li0CLZuhdmzi04iSZIktQTLnuqjWoUDD4RTTik6iSRJktQSLHvKX0pQqcCMGdDeXnQaSZIkqSVY9pS/FStg7VqncEqSJEl1ZNlT/txyQZIkSao7y57yV6nAhAkwZkzRSSRJkqSWYdlTvjZtgptucgqnJEmSVGeWPeXrhhugo8P99SRJkqQ6s+wpX5UKHHAATJtWdBJJkiSppVj2lJ/uLRemT4ehQ4tOI0mSJLUUy57ys2oVPPigUzglSZKkAlj2lJ9KJTta9iRJkqS6s+wpP9UqHHMMjBtXdBJJkiSp5Vj2lI8tW7KVOB3VkyRJkgph2VM+Fi/OCp9lT5IkSSqEZU/5qFZhv/3g9NOLTiJJkiS1JMue8lGtwmmnwf77F51EkiRJakmWPQ28hx6C++5zCqckSZJUIMueBl61mh1nzy42hyRJktTCLHsaeJUKjB2bbbsgSZIkqRCWPQ2sjg5YtCibwhlRdBpJkiSpZVn2NLBuvhk2b3YKpyRJklQwy54GVqUC7e3w2tcWnUSSJElqaZY9DaxqFU49FYYPLzqJJEmS1NIsexo4jzwCd9/tlguSJElSA7DsaeB0b7lg2ZMkSZIKZ9nTwKlW4cgjYeLEopNIkiRJLc+yp4HR2QkLF7rlgiRJktQgLHsaGL/+NWzc6BROSZIkqUFY9jQwqlVoa4MZM4pOIkmSJAnLngZKtQpTp8LBBxedRJIkSRKWPQ2EJ56AZcucwilJkiQ1EMue9t2CBdlx9uxic0iSJEn6I8ue9l21CoceCiecUHQSSZIkSTW5lr2IODMi7o+IVRHxyV4e/+uIuDsi7oyI/xcRx+WZRznYuhXmz4dZs2CQ/3YgSZIkNYrcfjuPiDbgcmA2cBxwbi9l7qqU0vEppROALwH/nFce5WTZMnjySadwSpIkSQ0mz6GYKcCqlNLqlFIHcDVwds8LUkrP9Dg9AEg55lEeqtVsE/Uzzig6iSRJkqQeBuf43EcCa3qcrwVO3vGiiPgQ8LfAEGB6jnmUh0oFTjoJRo4sOokkSZKkHvo0shcR+0fEsf187ujlvp1G7lJKl6eUjgY+AVyyi9e/ICKWRsTS9evX9zOGcvPkk3DrrU7hlCRJkhrQHsteRLwRuBOo1s5PiIh5fXjutcCYHuejgUd3c/3VwJt6eyCldEVKaXJKafKoUaP68NKqi+uug64u99eTJEmSGlBfRvbmkr3/7mmAlNKdwNg+fN5twPiIOCoihgBzgO1KYkSM73H6euC3fXheNYpKBQ45JJvGKUmSJKmh9OU9e50ppY0Rvc3K3LWUUmdEXAjMB9qAb6eU7o2Iy4ClKaV5wIURMQN4AdgAnN+/+CpMV1e2OMvMmdDWVnQaSZIkSTvoS9m7JyLeAbTVRuIuAm7py5OnlK4Frt3hvk/3uP03/ciqRrJ8OTzxhFM4JUmSpAbVl2mcHwYmAM8DVwEbgYvzDKUmUKlkx1mzis0hSZIkqVe7HdmrbYz+mZTSx4F/qE8kNYVqFSZNgsMOKzqJJEmSpF7sdmQvpbQVeFWdsqhZbNwIt9ziFE5JkiSpgfXlPXt31LZa+D/A5u47U0rX5JZKjW3RIujstOxJkiRJDawvZe8Q4Elgeo/7EmDZa1XVKhx4IEydWnQSSZIkSbuwx7KXUnpPPYKoSaSUlb0ZM6C9veg0kiRJknZhj6txRsToiPhZRKyLiCci4qcRMboe4dSAVqyANWucwilJkiQ1uL5svfAdYB5wBHAk8IvafWpF1Wp2tOxJkiRJDa0vZW9USuk7KaXO2sd3gVE551KjqlZhwgQYM6boJJIkSZJ2oy9l7/cRcV5EtNU+ziNbsEWtZtMmWLwYZs8uOokkSZKkPehL2Xsv8HbgceAx4K21+9RqbrgBOjqcwilJkiQ1gb6sxvk74Kw6ZFGjq1bhgANg2rSik0iSJEnag76sxvm9iDi4x/mIiPh2vrHUcFKCSgWmT4ehQ4tOI0mSJGkP+jKN8xUppae7T1JKG4BJ+UVSQ1q1ClavdgqnJEmS1CT6UvYGRcSI7pOIOIQ+TP9UyVQq2dGyJ0mSJDWFvpS2fwJuiYif1M7fBnw+v0hqSNUqHHMMjBtXdBJJkiRJfbDHkb2U0veBtwBPAOuAc1JKP8g7mBrIli3ZSpyO6kmSJElNY48jexFxNPBASmlFRJwOzIiIR3u+j08lt3hxVvgse5IkSVLT6Mt79n4KbI2IPwGuBI4Crso1lRpLtQr77Qenn150EkmSJEl91Jey15VS6gTOAb6WUvoIcHi+sdRQqlU47TTYf/+ik0iSJEnqo76UvRci4lzg3cAva/e15xdJDeWhh+C++5zCKUmSJDWZvpS99wBTgc+nlB6MiKOA/51vLDWMajU7zp5dbA5JkiRJ/bLHBVpSSiuAi3qcPwh8Ic9QaiDVKowdm227IEmSJKlp9GVkT62qowMWLcqmcEYUnUaSJElSP1j2tGs33wybNjmFU5IkSWpClj3tWrUK7e3w2tcWnUSSJElSP+2y7EVEW0S8PyI+GxGv3uGxS/KPpsJVKjBtGgwfXnQSSZIkSf20u5G9bwKnAU8C/xIR/9zjsXNyTaXiPfII3H23UzglSZKkJrW7sjclpfSOlNJXgZOBF0XENRExFHC1jrKbPz87ur+eJEmS1JR2V/aGdN9IKXWmlC4A7gR+Bbwo72AqWKUCRx4JEycWnUSSJEnSXthd2VsaEdsN66SULgO+A4zNM5QK1tkJCxe65YIkSZLUxHZZ9lJK56WUqr3cf2VKqT3fWCrUkiWwcaNTOCVJkqQmtsetFyKirR5B1EAqFWhrgxkzik4iSZIkaS/ttuxFxHDg53XKokZRrcLUqXDwwUUnkSRJkrSXdrfP3uHAdcAV9Yujwq1bB8uWOYVTkiRJanKDd/PYTcDHU0rz6hVGDcAtFyRJkqRS2N00zg3AkfUKogZRrcKhh8KkSUUnkSRJkrQPdlf2TgdmR8SH6pRFRdu6NRvZmzULBu1x7R5JkiRJDWx3Wy9sBs4CHOJpFcuWwZNPwuzZRSeRJEmStI929549Ukpbgb+sUxYVrVrNNlE/44yik0iSJEnaR/2eqxcRbRHxzjzCqGDVKpx0EowcWXQSSZIkSftod1svHBgRfx8R/xYRMyPzYWA18Pb6RVRdPPUULFniFE5JkiSpJHY3jfMHZCty/hfZVM6PA0OAs1NKd9Yhm+pp4ULo6nLLBUmSJKkkdlf2xqWUjgeIiCuB3wMvTSk9W5dkqq9qFQ45JJvGKUmSJKnp7e49ey9036gt1PKgRa+kurqysjdzJrS1FZ1GkiRJ0gDY3cjeKyPimdrtAPavnQeQUkoH5p5O9bF8OTz+uFM4JUmSpBLZZdlLKTnE0yqq1ew4a1axOSRJkiQNmH5vvaASqlRg0iQ47LCik0iSJEkaIJa9VrdxI9xyi1M4JUmSpJKx7LW6X/0KOjste5IkSVLJWPZaXaUCBx4IU6cWnUSSJEnSALLstbKUssVZZsyA9vai00iSJEkaQJa9VrZiBaxZ4xROSZIkqYQse62se8sFy54kSZJUOpa9VlatwoQJMGZM0UkkSZIkDTDLXqvatAkWL3ZUT5IkSSopy16ruuEG6OiA2bOLTiJJkiQpB5a9VlWtwrBhMG1a0UkkSZIk5cCy14pSyvbXmz4dhg4tOo0kSZKkHFj2WtHZUzh5AAAWwElEQVSqVbB6tVM4JUmSpBKz7LUit1yQJEmSSs+y14oqFRg/HsaNKzqJJEmSpJxY9lrNli3ZSpxO4ZQkSZJKzbLXam66KSt8TuGUJEmSSs2y12oqlWwFztNOKzqJJEmSpBxZ9lpNtQqnn57tsSdJkiSptHItexFxZkTcHxGrIuKTvTz+txGxIiKWR8SiiHhZnnla3kMPwX33OYVTkiRJagG5lb2IaAMuB2YDxwHnRsRxO1x2BzA5pfQK4CfAl/LKI7ZtueDiLJIkSVLp5TmyNwVYlVJanVLqAK4Gzu55QUrp+pTSH2qnvwZG55hH1SqMHQvHHFN0EkmSJEk5y7PsHQms6XG+tnbfrrwPqPT2QERcEBFLI2Lp+vXrBzBiC+nogEWLsimcEUWnkSRJkpSzPMteb40i9XphxHnAZODLvT2eUroipTQ5pTR51KhRAxixhdx8M2za5BROSZIkqUUMzvG51wJjepyPBh7d8aKImAH8A3BaSun5HPO0tmoV2tvhta8tOokkSZKkOshzZO82YHxEHBURQ4A5wLyeF0TEJOCbwFkppXU5ZlG1CtOmwfDhRSeRJEmSVAe5lb2UUidwITAfWAn8OKV0b0RcFhFn1S77MvAi4P9ExJ0RMW8XT6d98cgjsHy5UzglSZKkFpLnNE5SStcC1+5w36d73J6R5+urZv787Oj+epIkSVLLyHVTdTWIahWOPBImTiw6iSRJkqQ6seyVXWcnLFzolguSJElSi7Hsld2SJfD0007hlCRJklqMZa/sqlVoa4MZvj1SkiRJaiWWvbKrVGDqVDj44KKTSJIkSaojy16ZrVsHy5Y5hVOSJElqQZa9MluwIDta9iRJkqSWY9krs0oFDj0UJk0qOokkSZKkOrPsldXWrdlm6rNmwSD/mCVJkqRWYwsoq9tvhyefdAqnJEmS1KIse2VVqWSbqM+cWXQSSZIkSQWw7JVVtQonnQQjRxadRJIkSVIBLHtl9NRTsGSJUzglSZKkFmbZK6OFC6GrC2bPLjqJJEmSpIJY9sqoWoURI7JpnJIkSZJakmWvbLq6srI3cya0tRWdRpIkSVJBLHtls3w5PP64UzglSZKkFmfZK5tqNTu65YIkSZLU0ix7ZVOpwAknwOGHF51EkiRJUoEse2WycSPccotTOCVJkiRZ9krlV7+Czk7315MkSZJk2SuVSgUOPBCmTi06iSRJkqSCWfbKIqVscZYZM6C9veg0kiRJkgpm2SuLlSthzRqncEqSJEkCLHvlUalkR8ueJEmSJCx75VGtwoQJMGZM0UkkSZIkNQDLXhls3gyLFzuqJ0mSJOmPLHtlcP310NHh/nqSJEmS/siyVwbVKgwbBtOmFZ1EkiRJUoOw7JVBtQrTp8PQoUUnkSRJktQgLHvN7re/hQcecAqnJEmSpO1Y9ppdtZodXZxFkiRJUg+WvWZXrcL48TBuXNFJJEmSJDUQy14ze+65bCVOp3BKkiRJ2oFlr5ktXgxbtjiFU5IkSdJOLHvNrFrNVuA87bSik0iSJElqMJa9ZlapwOmnZ3vsSZIkSVIPlr1m9dBDcN99TuGUJEmS1CvLXrOaPz87WvYkSZIk9cKy16wqFRg7Fo49tugkkiRJkhqQZa8ZdXTAokXZqF5E0WkkSZIkNSDLXjO65RbYtMkpnJIkSZJ2ybLXjCoVaG+H6dOLTiJJkiSpQVn2mlG1CtOmwfDhRSeRJEmS1KAse83m0Udh+XKncEqSJEnaLctes6lWs+Ps2cXmkCRJktTQLHvNplqFI46AiROLTiJJkiSpgVn2mklnJyxc6JYLkiRJkvbIstdMliyBp592CqckSZKkPbLsNZNqFdraYMaMopNIkiRJanCWvWZSqcDUqXDwwUUnkSRJktTgLHvNYt06WLbMLRckSZIk9Yllr1ksWJAdLXuSJEmS+sCy1ywqFTj0UJg0qegkkiRJkpqAZa8ZbN0K8+fDrFkwyD8ySZIkSXtmc2gGt98OTz7pFE5JkiRJfWbZawaVSraJ+syZRSeRJEmS1CQse82gWoWTToKRI4tOIkmSJKlJWPYa3VNPwZIlTuGUJEmS1C+WvUa3cCF0dcHs2UUnkSRJktRELHuNrlqFESOyaZySJEmS1EeWvUaWUlb2Zs6Etrai00iSJElqIpa9RnbXXfD4407hlCRJktRvlr1GVq1mR7dckCRJktRPlr1GVq3CCSfA4YcXnUSSJElSk8m17EXEmRFxf0SsiohP9vL4ayLi9ojojIi35pml6TzzDNx8s1M4JUmSJO2V3MpeRLQBlwOzgeOAcyPiuB0u+x3wF8BVeeVoWosWQWen++tJkiRJ2iuDc3zuKcCqlNJqgIi4GjgbWNF9QUrpodpjXTnmaE7VKhx4IEydWnQSSZIkSU0oz2mcRwJrepyvrd3XbxFxQUQsjYil69evH5BwDS0lqFRgxgxoby86jSRJkqQmlGfZi17uS3vzRCmlK1JKk1NKk0eNGrWPsZrAypWwZo1TOCVJkiTttTzL3lpgTI/z0cCjOb5eeXRvuWDZkyRJkrSX8ix7twHjI+KoiBgCzAHm5fh65VGpwIQJMGbMnq+VJEmSpF7kVvZSSp3AhcB8YCXw45TSvRFxWUScBRARJ0XEWuBtwDcj4t688jSNzZth8WJH9SRJkiTtkzxX4ySldC1w7Q73fbrH7dvIpneq2w03QEeHZU+SJEnSPsl1U3XthUoFhg2DU08tOokkSZKkJmbZazTVKkyfDkOHFp1EkiRJUhOz7DWSVavggQecwilJkiRpn1n2Gkmlkh1nzy42hyRJkqSmZ9lrJNUqjB8P48YVnUSSJElSk7PsNYrnnoPrr3dUT5IkSdKAsOw1isWLYcsW368nSZIkaUBY9hpFtZqtwHnaaUUnkSRJklQClr1GUa3C6adne+xJkiRJ0j6y7DWChx+GlSudwilJkiRpwFj2GkG1mh0te5IkSZIGiGWvEVSrMHYsHHts0UkkSZIklYRlr2gdHXDdddmoXkTRaSRJkiSVhGWvaLfcAps2OYVTkiRJ0oCy7BWtUoH2dpg+vegkkiRJkkrEsle0ahWmTYPhw4tOIkmSJKlELHtFevRRWL7cKZySJEmSBpxlr0jdWy7Mnl1sDkmSJEmlY9krUrUKRxwBEycWnUSSJElSyVj2itLZCQsXuuWCJEmSpFxY9oqyZAk8/bRTOCVJkiTlwrJXlGoV2tpgxoyik0iSJEkqIcteUapVOOUUOPjgopNIkiRJKiHLXhHWrYOlS53CKUmSJCk3lr0iLFiQHd1fT5IkSVJOLHtFqFbh0ENh0qSik0iSJEkqKctevXV1wfz5MGsWDPLbL0mSJCkfto16W7YMfv97p3BKkiRJypVlr96q1WwT9Zkzi04iSZIkqcQse/VWqcBJJ8HIkUUnkSRJklRilr16euopWLLEKZySJEmScmfZq6frrssWaHF/PUmSJEk5s+zVU6UCI0Zk0zglSZIkKUeWvXpJKVucZeZMaGsrOo0kSZKkkrPs1cvy5fD4407hlCRJklQXlr16qVSyo1suSJIkSaoDy169VKtwwglw+OFFJ5EkSZLUAix79fDMM3DzzU7hlCRJklQ3lr16WLQIOjvdX0+SJElS3Vj26qFahQMPhKlTi04iSZIkqUVY9vLWveXCjBnQ3l50GkmSJEktwrKXt5Ur4Xe/cwqnJEmSpLqy7OWtWs2Olj1JkiRJdWTZy1u1ChMmwJgxRSeRJEmS1EIse3navBluvNFRPUmSJEl1Z9nL0w03QEeHZU+SJElS3Vn28lStwrBhcOqpRSeRJEmS1GIse3mqVGD6dBg6tOgkkiRJklqMZS8vq1bBAw84hVOSJElSISx7eenecmH27GJzSJIkSWpJlr28VCowfjyMG1d0EkmSJEktyLKXh+eeg+uvdwqnJEmSpMJY9vKweDFs2eIUTkmSJEmFsezloVrNVuA87bSik0iSJElqUZa9PFSrWdEbNqzoJJIkSZJalGVvoD38MKxc6RROSZIkSYWy7A207i0XXJxFkiRJUoEsewOtWoWXvQyOPbboJJIkSZJamGVvIHV0wHXXZVM4I4pOI0mSJKmFWfYG0i23wKZNTuGUJEmSVDjL3kCqVqG9HaZPLzqJJEmSpBZn2RtIlQpMmwbDhxedRJIkSVKLs+wNlEcfheXLncIpSZIkqSFY9gbK/PnZ0f31JEmSJDUAy95AqVTgiCNg4sSik0iSJElSvmUvIs6MiPsjYlVEfLKXx4dGxH/UHl8SEWPzzJObzk6YNy+bwumWC5IkSZIaQG5lLyLagMuB2cBxwLkRcdwOl70P2JBS+hPgK8AX88qTq1tvheefdwqnJEmSpIaR58jeFGBVSml1SqkDuBo4e4drzga+V7v9E+DPIppwaKxSyY4zZhSbQ5IkSZJq8ix7RwJrepyvrd3X6zUppU5gI/DiHZ8oIi6IiKURsXT9+vU5xd0Lc+dm0zY/97nsfMSI7Hzu3CJTSZIkSVKuZa+3Ebq0F9eQUroipTQ5pTR51KhRAxJuQMydCyllH7DttmVPkiRJUsHyLHtrgTE9zkcDj+7qmogYDBwEPJVjJkmSJElqCXmWvduA8RFxVEQMAeYA83a4Zh5wfu32W4FfpZR2GtlrCpdeWnQCSZIkSfqjwXk9cUqpMyIuBOYDbcC3U0r3RsRlwNKU0jzgW8APImIV2YjenLzy5M6pm5IkSZIaSG5lDyCldC1w7Q73fbrH7eeAt+WZQZIkSZJaUa6bqkuSJEmSimHZkyRJkqQSsuxJkiRJUglZ9iRJkiSphCx7kiRJklRClj1JkiRJKiHLniRJkiSVkGVPkiRJkkrIsidJkiRJJWTZkyRJkqQSsuxJkiRJUglZ9iRJkiSphCx7kiRJklRClj1JkiRJKqFIKRWdoV8iYj3wcNE5ejES+H3RIVRa/nwpT/58KW/+jClP/nwpT4368/WylNKoPV3UdGWvUUXE0pTS5KJzqJz8+VKe/PlS3vwZU578+VKemv3ny2mckiRJklRClj1JkiRJKiHL3sC5ougAKjV/vpQnf76UN3/GlCd/vpSnpv758j17kiRJklRCjuxJkiRJUglZ9iRJkiSphCx7AyAizoyI+yNiVUR8sug8Ko+IGBMR10fEyoi4NyL+puhMKp+IaIuIOyLil0VnUblExMER8ZOIuK/2/7GpRWdSeUTER2p/N94TET+KiP2KzqTmFhHfjoh1EXFPj/sOiYiFEfHb2nFEkRn7y7K3jyKiDbgcmA0cB5wbEccVm0ol0gl8NKX0p8ApwIf8+VIO/gZYWXQIldLXgGpK6eXAK/HnTAMkIo4ELgImp5QmAm3AnGJTqQS+C5y5w32fBBallMYDi2rnTcOyt++mAKtSSqtTSh3A1cDZBWdSSaSUHksp3V67/SzZL0pHFptKZRIRo4HXA1cWnUXlEhEHAq8BvgWQUupIKT1dbCqVzGBg/4gYDAwDHi04j5pcSmkx8NQOd58NfK92+3vAm+oaah9Z9vbdkcCaHudr8Zdx5SAixgKTgCXFJlHJfBX4O6Cr6CAqnXHAeuA7tWnCV0bEAUWHUjmklB4B/ifwO+AxYGNKaUGxqVRSL0kpPQbZP8IDhxacp18se/suernP/Sw0oCLiRcBPgYtTSs8UnUflEBFvANallJYVnUWlNBg4Efj3lNIkYDNNNv1Jjav2vqmzgaOAI4ADIuK8YlNJjceyt+/WAmN6nI/GaQQaQBHRTlb0fphSuqboPCqVVwNnRcRDZFPQp0fE/y42kkpkLbA2pdQ9G+EnZOVPGggzgAdTSutTSi8A1wD/reBMKqcnIuJwgNpxXcF5+sWyt+9uA8ZHxFERMYTszcHzCs6kkoiIIHu/y8qU0j8XnUflklL6+5TS6JTSWLL/d/0qpeS/jGtApJQeB9ZExLG1u/4MWFFgJJXL74BTImJY7e/KP8MFgJSPecD5tdvnAz8vMEu/DS46QLNLKXVGxIXAfLKVoL6dUrq34Fgqj1cD7wLujog7a/f995TStQVmkqS++jDww9o/hq4G3lNwHpVESmlJRPwEuJ1s5eo7gCuKTaVmFxE/Ak4HRkbEWuBS4AvAjyPifWT/yPC24hL2X6Tk28skSZIkqWycxilJkiRJJWTZkyRJkqQSsuxJkiRJUglZ9iRJkiSphCx7kiRJklRClj1JUq4iIkXEP/U4/1hEzB2g5/5uRLx1IJ5rD6/ztohYGRHX5/w6YyPintrtyRHxL3u4/r/nmUeS1Nwse5KkvD0PnBMRI4sO0lNEtPXj8vcBH0wpvbYOrwVASmlpSumiPVzWr7IXGf/ul6QW4f/wJUl56yTb7PgjOz6w48hcRGyqHU+PiBsj4scR8ZuI+EJEvDMibo2IuyPi6B5PMyMibqpd94ba57dFxJcj4raIWB4R7+/xvNdHxFXA3b3kObf2/PdExBdr930amAZ8IyK+vMP1p0fE4oj4WUSsiIhvdJepiNgUEZdFxBJgakS8qvY1LYuI+RFxeO26V0XEXRHxX8CHdnjuX9ZuvygivlPLtjwi3hIRXwD2j4g7I+KHtev+tpb9noi4uHbf2Nqo5NfJNqAeU/u+31N7vp3+XCRJ5TC46ACSpJZwObA8Ir7Uj895JfCnwFPAauDKlNKUiPgb4MPAxbXrxgKnAUcD10fEnwDvBjamlE6KiKHAzRGxoHb9FGBiSunBni8WEUcAXwReBWwAFkTEm1JKl0XEdOBjKaWlveScAhwHPAxUgXOAnwAHAPeklD4dEe3AjcDZKaX1EfHnwOeB9wLfAT6cUrpxxzLZw6dqX8/xtawjUko/jYgLU0on1O57FfAe4GQggCURcWPtazkWeE9K6YO1645MKU2sfd7Bu/wTkCQ1NUf2JEm5Syk9A3wf2NO0xJ5uSyk9llJ6HngA6C5rd5MVvG4/Til1pZR+S1YKXw7MBN4dEXcCS4AXA+Nr19+6Y9GrOQm4IaW0PqXUCfwQeE0fct6aUlqdUtoK/IhsFBBgK/DT2u1jgYnAwlqmS4DREXEQcHBK6cbadT/YxWvMICvMAKSUNvRyzTTgZymlzSmlTcA1wKm1xx5OKf26dns1MC4i/jUizgSe6cPXKElqQo7sSZLq5atk0wi/0+O+Tmr/8BgRAQzp8djzPW539TjvYvu/v9IOr5PIRrY+nFKa3/OBiDgd2LyLfLHHr6B3vb0+wHO1Atj93PemlKbukOfgXj5/V9n2dN3u8v/xa04pbYiIVwKzyKaNvp1shFGSVDKO7EmS6iKl9BTwY7LFTro9RDZtEuBsoH0vnvptETGo9j6+ccD9wHzgA7Xpk0TEMRFxwB6eZwlwWkSMrC2oci7Z1Ms9mRIRR9Xeq/fnwP/r5Zr7gVERMbWWpz0iJqSUngY2RkT3aOA7d/EaC4ALu08iYkTt5gvdXyOwGHhTRAyrfa1vBm7a8YlqC+UMSin9lGx66Il9+BolSU3IsidJqqd/Anquyvm/yArWrWTvNdvVqNvu3E9WyirAX6eUngOuBFYAt0e2lcE32cNslpTSY8DfA9cDdwG3p5R+3ofX/y/gC8A9wIPAz3p57g7grcAXI+Iu4E7gv9Uefg9weW2Bli27eI3PASNqi6rcBXSvCnoF2Xshf5hSuh34LnArWXG9MqV0Ry/PdSRwQ2066XdrX7MkqYQipb7MHpEkSTuqTQv9WErpDUVnkSRpR47sSZIkSVIJObInSZIkSSXkyJ4kSZIklZBlT5IkSZJKyLInSZIkSSVk2ZMkSZKkErLsSZIkSVIJ/X/vlnAZpSuJgAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "plt.rcParams['figure.figsize'] = (15,8)\n", "plt.plot(range(len(best_score)), best_score, 'r+-')\n", "plt.xlabel('Number of predictors')\n", "plt.ylabel('R^2 score')\n", "\n", "for i,s in enumerate(best_score):\n", " print('\\nScore of model with %i predictors has score %6.4f' % (i,s))\n", " print('\\tSelected predictors', best_model[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Observation**:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Observation**: We see that the $R^2$ score increases with the number of predictors.\n", "The value improves only insignificantly for $p > 5$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1.2 - Implementing step 3 of the *best subset selection* algorithm\n", "In the following tasks, you should implement step 3 of the *best subset selection* algorithm using the training data from the diabetes data set." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**: Implement at least one of the performance criteria discussed in the lecture and above **as a function** that takes as arguments `y`, `yhat`, `d`, and `shat`, if necessary:\n", "* AIC\n", " $$ AIC = \\frac{1}{n \\hat \\sigma^2} (RSS + 2 d {\\hat \\sigma}^2)$$\n", "* BIC\n", " $$ BIC = \\frac{1}{n} (RSS + \\log(n) d {\\hat \\sigma}^2)$$\n", "\n", "* Adjusted R^2\n", " $$ R^2_{Adj} = 1 - \\frac{RSS / (n - d - 1)}{TSS / (n - 1)} $$\n", " \n", "with $\\hat{\\sigma}^2$ referring to an estimate of the variance associated with each response in the linear model (estimated on a model containing all predictors):\n", " $$ {\\hat \\sigma}^2 = \\frac{1}{n - p - 1} \\sum_{i=1}^n (y_i - \\hat y_i)^2. $$\n", " \n", "**Remember**: $TSS$ is the total sum of squares, which is defined by $\\sum_{i=1}^n (y_i - \\bar y)^2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def RSS(y, yhat): return np.power(y-yhat,2).sum()\n", "def TSS(y): return np.power(y-y.mean(),2).sum()\n", "\n", "# Put your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def RSS(y, yhat): return np.power(y-yhat,2).sum()\n", "def TSS(y): return np.power(y-y.mean(),2).sum()\n", "\n", "# Put your code here\n", "def sigmaHat(X, y):\n", " n,p = X.shape\n", " lr = LinearRegression()\n", " lrfit = lr.fit(X, y)\n", " yhat = lrfit.predict(X)\n", " return np.sqrt(1./(n-p-1) * RSS(y, yhat))\n", "\n", "def AIC(y, yhat, d, shat):\n", " n = len(y)\n", " return (RSS(y,yhat) + 2 * d * shat**2) / (n * shat**2)\n", "\n", "def BIC(y, yhat, d, shat):\n", " n = len(y)\n", " return (RSS(y,yhat) + np.log(n) * d * shat**2) / (n)\n", "\n", "def adjustedRSquare(y, yhat, d):\n", " n = len(y)\n", " return 1 - (RSS(y,yhat) / (n - d - 1)) / (TSS(y) / (n - 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Task**:\n", "1. Use the return values of the function `bestSubsetComputation` to perform step 3 of the *best subset selection* algorithm using the training data from the diabetes data set.\n", "2. Plot your chosen score against the number of predictors.\n", "3. Compute and mark the minimizer ($AIC$, $BIC$) or maximizer ($R^2$) in the plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'BIC')" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA38AAAHiCAYAAABY0HU7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3XmcXFWZ//Hvk05n6SSdpKtKREKIYARlDKABMgoDISaCsshvcEZAFheCoyDOOMoyahJcQMTRURFE2VQWGcQxMmzFJshmwjKBsGgMxISwdLqTNEknnV6e3x/3VlNpqtd033uq+vN+vfrVXbdO3Xq6kpzc7z3n3mPuLgAAAABAZRuRdgEAAAAAgKFH+AMAAACAYYDwBwAAAADDAOEPAAAAAIYBwh8AAAAADAOEPwAAAAAYBgh/FcTMrjazb6b03mZmV5nZejP7UwrvP83M3MxGJvR+88zsf5J4r8FiZl8wswvTrgPDE/0T/VNPzOxoM7sh7TpQGehvEu1vTjezHyTxXuXAzP7TzD6bdh09IfwNITN70cxeNbNxRds+Y2b3pVjWUDlI0lxJU9z9gLSL6Q8zO9jMNnX5cjP7xx5e9m1JFxbtY5qZ3WtmzWb2nJl9sIf3u9rMtnV5v6r4uVlmljezRjOrN7P/NrOdi177ZTN72sxeN7MXzOzLXfb9opltKdrvnUVPXy7pE2b2lv5+Rqg89E/lw8yqzOybZrY2/rf/hJlN6uElXfunb5jZU2bWZmYLe3kvM7PvmFlD/HWRmVnR8/ua2WNxX/eYme1b9NwkM7vGzF6LvxZ22Xe3dbj7Ykl/Z2Yz+vKZoLzQ35SP+Phnc9FxxM97aDtK0lclfTd+nDWzB+O+Y4OZPWxmHyhq/3dmdoeZrTOzNy00bmZnmNlSM2sxs6t7eN8FcZ09HWvta2YPmNlGM1tjZl8vrtvMbor/XrqZHdrltQvNrLXLcdruffyMvivpP+LPJkiEv6E3UtJZaRfRX4Uw0g+7SXrR3TcPwnubmSX2d9PdH3D38YUvSUdK2iTp9m7q21/SRHd/pGjz9ZKekJSR9B+SbjKzXA9ve1Hxe7p7e7x9sqKQNk3RZ/q6pKuK317SyXG7wyWdYWYf77Lvo4r2O6/o99wq6bb49YBE/zSQ9060f4otkvR+SX8vqVbSSZK2lmrYTf+0QtJXJP1vH95rvqSPStpH0gxF/eHp8b5HSfqdpF8p6oOukfS7ooOc70uqUdR/HSDpJDP7ZD/quD5+f1Qm+pv+v3ca/Y0k7VN0HPGZHtodI+k5d38pfrxJ0qck5RT1Ed+R9Ht7YxSyVdKNkj7dzf7WSvqmpCu7e0Mz20PScZJe7uV3uE7S/ZLqJB0i6V/M7Oii5/8o6ROSXunm9b/ucpy2ssvzJT8jd39Z0nOSjlagCH9D77uS/r3UWVorMTRvZveZ2Wfin0+Nz6B8Pz6DstLM3h9vXx2fWT2ly26zFo0cvW5mfzCz3Yr2vZe9Mar0vJn9U9FzV5vZpWZ2q5ltljS7RL1vM7PF8etXmNlp8fZPS/q5pL+Pz4AsKvHawu/yo/gszHNmNqfL7/0tM3tQUrOk3c1sopldYWYvm9lLFp35LoyQVZnZxfHZo5WSPtK3P44+OUXSTT103EdI+kNR7e+U9F5JC9x9i7v/RtJTknoaOSzJ3W9z9/929yZ3b5b0Y0kfKHr+Ind/3N3b3P15RQdiH+hufyXcp8H9rFDe6J8Udv9kZpMlfVHSae6+yiNPxydzStmuf5Ikd7/G3W9TdDKpN6dI+p67r4kP6L4n6dT4uUMVHcD/wN1b3P2Hik5IHRY/f5SiE1vN7v6ipCsUHQj2tY77RP9UyehvFHZ/MwDb9TfuvtXdn3f3DkV9Q7uiEFgXP/+8u18haXmpnbn7ze7+P5IaenjPH0s6W9K2XmqbJulad293978qCnt7x++zzd1/4O5/jGscbPcp4L6M8Df0lir6S/DvA3z9gZKWKRpRuk7SDZL2l/QORWcsfmxm44vanyjpG5Kykp6UdK0kWTTVIh/v4y2Sjpf0EzPbu+i1J0j6lqQJiv6RdHW9pDWS3qborMu3zWxO/A/5s5Iejs+ALOjhd1kZ17ZA0s1mVlf0/EmKzvpOkLRK0Vnltvh33U/SPEmFsyunKTojvZ+kmXE93TKzW8zsnJ7axO1q4n1d00Oz90h6vujx3pJWunvxAc3/xdu787n4P43HrOfppf+gbjpJMzNJB5d4/lqLpozeaWb7dHnuWUVn9AGJ/qnr7xJi//Se+H2OM7NXzOzPZvb5HnbXtX/qr70V9V8FxX3Z3pKWuXvxdK1l2r6vsy4//10/3vtZSdPMrLYfr0H5oL/Z/ncJsb8puD/ub242s2k9tCvZ35jZMkWzExZL+rm7v9bL+/WJmX1M0jZ3v7UPzX8g6WQzqzazPRXNnLirH293VHycttzM/qXE8z19RkEfaxH+kvF1SWdaz9MAu/OCu18VTwv8taRdJZ0fn3W9U9GZj3cUtf9fd7/f3VsUTT/8ezPbVVHH8GK8rzZ3f1zSb7R9J/E7d3/Q3Tu6nlWO93GQpLPjMztPKjq7dVI/fpfXFJ0xbnX3XyvqMIrPjFzt7svdvU3RWaIjJH3R3TfHHcf3JRWmOP5TvK/V7t4o6YKe3tjdj3T3vtzs5B8lrVOXM+ddTNL2Z67HS9rYpc1GRZ12KT+UNF3Rfzpfk3S1Fc2JL7Do2pevS/py1+diCxX9Gy6eFnqi3pgyeq+kO7qcZX1d0sRu9ofhif4pEmr/NEXRv9l3Snq7os9koZnN7aZ91/6pv7r2ZxsljY9PNvXW190u6Rwzm2Bm71A06lfTj/cu1N3T9Ywob/Q3kVD7GymaIjlN0l6KpmHeYt3fPKZkf+PuMxRNUT9BpcNzv8XB/tuKZkL0xS2K/ky3KJqGeYW7L+nja2+U9C5F01dPk/R1Mzu+6PnePqPXFXA/RvhLgLs/regvYa8jTyW8WvTzlnh/XbcVn+laXfS+myQ1KjoztZukA+PpEhvMbIOioPDWUq8t4W2SGruMbq2StEs/fpeXupwxXhXvt9T77yapWtLLRfX+VFFgKtRT3H5VP+roySmSftGlzq7Wa/tgt0lRJ1esVt0cgHk0bbMh/k/nVkVnI/9fcZv4wOk2SWe5+wNd92FmZyi6du8j8X9shX0/6NHU02Z3v0DSBkWjgwUT9OaDNwxj9E+dQu2ftsTfz4//bS9TNOLx4W7ad+2f+qtrf1YraVP82fTW130hrvcviqakF0ZH+qpQ94Z+1owyQX/TKdT+RnFg3ubuGxRdo/l2RUGolG77mzgYX6/ohNBgjIItkvRLd3+ht4bxKOrtks6XNEbRiYIPmdnn+vJG7v6Mu6/1aMroQ5L+S0UnB/rwGU1QwP0Y4S85CxSdPSjuHArXlBWfGS3ufAZi18IP8VmSOkVnJVZL+oO7Tyr6Gu/uxUPZPQWetZLqzKz4H/lUSS91076UXeKzx8WvX9vN+6+W1CIpW1RvrbsXpmW8rKLfNd7XDonP5h0q6Re9NF2m6Cx8wXJFc/KLP5t91M10zRJcRVOlLLou4S5J33D3X5ao81OK/uOc4+69HVhtt29FndP/ddMWwxf9U7j907IS799b+3f22qp7y7X9dKXivmy5pBldPqcZhefdvdHdT3T3t8afxQhJ/bnV/bsUjcg0Dbh6lAP6m3D7m1K6HkcU60t/Uy1p917a9MUcSV+Ip1q+ouh3vtHMzi7RdndJ7e7+i/hE+xr1fNKsNz19BqWeD/pYi/CXEHdfoWiawheKttUr6iw+EV+w+ylJe+zgW33YzA6y6O5r35D0qLuvVnSm7Z1mdlI8/7nazPY3s+7O5nStf7WkhyRdYGZj4imJn1Y8h76P3qLoH251PG/7XZJKztv26G5Jd0r6npnVmtkIM9vDzA6Jm9wY72uKRTdEGMhZxK5OkvSQRxcG9+RWRUP+hVr/rOh6ggXxZ3OsogOi35R6sZkdZ2bj499pnqJrFRbHz+0i6R5Jl7j7ZSVee6KiaQ9zvcudp8xsqpl9wKJbGI+xaBmIrKQHi5odomhEEehE/yQp0P4p7o8eUHTr8NHxZ/LPij6zUrbrnyQp/p3GKPo/f2T8GXV3B8NfSPo3M9vFzN4m6UuSro6fu0/RzRG+ENdyRrz9nvh99jCzTPz35QhF1yx1rrXWhzron4YB+htJgfY3Zra3RUskVMWB+XuK/lye7eYl2/U3Fi1XdVB8HDI2DmY7SXo0ft7iPmBU/HiMmY0uev3I+PkqSVXx84XplHMUXUO8b/y1VtGdiC8pUdef47c7If683qqo3+wMZHEfNiZ+WDhusvi5Y8xsclzvAYr+rv6uH59R0H0Z4S9Z50sa12XbaYqu6WpQdNH8Qzv4HtcpOqvWKOl9iqYyKJ6eME/RHPG1im5t+x1Jo0vvpqTjFc1xXivpt4rubpnvx+sfVXSt2zpFF1If5+493dHpZEUdxDOKphbcJKmw5t3PJN2h6B/y45JuLn6hmV1mZpcVPb7NzM7rpb6T1fONXiRF0zYlbTSzA4s2f1zRhdbrFa2vdVz8n5nM7EQzKx4FPEtRR7FB0d3PTnP3++LnPqPojNUCK1pfpui131R0sfuSoucLv+cESZfGNbykaCmIIwqfcdzJfbgvvyOGJfqncPun4xVN/WpQtEzC19z97lINu+mffqZoStzxiq592qL4+iSL1zktavtTSb9XdMfip+P3+2m8722KloE4WVH/9SlJH423S9Gf6VOKpoFeIOlEdy/u+7qto+j3/GkPnwMqB/1NmP3NToqCeZOiG9JMk3Sku7d20/73kvaKTxRJ0Wd4iaI/w5cUHXN8xN0Lo5q7Kfp3X+gXtmj7G8Z8Nd52jqIT41vibYovl3ml8KXoRNT6eErvdr9nPHvg/0n61/jzelJRf/atovd6Pt7/Loo+vy1xfVL0d2OFor7sF5K+4+6FY6cePyOL1mZ+t6T/6eYzS531fGkTMDjM7FRJn3H3g9KuZTBYNGL3OXf/aNq19JWZnSlpV3f/Stq1ACGhf0qfmR0l6SR3/6deGwNlrAL7m/mS3u3ufb0RS0Uzs+9J+qu7/yTtWrrT3d17APTAozuL3Zl2Hf3h7j9KuwYAQ69M+6ffKxpFAFBG3P3ytGsIibt/Ke0aesO0TwAAAAAYBpj2CQAAAADDACN/AAAAADAMEP4AAAAAYBgo+xu+ZLNZnzZtWtplABhEjz322Dp3z6Vdx46gbwIqTyX0TRL9E1CJ+to/lX34mzZtmpYuXZp2GQAGkZmtSruGHUXfBFSeSuibJPonoBL1tX9i2icAAAAADAOEPwAAAAAYBgh/AAAAADAMEP4AAAAAYBgg/AEAAADAMED4AwAAAIBhgPAHAAAAAMMA4Q8AAAAAhgHCHwAAAAAMA4Q/AAAAABgGCH8AAAAAMAwQ/gAAAABgGCD8AQAAAMAwQPgDAAAAgGGA8AcAAAAAwwDhD8PewoULB3+n114rTZsmN5OmTYseD6IhqRmp4c8TQKjon4DKQvhD2Riq/4AWLVo0KPtxd7W1tWnb1VfLTztNWrVKJkmrVsnnz1fHr341KO8jDV7NXfGffPI6Ojq0aNEitba2pl0KAGynubl5yP6/AZCOkWkXAPTVokWLeg0nHR0d2rRpk5qamrRx40Zt3Lix158lafbs2Wpra+v8am9v3+5xX57r6OiQJL0gaVqXuqy5WatOOkm7n3yyRo4c2aevqqqqbp+TpOOOO041NTWqqanR2LFje/ze03NjxoyRmfX5Mx6ohQsXEi5LuO222yRJjz/+uA488MCUqwGAN+y1115plwBgkBH+UBbuvfdeSdJXvvKVHgPd66+/LnfvcV9mplGjRqmlpaVz23333SdJ2n333bXnnnv2K4x13T71a18r+b5TJX31q1/tNlT2Fjb//Oc/68UXX+zc329+8xtJ0rhx4zRixAg1Nzervb19QJ9vIQhK0rve9S6NHz9+h79qamo6Q6U0tMGyHC1cuHC7M+qzZs2SJC1YsIDPCUCquvZPhb6c/gkof9bbgXLoZs6c6UuXLk27DAyR8847TxdccMGbttfV1Wn33XdXbW2tJk6cqIkTJ5b8udS2QlgqMLNeA2O/TJsmrVr15u277SYVhbcd0V3Nra2tam5u1pYtW7b7Xmpb4fttt92m+++//037mjp1qjKZjDZt2rTdV18/KzPTuHHjOsPgihUr+vPax9x9Zp8aB6qvfdPKlSu1xx576KqrrtKpp5469IUBGLBK6JukvvdP8+bNUz6fH9z/IwEMib72T4z8IVhLlizRTTfdJDPTl770JV188cXl8R/Qt74lzZ8vNTe/sa2mJto+xKqrqzvDbl+dc845nT/3FoTdXVu2bNHmzZvfFApLfd1555165JFHttu/xNnjYtlsVpJUX1+fciUAsL1C/wSgchD+EJy2tjZdeOGFWrRokd761rfq7rvv1uzZs3XxxRcPyfstWLBgcHd44onR9//4D/mqVbLddouCX2H7IBj0mvvIzDqvIczlcr227zptqCzCe8ImTJigESNGaN26dWmXAgDbyWazGj16dNplABhEhD8EZeXKlTrppJP00EMP6fjjj9cll1yiyZMnSxq6wDMkI1AnniideKKs95YDMlSjZmmFyuHMzLTzzjsz8gcgONlsVi0tLWptbVV1dXXa5QAYBCz1gCC4u6666irts88+Wr58ua699lpdd911ncFPYhmCJAzlZ0yw7F42myX8AQhOYdpnY2NjypUAGCyEP6SuoaFBxx13nD71qU/pfe97n5YtW6YTTjgh7bIwyAjv3cvlckz7BIYZM5tkZjeZ2XNm9qyZ/b2Z1ZlZ3sz+En+fHLc1M/uhma0ws2Vm9t6i/ZwSt/+LmZ0ymDUWwh/9E1A5CH9I1R133KH3vOc9+v3vf6+LLrpId999t6ZOnZp2WUCicrkcI3/A8PNfkm53970k7SPpWUnnSLrb3adLujt+LElHSJoef82XdKkkmVmdpAWSDpR0gKQFhcA4GDKZjCTCH1BJCH9IxZYtW3TWWWfp8MMP1+TJk/WnP/1JX/7yl1VVVZV2aUDimPYJDC9mVivpHyRdIUnuvs3dN0g6RtI1cbNrJH00/vkYSb/wyCOSJpnZzpI+JCnv7o3uvl5SXtLhg1UnI39A5SH8IXFPPPGEZs6cqR/+8Ic666yztHTpUu27775plwWkJpfLqampSdu2bUu7FADJ2F1SvaSrzOwJM/u5mY2TtJO7vyxJ8fe3xO13kbS66PVr4m3dbR8UhD+g8hD+kJj29nZddNFFOvDAA7V+/Xrdcccd+sEPfqCxY8emXRqQqsKyGRxgAcPGSEnvlXSpu+8nabPemOJZSqmbR3sP29+8A7P5ZrbUzJb2daZBYdpnQ0NDn9oDCB/hD4lYtWqV5syZo7PPPltHHXWUnnrqKc2bNy/tsoAgsNA7MOyskbTG3R+NH9+kKAy+Gk/nVPz9taL2uxa9foqktT1sfxN3v9zdZ7r7zL6s0ypJY8aM0fjx4zkxBVQQwh+G3LXXXqsZM2boscce01VXXaWbbrqp82wiAEb+gOHG3V+RtNrM9ow3zZH0jKTFkgp37DxF0u/inxdLOjm+6+csSRvjaaF3SJpnZpPjG73Mi7cNmmw2S98EVJBEF3k3s8MV3d2qStLP3f3CLs+fKum7kl6KN/3Y3X+eZI0YPOvXr9fnPvc53XDDDXr/+9+vX/7yl9p9993TLgsITiH8MfIHDCtnSrrWzEZJWinpk4pOyt9oZp+W9DdJH4vb3irpw5JWSGqO28rdG83sG5KWxO3Od/dBXZSP8AdUlsTCn5lVSbpE0lxF0xSWmNlid3+mS9Nfu/sZSdWFwbdw4UIdcsghOvnkk/XKK6/om9/8ps4++2yNHJnouQagbDDtExh+3P1JSTNLPDWnRFuX9Plu9nOlpCsHt7o3EP6AypLk0fgBkla4+0pJMrMbFN26uGv4QxlraWnRokWLZGaaPn26HnroIe2///5plwUELZPJyMw4wAIQnEwmo+effz7tMgAMkiSv+evr7Yj/0cyWmdlNZrZriecRoJaWFl1//fWdQe/000/X448/TvAD+qCqqkp1dXWM/AEIDiN/QGVJMvz15XbEv5c0zd1nSLpLbyx0uv2OBnC7YgyNVatW6bzzzlNdXZ1OOOEEPfXUU5Kkyy67TOPHj9fChQvTLRAoEyz0DiBE2WxWr7/+ulpaWtIuBcAgSDL89Xo7YndvcPdC7/IzSe8rtaOB3K4Yg6ejo0O33367jj76aO2+++76zne+o7lz5+qOO+5Qe3u7JMnd5e6EP6CPcrkcZ9cBBKdwTTJr/QGVIclr/pZImm5mb1d0N8+PSzqhuIGZ7RzfuliSjpb0bIL1oRfr1q3TVVddpcsuu0wrV67UTjvtpPPOO0+nnXaapk6dmnZ5QFnL5XL685//nHYZALCd4vD3tre9LeVqAOyoxMKfu7eZ2RmK1p+pknSluy83s/MlLXX3xZK+YGZHS2qT1Cjp1KTqQ2nurkcffVQ/+clPdOONN6qlpUWHHHKIvv3tb+vYY4/VqFGj3vSaBQsWpFApUN6y2aweeuihtMsAgO0Uwh8zE4DKkOi99939VkVr1RRv+3rRz+dKOjfJmlDa5s2bdd111+nSSy/VE088oQkTJugzn/mM/uVf/kV77713j69lqifQf4Vpnx0dHRoxIskZ+QDQPcIfUFlYeA3bee6553TppZfqmmuu0caNGzVjxgxddtllOvHEEzV+/Pi0ywMqVi6XU3t7uzZs2KC6urq0ywEASYQ/oNJwenmYKh6da21t1U033aTDDjtM73rXu3TZZZfpyCOP1IMPPqgnn3xSp59+OsEPGGIcYAEIUSaTkUTfBFQKwt8wtWjRIq1Zs0YLFizQbrvtpo997GN64YUXdOGFF2rNmjX61a9+pfe///0yK7VCB1A5zGxXM7vXzJ41s+VmdlaJNmZmPzSzFfE6pO8d7DoKdy5muQcAIamurlZtbS3hD6gQTPschh599FFJ0rRp09TR0aEjjjhCP/vZz3T44Yerqqoq5eqAxLVJ+pK7P25mEyQ9ZmZ5d3+mqM0RkqbHXwdKujT+PmgIfwBCxULvQOVg5G8Y+cIXviAz06xZsyRJ7e3tcnftv//++shHPkLww7Dk7i+7++Pxz68rWmJmly7NjpH0C488ImmSme08mHUw7RNAqAh/QOUg/A0DDQ0N+uIXv6jLLrtMNTU1nUsxsBA7sD0zmyZpP0mPdnlqF0mrix6v0ZsD4g5h5A9AqLLZLIu8AxWC8FfBtm7dqosuukh77LGHfvSjH+nUU0/VihUrCHtACWY2XtJvJH3R3Zu6Pl3iJV5iH/PNbKmZLe1viBs7dqzGjRtH+AMQHEb+gMpB+KtAHR0duvbaa7Xnnnvq7LPP1kEHHaRly5bp8ssv1847RzPVWIgdeIOZVSsKfte6+80lmqyRtGvR4ymS1nZt5O6Xu/tMd59ZGMnrDw6wAISIvgmoHIS/CnPvvfdq//331yc+8Qlls1ndc889uuWWW960MDujf0DEolvaXiHpWXf/z26aLZZ0cnzXz1mSNrr7y4NdSy6XY+QPQHCy2aw2b96sLVu2pF0KgB1E+KsQzzzzjI488kgddthhWrdunX71q19pyZIlmj17dtqlAaH7gKSTJB1mZk/GXx82s8+a2WfjNrdKWilphaSfSfrcUBRC+AMQosINqbjuDyh/LPVQ5l5++WUtWLBAV1xxhSZMmKCLLrpIZ555psaMGZN2aUBZcPc/qvQ1fcVtXNLnh7qWbDarZ555pveGAJCg4oXep0yZknI1AHYE4a9Mbdq0SRdffLEuvvhibdu2TWeeeaa+9rWvdXbQAMoPI38AQsRSNEDlIPyVmba2Nl155ZVasGCBXnnlFX3sYx/TBRdcoD322CPt0gDsoFwup+bmZjU3N6umpibtcgBAEuEPqCRc81cm3F233HKL9tlnH51++unaY4899PDDD+vGG28k+AEVggMsACHimj+gchD+Ardw4UI99thjOuyww3TUUUeptbVVN998sx544AHNmjUr7fIADCIWegcQorq6OkmcmAIqAeEvYKtWrdKiRYs0c+ZMPf300/rxj3+s5cuX69hjj1V0d3oAlYTwByBEI0eO1OTJkwl/QAXgmr9AdXR06OCDD5YknXvuuTr77LM1ceLElKsCMJSY9gkgVCz0DlQGRv4CtHDhQlVVVWn16tWSpAsuuECTJk1iYXagwjHyByBUhD+gMjDyF6CFCxdqzJgxOvfccyVFN3sBUPkmTZqkqqoqwh+A4GQyGa1ZsybtMgDsIEb+ApXP5zVjxoy0ywCQIDPj7DqAINE3AZWBkb8ANTc3649//KPOPPNMHXvssWmXAyBBLPQOIESF8Ofu3HQOKGOEvwA98MAD2rZtm+bOnasPfehDaZcDIEGEPwAhymaz2rp1q5qbmzVu3Li0ywEwQEz7DFA+n9eoUaM67/YJYPhgahWAELHQO1AZCH8ByufzOuigg1RTU5N2KQASxsgfgBCxFA1QGQh/gXn11Ve1bNkyzZ07N+1SAKQgl8upsbFRbW1taZcCAJ0If0BlIPwF5q677pIkwh8wTBUOsBobG1OuBADeQPgDKgPhLzD5fF6ZTEb77bdf2qUASAELvQMIUSaTkUT4A8od4S8g7q58Pq85c+ZoxAj+aIDhiPAHIESTJ0+WmRH+gDJHwgjIs88+q7Vr1zLlExjGmFoFIERVVVWqq6ujbwLKHOEvIPl8XhLX+wHDGSN/AELFUjRA+SP8BSSfz2v69Onabbfd0i4FQEoKI3+EPwChyWazrPMHlDnCXyC2bdum++67j1E/YJirrq7WxIkTObsOIDiM/AHlj/AXiEceeUSbN28m/AFgoXcAQSL8AeWP8BeIfD6vqqoqzZ49O+1SAKSM8AcgRIXw5+5plwJggAh/gcjn8zrggAM0ceLEtEsBkDLOrgMIUSaT0bZt27Rp06a0SwEwQIS/AKxfv15LlixhyicASYz8AQgTS9EA5Y/wF4B7771XHR0dhD8AkqIDrPr6eqZWAQgK4Q8of4S/AOTzeU2YMEEHHnhupbeMAAAgAElEQVRg2qUACEAul1Nra6tef/31tEsBgE6EP6D8Ef4CkM/ndeihh6q6ujrtUgAEgIXeAYSoEP5Y6w8oX4S/lL3wwgv661//ypRPAJ04uw4gRPRNQPkj/KUsn89LEuEPQCdG/gCEaOLEiaqqqiL8AWWM8JeyfD6vKVOmaM8990y7FACBIPwBCNGIESOUyWQIf0AZI/ylqL29XXfffbfmzp0rM0u7HACBYGoVgFAR/oDyRvhL0eOPP67169cz5RPAdsaPH6/Ro0cz8gcgONlslvAHlDHCX4oK1/vNmTMn5UoAhMTMWOgdQJAIf0B5I/ylKJ/Pa99999Vb3vKWtEsBEBgOsACEiL4JKG+Ev5Rs3rxZDz74IFM+AZTEyB+AEGWzWTU0NMjd0y4FwAAQ/lJy//33q7W1lfAHoCTCH4AQZbNZtbW1qampKe1SAAwA4S8l+Xxeo0eP1kEHHZR2KQACxNQqACHibsRAeSP8pSSfz+vggw/W2LFj0y4FQIByuZyamprU0tKSdikA0InwB5Q3wl8KXn75ZT399NNM+QTQrcJC7xxgAQgJ4Q8ob4S/FNx1112SRPgD0C0OsIDKZmYvmtlTZvakmS2Nty00s5fibU+a2YeL2p9rZivM7Hkz+1DR9sPjbSvM7JyhrjuTyUiibwLK1ci0CxiO8vm8crmc9tlnn7RLARCowsgfN30BKtpsd++aor7v7hcXbzCzd0v6uKS9Jb1N0l1m9s746UskzZW0RtISM1vs7s8MVcGcmALKG+EvYe6uu+66S3PmzNGIEQy8AiiN8AegyDGSbnD3FkkvmNkKSQfEz61w95WSZGY3xG2HLPzV1tZq5MiRhD+gTJE+ErZ8+XK9/PLLTPkE0CPOrgMVzyXdaWaPmdn8ou1nmNkyM7vSzCbH23aRtLqozZp4W3fb38TM5pvZUjNbuiMnlcyMuxEDZYzwl7B8Pi+J6/0A9Kyurk5mxsgfULk+4O7vlXSEpM+b2T9IulTSHpL2lfSypO/Fba3E672H7W/e6H65u89095mFmQUDVVjoHUD5IfwlLJ/Pa88999Suu+6adikAAlZVVaVMJkP4AyqUu6+Nv78m6beSDnD3V9293d07JP1Mb0ztXCOp+MBhiqS1PWwfUoz8AeWL8JeglpYW/eEPf2DUD0CfcIAFVCYzG2dmEwo/S5on6Wkz27mo2bGSno5/Xizp42Y22szeLmm6pD9JWiJpupm93cxGKbopzOKhrp++CShf3PAlQQ8//LCam5sJfwD6JJfLMfIHVKadJP3WzKToWOw6d7/dzH5pZvsqmrr5oqTTJcndl5vZjYpu5NIm6fPu3i5JZnaGpDskVUm60t2XD3XxhD+gfBH+EpTP51VVVaVDDz007VIAlIFcLqfnnnsu7TIADLL47pxvWu/J3U/q4TXfkvStEttvlXTroBbYi0wmo4aGBnV0dHDncqDM8C82Qfl8XrNmzVJtbW3apQAoA5xdBxCibDarjo4ObdiwIe1SAPQT4S8hjY2NWrp0KVM+AfRZLpfrPLsOAKFgKRqgfBH+EnLPPffI3Ql/APosl8upvb1d69evT7sUAOhE+APKF+EvIfl8XrW1tTrggAN6bwwA4gALQJgKfRNr/QHlh/CXkHw+r9mzZ2vkSO6xA6BvCgsxc8dPACHhxBRQvgh/CfjrX/+qF154gSmfAPqF8AcgRIQ/oHwR/hKQz+clifAHoF84wAIQonHjxmn06NH0TUAZSjT8mdnhZva8ma0ws3N6aHecmbmZzUyyvqGSz+c1depUTZ8+Pe1SAJQRRv4AhMjMlMlkCH9AGUos/JlZlaRLJB0h6d2Sjjezd5doN0HSFyQ9mlRtQ6m9vV333HOP5s6dKzNLuxwAZWTMmDEaP3484Q9AcFiHFChPSY78HSBphbuvdPdtkm6QdEyJdt+QdJGkrQnWNmSWLl2qDRs2MOUTwIBwgAUgRPRNQHlKMvztIml10eM18bZOZrafpF3d/ZaedmRm881sqZktDf2MeOF6vzlz5qRcCYBylMvlGPkDEBzCH1Cekgx/peY8eueTZiMkfV/Sl3rbkbtf7u4z3X1m4ZqYUOXzee23336dN24AgP4g/AEIUTabZZ0/oAwlGf7WSNq16PEUSWuLHk+Q9HeS7jOzFyXNkrS4nG/6smnTJj388MNM+QQwYJxdBxCibDarxsZGtbe3p10KgH5IMvwtkTTdzN5uZqMkfVzS4sKT7r7R3bPuPs3dp0l6RNLR7r40wRoH1R/+8Ae1trYS/gAMGCN/AEKUzWbl7lq/fn3apQDoh8TCn7u3STpD0h2SnpV0o7svN7PzzezopOpIUj6f15gxY3TQQQelXQqAMpXL5bRlyxZt3rw57VIAoBPrkALlaWSSb+but0q6tcu2r3fT9tAkahpK+XxeBx98sMaMGZN2KQDKVPEB1rhx41KuBgAimUxGEuEPKDeJLvI+nLz00kt65plnmPIJlAEzu9LMXjOzp7t5fqKZ/d7M/s/MlpvZJ5OqjYXeAYSIkT+gPBH+hshdd90lSYQ/oDxcLenwHp7/vKRn3H0fSYdK+l587fKQI/wBCBHhDyhPhL8hks/nlcvlNGPGjLRLAdALd79fUmNPTSRNMDOTND5u25ZEbRxgAQgRfRNQnhK95m+4cHfddddd+uAHP6gRI8jXQAX4saK7E69VtCzNP7t7RxJvzMgfgBDV1NRo7NixrPUHlBmSyRB46qmn9OqrrzLlE6gcH5L0pKS3SdpX0o/NrLZrIzObb2ZLzWzpYIW1iRMnauTIkZxdBxAc1iEFyg/hbwjk83lJXO8HVJBPSrrZIyskvSBpr66N3P1yd5/p7jMLI3Y7ysyUzWYZ+QMQHMIfUH4If0Mgn89rr7320pQpU9IuBcDg+JukOZJkZjtJ2lPSyqTenIXeAYSI8AeUH8LfINu6davuv/9+Rv2AMmJm10t6WNKeZrbGzD5tZp81s8/GTb4h6f1m9pSkuyWd7e6JHfHkcjkOsAAEJ5PJ0DcBZYYbvgyyhx56SFu2bCH8AWXE3Y/v5fm1kuYlVM6bZLNZPfnkk2m9PQCUxMgfUH4Y+Rtk+XxeI0eO1KGHHpp2KQAqBNM+AYQom81qw4YNam1tTbsUAH1E+Btk+Xxes2bN0oQJE9IuBUCFyOVyWr9+vdraEllaEAD6pLDWX2NjT8ukAggJ4W8QrVu3To8//jhTPgEMqsIBFutpAQgJC70D5YfwN4juvvtuuTvhD8CgYqF3ACHixBRQfgh/gyifz2vixInaf//90y4FQAUphD/OrgMICSN/QPkh/A0Sd1c+n9fs2bM1ciQ3UQUweAoHWIz8AQgJ4Q8oP4S/QfKXv/xFf/vb35jyCWDQMe0TQIgymYwkwh9QTgh/gySfz0sS4Q/AoOMAC0CIxowZo3HjxtE3AWWE8DdI8vm8dtttN73jHe9IuxQAFaa6ulqTJk1i5A9AcFjoHSgvhL9B0NbWpnvvvVdz586VmaVdDoAKxELvAEJE+APKC+FvEPzpT39SU1MTUz4BDBkOsACEiL4JKC+Ev0GQz+dlZpozZ07apQCoUIz8AQhRNptlnT+gjBD+BsEVV1yh9773vZ03ZQCAwUb4AxAiRv6A8kL420FNTU1avXo1Uz4BDKnCAZa7p10KAHTKZrNqamrStm3b0i4FQB8Q/nbQww8/LEn64Ac/mHIlACpZLpdTa2urmpqa0i4FADoVFnpn6idQHgh/A7Rw4UKZmQ4//HBJUfgzMy1cuDDdwgBUJBZ6BxAi1iEFygvhb4AWLlwod9eFF14oSdq0aZPcnfAHYEgUzq5zgAUgJPRNQHkh/O2gwln4mpqalCsBUMkY+QMQIsIfUF4Ifzuovr5etbW1LO4OYEgR/gCEiPAHlBfC3w6qr6/XO97xjrTLAFDhOMACEKLCNX/c8AUoD4S/HVRfX995Rh4Ahsq4ceM0ZswYRv4ABGXUqFGqra3lxBRQJgh/O4jwByAJZsZC7wCCxELvQPkg/O0gwh+ApHCABSBE9E1A+SD87YDm5mY1NzcT/gAkgpE/ACHKZDKEP6BMEP52QOEgjPAHIAmEPwAhYuQPKB+Evx1A+AOQJA6wAISIvgkoH4S/HUD4A5CkXC6n119/XS0tLWmXAgCdstmsNm/erK1bt6ZdCoBeEP52AOEPQJJY6B1AiArrkLLWHxA+wt8OIPwBSBILvQMIEX0TUD4Ifzugvr5e1dXVmjhxYtqlABgGGPkDECLCH1A+CH87oL6+XtlsVmaWdikAhgHCH4AQEf6A8kH42wEs8A4gSRxgAQhRJpORRN8ElAPC3w4g/AFIUl1dnUaMGMHIH1ABzOxFM3vKzJ40s6Xxtjozy5vZX+Lvk+PtZmY/NLMVZrbMzN5btJ9T4vZ/MbNT0vhd6urqJBH+gHJA+NsBhD8ASRoxYoQymQzhD6gcs919X3efGT8+R9Ld7j5d0t3xY0k6QtL0+Gu+pEulKCxKWiDpQEkHSFpQCIxJqq6u1qRJkwh/QBkg/O0Awh+ApLGYMlDRjpF0TfzzNZI+WrT9Fx55RNIkM9tZ0ock5d290d3XS8pLOjzpoiX6JqBcEP4GqKWlRU1NTYQ/AInK5XKM/AGVwSXdaWaPmdn8eNtO7v6yJMXf3xJv30XS6qLXrom3dbc9cYQ/oDyMTLuAclXo4Ah/AJKUy+X07LPPpl0GgB33AXdfa2ZvkZQ3s+d6aFvqtuLew/Y37yAKmPMlaerUqf2ttVfZbFYvvfTSoO8XwOBi5G+ACH8A0pDNZhn5AyqAu6+Nv78m6beKrtl7NZ7Oqfj7a3HzNZJ2LXr5FElre9he6v0ud/eZ7j5zKI5dGPkDygPhb4AKB1+EPwBJyuVyamhoUEdHR9qlABggMxtnZhMKP0uaJ+lpSYslFe7YeYqk38U/L5Z0cnzXz1mSNsbTQu+QNM/MJsc3epkXb0sc4Q8oD0z7HCDCH4A05HI5dXR0aP369Z1rawEoOztJ+q2ZSdGx2HXufruZLZF0o5l9WtLfJH0sbn+rpA9LWiGpWdInJcndG83sG5KWxO3Od/fG5H6NN2QyGW3ZskXNzc2qqalJowQAfUD4GyDCH4A0FBZ6r6+vJ/wBZcrdV0rap8T2BklzSmx3SZ/vZl9XSrpysGvsr0LftG7duiG5phDA4GDa5wDV19fLzDR5cuLL6QAYxgonnLjuD0BIisMfgHAR/gaocNa9qqoq7VIADCOF8McBFoCQEP6A8kD4GyAWeAeQhuJpnwAQCsIfUB4IfwNE+AOQBqZ9AghRIfw1NDSkXAmAnhD+BojwByANo0eP1oQJEzi7DiAokydPlpnRNwGBI/wNEOEPQFpY6B1AaKqqqlRXV0f4AwJH+BuA9vZ2NTY2Ev4ApCKXyxH+AAQnk8kQ/oDAEf4GoKGhQe5O+AOQilwuxwEWgOBks1n6JiBwhL8BYIF3AGli2ieAEBH+gPAR/gaA8AcgTYVpn+6edikA0InwB4SP8DcAhD8Aacrlctq6dauam5vTLgUAOhXCHyemgHAR/gaA8AcgTSz0DiBE2WxW27Zt0+bNm9MuBUA3CH8DUDjgKhyAAUCSWOgdQIgKx0VM/QTCRfgbgPr6ek2aNEnV1dVplwJgGCqEPw6wAISE8AeEj/A3ACzwDiBNTPsEECLCHxA+wt8AEP4ApIlpnwBClMlkJBH+gJAlGv7M7HAze97MVpjZOSWe/6yZPWVmT5rZH83s3UnW11eEPwBpqq2tVXV1NQdYAILCyB8QvsTCn5lVSbpE0hGS3i3p+BLh7jp3f4+77yvpIkn/mVR9/bFu3TrCH4DUmBkLvQMIzqRJkzRixAjCHxCwJEf+DpC0wt1Xuvs2STdIOqa4gbs3FT0cJym4hWLcnfAHIHWFhd4BIBQjRoxQJpMh/AEBG5nge+0iaXXR4zWSDuzayMw+L+nfJI2SdFgypfXdhg0b1NbWRvgDkKpcLscBFoDgZLNZNTQ0pF0GgG4kOfJnJba9aWTP3S9x9z0knS3pqyV3ZDbfzJaa2dKkz3yzwDuAEDDtE0CIstksJ6aAgCUZ/tZI2rXo8RRJa3tof4Okj5Z6wt0vd/eZ7j4z6RBG+AMQAqZ9AggR4Q8IW5Lhb4mk6Wb2djMbJenjkhYXNzCz6UUPPyLpLwnW1yeEPwAhyOVy2rBhg1pbW9MuBQA6Ef6AsCV2zZ+7t5nZGZLukFQl6Up3X25m50ta6u6LJZ1hZh+U1CppvaRTkqqvrwrhr3A7YwBIQ6EPamho0Fvf+taUqwGASOGGL+4us1JX/ABIU5I3fJG73yrp1i7bvl7081lJ1jMQjPwBlcXMrpR0pKTX3P3vumlzqKQfSKqWtM7dD0muwtKKF3on/AEIRTabVVtbm5qamjRx4sS0ywHQRaKLvFeC+vp6jRs3TmPHjk27FACD42pJh3f3pJlNkvQTSUe7+96SPpZQXT1iMWUAIaJvAsJG+Oun+vp6Rv2ACuLu90tq7KHJCZJudve/xe1fS6SwXhSP/AFAKAh/QNgIf/1E+AOGnXdKmmxm95nZY2Z2ctoFSYQ/AGEi/AFhS/Sav0pQX1+vnXfeOe0yACRnpKT3SZojaaykh83sEXf/c9eGZjZf0nxJmjp16pAWVVdXJ4kDLABhKb4ZFYDwMPLXT4z8AcPOGkm3u/tmd18n6X5J+5RqmOQapNXV1Zo8eTIjfwCCwsgfEDbCXz+4O+EPGH5+J+lgMxtpZjWSDpT0bMo1SWKhdwDhqa2t1ciRIwl/QKCY9tkPmzZtUktLC+EPqCBmdr2kQyVlzWyNpAWKlnSQu1/m7s+a2e2SlknqkPRzd386rXqLsZgygNCYWedafwDCQ/jrB9b4AyqPux/fhzbflfTdBMrpl1wup5UrV6ZdBgBshxNTQLiY9tkPhD8AIWHaJ4AQEf6AcBH++oHwByAkhQMsd0+7FADoRPgDwkX46wfCH4CQ5HI5tbW1aePGjWmXAgCdCH9AuAh//UD4AxCSQl/EQRaAkGSzWTU2NqqjoyPtUgB0Qfjrh/r6eo0ePVrjx49PuxQA6FxPi+v+AIQkm82qvb2dWQlAgAh//bBu3TrlcjmZWdqlAEDnyB/hD0BIWOgdCBfhrx9Y4B1ASJj2CSBEmUxGEn0TECLCXz8Q/gCEhGmfAELEyB8QLsJfPxD+AIRk3LhxGjt2LOEPQFAIf0C4CH/9QPgDEJpcLscBFoCgEP6AcBH++mjr1q3atGkT4Q8IjJmNMbM3/cM0s7eY2Zg0akpSNptl5A9IwXDve3oyfvx4jRo1ivAHBIjw10es8QcE64eSDi6xfa6k7ydcS+JyuRzhD0jHsO57emJmymazamhoSLsUAF30Gv7MbK6Z/czM9o0fzx/6ssJD+AOCdZC739x1o7tfK+kfUqgnUUz7BFIzrPue3mSzWfomIEAj+9Dmc5I+KemrZlYnad+hLSlMhfBXmMcOIBg9LbxZ8bMbmPYJpGZY9z29IfwBYepL51Tv7hvc/d8lzZO0/xDXFCRG/oBgvWZmB3TdaGb7S6r4VJTL5bRp0yZt3bo17VKA4WZY9z29yWQyhD8gQH0Z+fvfwg/ufo6ZnTmE9QSL8AcE68uSbjSzqyU9Fm+bKelkSR9Pq6ikFC/0PmXKlJSrAYaVYd339IaRPyBMvY78ufvvujz+0dCVE676+npVVVVp0qRJaZcCoIi7/0nSgYqmYJ0af5mkA9390fQqSwYLvQPpGO59T2+y2awaGxvV3t6edikAivRl5E9mdpKk/5TUIuk8d/+Fmc2SdKSkI9z9fUNYYxDq6+uVzWY1YsSwn8YPBMfdX5W0IO060lAY+SP8Ackbzn1Pb7LZrNxd69ev534JQED6FP4kfV3ShyW9IOkMM8tL2kvS9ZK+OES1BYUF3oEwmdlTkrzUU5I63H2fhEtKVPG0TwDJ6aXvcXefkXBJQSle6J3wB4Sjr+Fvk7svkSQzWyTpVUnvdPcNQ1ZZYAh/QLCOLLHNJE2RdF7CtSSOaZ9Aakr1PYgVhz8A4ehr+HtrvL7f8/HXmuEU/KTowGq//fZLuwwAXbj7qsLP8XqkJ0j6J0UzFX6TVl1JmTx5sqqqqgh/QMKK+54CM8tKanD3UiOCw0oh/LHQOxCWvoa/BZJmSDpR0nskTTCzuyQ9IekJd79uiOoLBiN/QJjM7J2K7qx3vKQGSb+WZO4+O9XCEjJixAhuqQ6kIL73wYWSGiV9Q9IvJWUljTCzk9399jTrSxsjf0CY+hT+3P3y4sdmNkVRGHyPpCMkVXT4a21t1YYNGwh/QJiek/SApKPcfYUkmdm/pltSsljoHUjFjxVNLZ8o6R5FN8B7xMwK90QY1uEvk8lIIvwBoenryN923H2NpDWSbh3ccsJU6LgIf0CQ/lHRyN+9Zna7pBsUXfM3bORyOcIfkLyR7n6nJJnZ+e7+iCS5+3Nmw6oLKqmmpkZjxowh/AGBYd2CPmCBdyBc7v5bd/9nRXcgvk/Sv0raycwuNbN5qRaXkFwuxwEWkLyOop+3dHlu2F/zZ2Ys9A4EiPDXB4z8AeFz983ufq27H6noTp9PSjon5bISwbRPIBX7mFmTmb0uaUb8c+Hxe/qyAzOrMrMnzOyW+PHVZvaCmT0Zf+0bbzcz+6GZrTCzZWb23qJ9nGJmf4m/ThmKX3SgCH9AeAY07XO4YeQPKC/u3ijpp/FXxcvlcmpoaFB7e7uqqqrSLgcYFtx9MP6xnSXpWUm1Rdu+7O43dWl3hKTp8deBki6VdKCZ1Sm6Kd9MRaONj5nZYndfPwi17TDCHxAeRv76gPAHIGS5XE7urvXrgzjeA9AH8c3zPiLp531ofoykX3jkEUmTzGxnSR+SlHf3xjjw5SUdPmRF9xPhDwgP4a8P6uvrZWadd64CgJCw0DtQln4g6Sva/tpBSfpWPLXz+2Y2Ot62i6TVRW3WxNu62x6EbDbLOn9AYAh/fVBfX6+6ujqmUwEIUmFWAuEPKA9mdqSk19z9sS5Pnavo5lX7S6qTdHbhJSV24z1sL/We881sqZktTaqvyGazWr9+vdra2hJ5PwC9I/z1AQu8AwhZoX9iehVQNj4g6Wgze1HR8jSHmdmv3P3leGpni6SrJB0Qt18jadei10+RtLaH7W/i7pe7+0x3n5nUMU1hVkJjY2Mi7wegd4S/PiD8AQgZ0z6B8uLu57r7FHefpmid0nvc/RPxdXyyaKHAj0p6On7JYkknx3f9nCVpo7u/LOkOSfPMbLKZTZY0L94WBBZ6B8LD3T77oL6+XnvttVfaZQBASYQ/oGJca2Y5RdM5n5T02Xj7rZI+LGmFpGZJn5SiOxub2TckLYnbnR/f7TgIhb6J8AeEg/DXB/X19Tr44IPTLgMASho9erRqa2s5wALKkLvfJ+m++OfDumnjkj7fzXNXSrpyiMrbIYQ/IDxM++xFR0eHGhoamPYJIGgs9A4gNIQ/IDyEv140Njaqo6OD8AcgaLlcjvAHIChc8weEh/DXCxZ4B1AOcrkcB1gAgjJ27FiNGzeOtf6AgBD+elEIf4WpCwAQIqZ9AghRNpvlxBQQEMJfLxj5A1AOCtM+o/tCAEAYCH9AWAh/vSD8ASgHuVxOLS0t2rx5c9qlAECnTCZD+AMCQvjrBdM+AZQD1voDECJG/oCwEP56UV9fr9raWo0ePTrtUgCgW4XZCYQ/ACEh/AFhIfz1or6+nimfAIJX6Kc4yAIQkmw2q6amJm3bti3tUgCI8Ncrwh+AcsC0TwAhKvRNLPcAhIHw1wvCH4BywMgfgBAR/oCwEP56QfgDUA4mTJigUaNGMfIHICiF8MeJKSAMhL8euLvWrVtH+AMQPDNjoXcAwSH8AWEh/PWgqalJra2thD8AZSGXy3GABSAomUxGEuEPCAXhrwcs8A6gnDDyByA0hD8gLIS/HhD+AJSTXC5H+AMQlNGjR2vChAmEPyAQhL8eEP4AlBOmfQIIEQu9A+Eg/PWA8AegnGSzWW3YsEGtra1plwIAnQh/QDgIfz0g/AEoJ6z1ByBEhD8gHIS/HtTX16umpkY1NTVplwIAvSL8AQhRNptlkXcgEIS/HrDAO4ByUlhPi5u+AAgJI39AOAh/PSD8ASgnhf6K8AcgJJlMRps2bdLWrVvTLgUY9gh/PSD8ASgnTPsEEKLCrASmfgLpI/z1gPAHoJzU1dVJYuQPQFgK4Y8TU0D6CH89IPwBKCcjR45UXV0d4Q9AUAh/QDgSDX9mdriZPW9mK8zsnBLP/5uZPWNmy8zsbjPbLcn6im3evFlbtmwh/AEoKyz0DiA0hD8gHImFPzOrknSJpCMkvVvS8Wb27i7NnpA0091nSLpJ0kVJ1ddV4cx5ocMCgHKQzWYZ+QMQFMIfEI4kR/4OkLTC3Ve6+zZJN0g6priBu9/r7s3xw0ckTUmwvu2wwDuAcpTL5Qh/AIJSuB6ZG74A6Usy/O0iaXXR4zXxtu58WtJtQ1pRDwh/AMoR0z4BhKa6ulqTJk2ibwICMDLB97IS27xkQ7NPSJop6ZBunp8vab4kTZ06dbDq2w7hD0A5Kiym7O4yK9XtAkDyMpkM4Q8IQJIjf2sk7Vr0eIqktV0bmdkHJf2HpKPdvaXUjtz9cnef6e4zhyqcEf4AlKNcLqe2tjZt2LAh7VIAoFPhxBSAdCUZ/pZImm5mbzezUZI+LmlxcQMz20/STxUFv9cSrO1N6uvrVV1drdra2jTLAJAAM7vSzF4zs6d7abe/mbWb2XFJ1dZfLPQOIESEPyAMiYU/d2+TdIakOyQ9K+lGd19uZueb2dFxs+9KGi/pv83sSTNb3M3uhlxhjT+mTQHDwtWSDu+pQXzH4u8o6sOCVbirHjd9ARASwh8QhiSv+ZO73yrp1i7bvl708weTrKcn69atY8onMEy4+/1mNq2XZmdK+o2k/Ye8oB1Q6LcIfwBCQvgDwpDoIu/lpDDyBwBmtqInemMAACAASURBVIukYyVdlnYtvWHaJ4AQZbNZbdmyRc3Nzb03BjBkCH/dIPwBKPIDSWe7e3tPjcxsvpktNbOlaY28Me0TQIgKfRNr/QHpIvx1g/AHoMhMSTeY2YuSjpP0EzP7aNdGSdyJuDc1NTWqqakh/AEISiH8MSsBSFei1/yVi5aWFjU1NRH+AEiS3P3thZ/N7GpJt7j7/6RXUc9Y6B1AaDKZjCTCH5A2wl8JhY6J8AcMD2Z2vaRDJWXNbI2kBZKqJcndg7/Or6tsNsvIH4CgMPIHhIHwVwILvAPDi7sf34+2pw5hKf+/vTuPr6K+9z/++mQDjJhAEsACAoJLAUEWWa8gcKZga7Xeq5ZWLVZburhRqte9SRRrrbVK708Ud6wLUmuvtlosENBYUFmEsLigIBrZwiqyBJJ8f39kkhtCgADJmXPOvJ+PRx6cM2eSvMMDvjnvme/Mt0Hk5OSo/IlITFH5E4kNuuavDip/IhLPNO1TRGJNixYtMDONTSIBU/mrg8qfiMQzTfsUkViTkpJCixYtVP5EAqbyVweVPxGJZzk5OezcuZPdu3cHHUVEpJoWehcJnspfHUpKSkhOTqZFixZBRxEROWJa6F1EYlF2drbW+RMJmMpfHUpKSsjKyiIpSX89IhJ/tNC7iMQinfkTCZ7aTR20wLuIxLOq8UvlT0RiicqfSPBU/uqg8ici8UzTPkUkFmVlZbFp0yacc0FHEQktlb86qPyJSDzTtE8RiUXZ2dmUlpayc+fOoKOIhJbKXx1U/kQknmVmZpKcnKzyJyIxRQu9iwRP5a+WsrIytmzZovInInErKSlJ19aISMxR+RMJnspfLVW3IK4aoERE4pEWeheRWKPyJxI8lb9atMC7iCSCnJwcvcESkZii8icSPJW/WlT+RCQR5OTk6MyfiMSUqvKnhd5FgqPyV4vKn4gkAk37FJFYk5mZSVJSks78iQRI5a+WqgFJ5U9E4llOTg5btmyhvLw86CgiIkDlzahatmyp8icSIJW/WqqOlGdlZQWcRETk6OXk5OCcY8uWLUFHERGppjsRiwRL5a+WkpISWrRoQWpqatBRRESOmhZ6F5FYpPInEiyVv1q0wLuIJIKqcUxvskQklqj8iQRL5a8WlT8RSQRV45jO/InELjNLNrP3zewf/vNOZvauma00sxfNLM3f3sR//on/escaX+MWf/tHZjYymJ+k/lT+RIKl8leLyp+IJIKqaZ9PPfVUwElE5BCuBz6o8fxe4AHn3CnAVuAqf/tVwFbnXBfgAX8/zKwrMBroBowCJplZcpSyH5Wq8uecCzqKSCip/NWi8iciiaCq/L322msBJxGRuphZO+A7wOP+cwOGAy/5u0wBvuc/vsB/jv/6CH//C4CpzrlS59xq4BOgX3R+gqOTnZ1NWVkZO3bsCDqKSCip/NVQUVHBpk2bVP5EJO6lpaWRkZERdAwRObgHgf8GKvznWcA251yZ/7wYaOs/bgt8AeC/vt3fv3p7HZ8Tk6oOTGnqp0gwVP5q2LZtG+Xl5Sp/IhLX8vLyMDO2b98OgJlhZuTl5QUbTEQAMLPzgI3OuYU1N9exqzvMa4f6nNrfc6yZLTCzBUFeC1y1lJbKn0gwVP5qqBoMVf5EJJ7l5eXhnOPWW28FYPv27TjnVP5EYsdg4Hwz+wyYSuV0zweBTDNL8fdpB6z1HxcD7QH81zOALTW31/E5+3HOPeqc6+uc6xvk+xyd+RMJlspfDSp/IpJIIpEIAHPmzAk2iIjsxzl3i3OunXOuI5U3bClwzl0KzAYu8ncbA7ziP37Vf47/eoGrvGPKq8Bo/26gnYBTgPei9GMcFZU/kWCp/NWg8iciiWTQoEGkpKQwc+bMoKOISP3cBIw3s0+ovKbvCX/7E0CWv308cDOAc245MA1YAUwHrnbOlUc99RFQ+RMJVsrhdwkPlT8RSSRNmjQhEokwY8aMoKOIyEE45+YAc/zHq6jjbp3OuT3AxQf5/LuBuxsvYcPKyMggOTlZ5U8kIDrzV4PKn4gkGs/z+PDDDykuLg46iogIZqaF3kUCpPJXQ0lJCc2bN6dJkyZBRxERaRCe5wHo7J+IxIzs7Gw2b94cdAyRUFL5q0ELvItIounevTutW7dW+RORmKEzfyLBUfmrQeVPRBKNmRGJRJg5cyYVFRWH/wQRkUaWlZXFihUrgo4hEkoqfzWo/IlIIvI8j5KSEpYuXRp0FBERnfkTCZDKXw0qfyKSiKrW+9PUTxGJBVXLPZSXx/SqFCIJSeXP55yjpKSkekASEUkUbdu2pWvXrip/IhKovLw8zIzf/va3AKSkpGBm5OXlBRtMJERU/nxff/01e/fu1Zk/EUlInufx1ltvsWfPnqCjiEhI5eXl4Zxj/fr1APz2t7/FOafyJxJFKn8+rfEnIoksEomwZ88e/v3vfwcdRURCrnXr1gDMnDkz4CQi4aPy51P5E5FENnToUFJSUvRmS0RiwoABA3j77bfZtWtX0FFEQkXlz6fyJyKJrHnz5gwcOFDX/YlITMjNzWXv3r28/fbbQUcRCRWVP5/Kn4gkOs/zWLRoEZs3bw46ioiE3Nlnn01aWpoOSIlEmcqfT+VPRBKd53k455g1a1bQUUQk5NLT0xk0aJDKn0iUqfz5SkpKaNq0Kenp6UFHERFpFH379iUjI0NvtkQkJniex5IlS9i4cWPQUURCQ+XPV7XAu5kFHUVEpFGkpKQwbNgwZsyYgXMu6DgiEnKRSARAsxFEokjlz1dV/kREEpnneaxZs4ZPP/006CgiEnJ9+vQhMzNTdyEWiSKVP5/Kn4iEged5AJr6KSKBS05OZvjw4ZqNIBJFKn8+lT8RCYMuXbrQoUMHlT8RiQme5/HFF1+wcuXKoKOIhILKn0/lT0TCwMzwPI+CggLKysqCjiMiIafZCCLRpfIH7N69m507d6r8iUgoRCIRtm/fzoIFC4KOIiIhd/LJJ9OxY0eVP5EoUflDa/yJSLiMGDECM9NNFkQkcFWzEWbPnq3ZCCJRoPKHyp+IhEt2dja9evXSkXYRiQmRSISvvvqK+fPnBx1FJOGp/KHyJyLh43ke8+bN4+uvvw46ioiE3PDhwzUbQSRKVP5Q+ROR8PE8j3379vHmm28GHUVEQk6zEUSiR+UPlT8RCZ/BgwfTtGlTvdkSkZhQNRthx44dQUcRSWgqf1SWv9TUVDIyMoKOIiISFU2bNuXss89W+RORmOB5HmVlZbz11ltBRxFJaCp/VJa/7OxszCzoKCIiUeN5HitWrGDt2rVBRxGRkNNsBJHoUPlDC7yLSDhVLa6smyyISNCqZiNoPBJpXCp/wKZNm1T+RCR0evToQU5Ojo60i0hMiEQiLF++XLMRRBqRyh//N+1TRCRMkpKSGDFiBDNnzsQ5F3QcEQm5qtkIs2bNCjiJSOJS+UPTPkUkvDzPY/369SxbtizoKCIScj179iQ7O1uzEUQaUejL3759+9i2bZvKn4iEkq77E5FYodkIIo0v9OVv06ZNgNb4E5Fwat++PaeddpqOtItITPA8j3Xr1rFixYqgo4gkpNCXPy3wLiJh53keb775JqWlpUFHEZGQi0QiADogJdJIolr+zGyUmX1kZp+Y2c11vD7EzBaZWZmZXRSNTCp/IhJ2nuexa9cu5s2bF3QUEQm5Dh06cMopp2gqukgjiVr5M7Nk4CHgXKAr8AMz61prt8+BK4Dno5VL5U9Ewm7o0KEkJyfrSLuIxATP85gzZw579+4NOopIwonmmb9+wCfOuVXOub3AVOCCmjs45z5zzhUBFdEKpfInImGXkZFB//79Vf5EJCZEIhF27tzJO++8E3QUkYQTzfLXFviixvNif9sRM7OxZrbAzBZUlbejVVJSgpnRsmXLY/o6IiLxzPM8FixYwNatW4OOIiIhN2zYMJKSkjT1U6QRRLP8WR3bjuo+vs65R51zfZ1zfY/1jF1JSQlZWVkkJycf09cREYlnnufhnKOgoCDoKCIScpmZmZx11lmajSDSCKJZ/oqB9jWetwPWRvH710kLvIuEm5k9aWYbzazOVc7N7FIzK/I/5ppZz2hnjIZ+/frRvHlzvdkSkZjgeR7vvfce27dvDzqKSEKJZvmbD5xiZp3MLA0YDbwaxe9fJ5U/kdB7Ghh1iNdXA0Odcz2Au4BHoxEq2lJTUznnnHNU/kQkJkQiESoqKpg9e3bQUUQSStTKn3OuDLgGeAP4AJjmnFtuZnea2fkAZnaWmRUDFwOTzWx5Y+dS+RMJN+fcW8CWQ7w+1zlXdSHcO1TOWkhInuexatUqVq1aFXQUEQm5gQMHkp6eruv+RBpYVNf5c8697pw71TnX2Tl3t7/tN865V/3H851z7Zxz6c65LOdct8bOpPInIkfgKuCfB3uxIW9GFQTP8wAtriwiwUtLS2Po0KEaj0QaWFTLX6wpLy9n8+bNKn8iclhmNozK8nfTwfZpyJtRBeG0006jXbt2OtIuIjEhEonw8ccf8/nnnwcdRSRhhLr8bdmyBeecyp+IHJKZ9QAeBy5wzm0OOk9jMTM8z2PWrFmUl5cHHUdEQq5qNoIOSIk0nFCXPy3wLiKHY2YnAS8DlzvnPg46T2OLRCJs3bqVRYsWBR1FREKuW7dutGnTRlM/RRqQyh8qfyJhZmYvAPOA08ys2MyuMrOfm9nP/V1+A2QBk8xssZktCCxsFEQiEUDX/YlI8MyMSCTCrFmzqKioCDqOSEJQ+UPlTyTMnHM/cM6d6JxL9W849YRz7hHn3CP+6z9xzrVwzp3pf/QNOnNjatWqFT179lT5E5GY4HkeJSUlFBUVBR1FJCGo/KHyJyJSk+d5zJ07l507dwYdRURCbsSIEYBmI4g0lFCXv02bNgGQnZ0dcBIRkdjheR579+6lsLAw6CgiEnJt27ala9euuumLSAMJdfkrKSkhMzOT1NTUoKOIiMSMs88+myZNmuhIu4jEhEgkwltvvcWePXuCjiIS90Jf/jTlU0Rkf82aNWPw4MEqfyISEzzPY8+ePcydOzfoKCJxT+VP5U9E5ACe57F06VLWr18fdBQRCbmhQ4eSkpKiA1IiDSD05U/X+4mIHEiLK4tIrGjevDkDBgzQeCTSAEJf/nTmT0TkQL169SIrK0tvtkQkJniex8KFC9m8eXPQUUTiWmjLn3OOTZs2qfyJiNQhKSmJESNGMGPGDJxzQccRkZCLRCI45ygoKAg6ikhcC2352759O/v27VP5ExE5iEgkwtq1a/nggw+CjiIiIdevXz+aN2+u2Qgixyi05U8LvIuIHFrVdX+6yYKIBC0lJYVhw4ZpPBI5Rip/Kn8iInXq2LEjXbp00ZstkQZmZk3N7D0zW2Jmy80s39/+tJmtNrPF/seZ/nYzsz+Z2SdmVmRmvWt8rTFmttL/GBPUzxQNnuexevVqVq1aFXQUkbil8qfyJyJyUJ7nMWfOHPbt2xd0FJFEUgoMd871BM4ERpnZAP+1G51zZ/ofi/1t5wKn+B9jgYcBzKwlkAv0B/oBuWbWIoo/R1RpNoLIsVP5U/kTETkoz/PYuXMn77zzTtBRRBKGq/S1/zTV/zjUnZUuAJ7xP+8dINPMTgRGAjOcc1ucc1uBGcCoxswepFNPPZV27dqp/IkcA5U/lT8RkYMaNmwYSUlJerMl0sDMLNnMFgMbqSxw7/ov3e1P7XzAzJr429oCX9T49GJ/28G2JyQzw/M8CgoKKC8vDzqOSFwKdflLT0+nWbNmQUcREYlZmZmZnHXWWSp/Ig3MOVfunDsTaAf0M7PuwC3A6cBZQEvgJn93q+tLHGL7AcxsrJktMLMFVQfA41EkEmHr1q0sWrQo6CgicSnU5U9n/UREDs/zPN577z22bdsWdBSRhOOc2wbMAUY559b5UztLgaeovI4PKs/ota/xae2AtYfYXtf3edQ519c51zee3/+MGDECQEs+iBwllT8RETkkz/OoqKhg9uzZQUcRSQhmlmNmmf7jZkAE+NC/jg8zM+B7wDL/U14FfuTf9XMAsN05tw54A/iWmbXwb/TyLX9bwmrdujU9evTQbASRo6TyJyIihzRgwADS09N1pF2k4ZwIzDazImA+ldf8/QN4zsyWAkuBbGCCv//rwCrgE+Ax4JcAzrktwF3+15gP3OlvS2ie5/Hvf/+bXbt2BR1FJO6kBB0gKCUlJZxxxhlBxxARiXlpaWmcc845OtIu0kCcc0VArzq2Dz/I/g64+iCvPQk82aABY5znedx///0UFhYycuTIoOOIxJVQnvlzzunMn4jIEYhEIqxcuZI1a9YEHUVEQu7ss88mLS1NB6REjkIoy9/OnTvZs2ePyp+ISD1pcWURiRXHHXccgwcP1lR0kaMQyvKnNf5ERI5M165d+cY3vqHyJyIxIRKJsGTJEjZs2BB0FJG4Esryt2nTJkDlT0SkvsyMSCTCrFmzqKioCDqOiIRc1WyEgoKCgJOIxJdQlj+d+RMROXKe57F582YWL14cdBQRCbnevXvTokULzUYQOUIqfyIiUi9ViyvrzZaIBC05OZnhw4czY8YMKm+GKiL1ofInIiL1cuKJJ9K9e3eVPxGJCZ7nUVxczMcffxx0FJG4Edry16RJE44//vigo4iIxBXP83j77bfZvXt30FFEJOQikQig2QgiRyK05S8nJwczCzqKiEhc8TyP0tJSxo4dG3QUEQm5zp0706lTJy35IHIEQlv+srOzg44hIhJ3hgwZQlpaGs8++2zQUUREiEQizJ49m7KysqCjiMSF0JY/Xe8nInLk0tPTGTRoEAC7du0KOI2IhJ3neXz11VfMnz8/6CgicUHlT0RE6iUvLw8zY86cOUBlETQzbr/99mCDiUhoDR8+HDPTdX8i9aTyJyIi9ZKXl4dzrvq26kOGDAFgypQpPProo+zbty/IeCISQllZWfTu3VvlT6SeQlf+SktL2bFjh8qfiMgxmjNnDjNmzKBdu3b87Gc/4/TTT2fKlCm69kZEosrzPN555x127NgRdBSRmBe68qc1/kREjl1ubi5mRiQSYe7cubz22mtkZGRwxRVX0L17d6ZOnUpFRUXQMUUkBCKRCGVlZbz55ptBRxGJeSp/IiJyxPLy8qofmxnf/va3WbhwIS+//DKpqan84Ac/oGfPnvztb3+rniYqItIYBg8eTNOmTbXkgyScmr9rG4rKn4iINAgz48ILL2TJkiW88MIL7N27l//8z/+kb9++vP766yqBItIomjZtypAhQ3TdnySc/Pz8Bv+aKn8iItKgkpKSGD16NMuXL+fpp59m69atfOc732HQoEHMmjVLJVBEGlwkEmHFihV8+eWXQUcROWa7du3itddea5SvrfInIiKNIiUlhTFjxvDRRx8xefJkiouLiUQiDBs2jMLCwqDjiUgC8TwPgFmzZgWcROTorFq1iv/5n/+hS5cupKenc9555wGVs2rMrMGmgIay/CUnJ5OZmRl0FBGRUEhNTWXs2LGsXLmSP/3pT3z00UcMGTKEkSNH8t577wUdT0QSQI8ePcjJyWHGjBmNcp2USEMrLS1l5syZjB8/ntNPP53OnTtz3XXXkZyczLhx4/jXv/4FUL3EksrfUSopKSE7O5ukpND96CIigWratCnXXnstn376Kffddx+LFi2if//+nH/++SxevLh6P71xE5EjlZSUxIgRI5g5c2ajXCcl0hCKi4t57LHHuPDCC8nOzsbzPCZNmkTHjh2ZOHEiK1eu5KOPPuKBBx6oPpvd0ELXgLTAu4hIsI477jhuuOEGVq1axYQJEygsLKRXr15cfPHFrFixQm/cROSoeJ7H+vXrg44hUq2srIzCwkJuueUWevbsSfv27Rk7diyLFi3isssu4+9//zubN29m+vTpXHfddXTp0mW/z8/NzW3wTCp/IiISiObNm3PbbbexevVq7rjjDqZPn0737t0BeP7559m1a1fACUUkXuTl5XHVVVdVP6+6TurWW28NMJWEQe3ZKhs3buSZZ55h9OjR5OTkMGTIEP7whz/QsmVLfv/737Ns2TI+++wzHn74Yc477zzS09Pr/bUbgsqfiIgEKjMzk6SkJL7++uvqO4FeeumlpKen07t3bwoLC3WHUBE5pLy8PJxzLF++HICTTjoJgAcffJDRo0fz6quvUlpaGmRESVD5+fnMnz+fvLw8+vXrR5s2bRgzZgxz5szhwgsv5C9/+QubNm1i9uzZ3HjjjXTr1g0zCyxv6Mrfpk2bVP5ERGJM1Ru3qpI3e/ZsrrjiCj7++GOGDBlCly5dyM/PZ/Xq1QEnFZFY1rVrVwBWr15NYWEhV1xxBbNmzeKCCy6gTZs2/PSnP2X27NmUl5cHnFTiWUlJCc899xyXX345AP369ePOO+8kOTmZ/Px8Fi5cyNq1a3nyySe56KKLyMjICDjx/wlV+SsrK2PLli0qfyIiMe6cc87hqaeeYsOGDTzzzDN06tSJ/Px8Tj75ZIYOHcqTTz7JV199FXRMEYlBubm5JCUl8R//8R9MmjSJtWvX8vrrr/Pd736XqVOnMnz4cNq3b8/48eOZP3++ZhbIYZWVlfH2229z++23c9ZZZ9GqVSsuu+wynn322ep9nHOMHDmSO+64g969e8fszSVjM1Uj2bx5M6A1/kREYlnNC9zT09O5/PLLmTlzJp999hl3330369ev56qrrqJNmzZcdtllzJgxQ0fxRaRa7eukUlNTOffcc3nmmWfYsGED06ZNo3///jz00EP069ePU089ldzcXD788MNgAktM+vzzz3nsscf4r//6L7Kzszn77LO55557aNKkCXfddRfvvfceZWVlQMMvx9CYLN6PdvTt29ctWLCgXvsuW7aMM844g2nTpnHxxRc3cjIROVpmttA51zfoHMfiSMYmOTLOOd59912mTJnC1KlT2bZtG23btuXyyy9nzJgxnH766UFHlASVCGMTaHyqsm3bNl5++WWef/55CgoKcM7Rq1cvfvjDHzJ69GjatWu33/55eXlx8eZejs7u3bspLCxk+vTpTJ8+nQ8++ACAdu3aMWrUKEaNGsWIESMOWCvczGLi7HF9x6dQnfkrKSkBdOZPRCSemRkDBgzg4YcfZt26dUybNo0zzzyT++67j29+85v079+fSZMmsWXLlv0+T2/aRKSmzMxMrrzySmbOnMmXX37Jgw8+SGpqKjfeeCMnnXQSQ4cOZfLkydUzx7QMTXyr/TvAOceHH37IxIkTOffcc2nZsiUjR45k0qRJtGvXjvvvv59ly5btdwawdvGDxlmOoTGF6szftGnT+P73v8+yZcvo1q1bIycTkaOVCEfXdWQ9+tavX8/zzz/P008/zdKlS0lLS+O73/0uY8aMYdSoUaSlpcXE0VmJX4kwNoHGp8P55JNPmDp1Ks8//zwffPABKSkpjBw5ktdee42FCxfSvXt30tLSgo4pR8jM2L59OwUFBdVn99asWQPAqaeeWn12b+jQoRx33HEBpz1y9R2fUqIRJlbozJ+ISOJq06YN48eP51e/+hWLFy9mypQpPP/88/z1r3+lVatWAOzdu1dv2kTkkLp06cLtt9/Obbfdxi9+8QsmT57Ma6+9BkCfPn0AOPHEEzn//PPp06cPffv2pVu3bhpbYtC+fftYsGABBQUFAGRlZVFWVsbxxx/PiBEjuPnmmxk5ciSdOnUKOGn0hK78mRktW7YMOoqIiDQSM6NXr1706tWL5s2bM2HCBDZu3AhAkyZNgMppOpoGKiKHYmY88sgjPPLIIzjnSEpK4sUXX2ThwoUsWLCAF198kcmTJwOQlpZGjx49qstgnz59VAgDUFFRwZIlSygoKKCgoIAZM2awb9++6terbtAybtw47rrrrqBiBip05a9FixakpITqxxYRCa277rqr+he8mdG2bVvWrVvHnj172LNnD02bNg04oYjEg6pFuS+55BIuueQSoPKasVWrVlWXwYULFzJ16tQDCmFVGTxYIdSNZI5e1XV7VWVvzpw51dd7n3baafzkJz9h+PDhDB06lFatWmnqPyEsf5ryKSISXsuXL+fXv/419957L6+88gpPPvkkAwcODDqWiMSB2jf2MDM6d+5M586dDyiEVWVw4cKFvPDCCzzyyCNAZSHs2bNndRns06cP+fn5Kn/15Jxj9erV1WVv9uzZrF+/HoAOHTpwwQUXMHz4cIYNG0bbtm0DThubVP5ERCQUcnNzycjI4PHHH+eSSy7hpz/9KYMHD2bcuHFMmDAhLi/wF5HoqU9Bq1kIv//97wOVheXTTz+tLoO1CyHATTfdxI033kh2dnZjxY8btc+Efvnll8yePbu68FXdpKVNmzYMHz68+uNw1+3F2105G0uo7vbZvXt3Tj31VF5++eVGTiUixyIR7qinu+nFvh07dnDTTTfx8MMP06VLF5544gmGDBkSdCyJYYkwNoHGp1iQm5vLnXfeecD2wYMH88orr5CVlRVAqthgZkybNq267H388ccAtGzZkmHDhlWXvdNOO616Oq5onb866cyfiIhUad68OZMmTaKgoICKigqGDh3Ktddey9dffx10NBFJcPn5+Tjnqq9BW7FiBaNHj2bu3Ll07NiR22677YC1ShNR1VnRP//5z/zyl7+kR48eQOW1lc899xynnHIK999/P++//z4lJSW89NJL/PKXv+T0009X8TtKoSl/FRUVbN68WeVPRET2M2zYMIqKirj++ut56KGHOOOMM5g1a1bQsSRKdK2VxIJvfvObvPDCCyxdupRvf/vb3HPPPXTs2JE77riDrVu3Bh2vwezevZvCwkJ+//vf873vfY82bdrQpUsXfvSjH/Hwww+zdOnS6n137NhB3759GT9+PGeeeSZJSaGpLY0qNH+LW7dupby8XOVPREQOkJ6ezoMPPkhhYSFpaWlEIhF+9rOfsX379qCjSSPLz88POoKEXM1r0bp168aLL75IUVERI0eOZMKEngHmLwAADttJREFUCXTs2JHc3Fy2bdsWYMqjU1xczLRp0xg3bhz9+/fnhBNOYMiQIdx0002sWLGCc889l8mTJ1NUVERZWVn1mdCqs6I6ONPwQnPDFy3wLiIihzN48GAWL15Mbm4u999/P6+//jqPPfYYo0aNCjpa6B3N7fArKirYtGkT69evr/7YsGHDfs+r9tNZBQlKXf+uu3fvzl/+8heKiorIz8/nzjvvZOLEifzqV79i3LhxZGRkRD9oDXX9f9y7dy+LFy9m3rx5zJ07l7lz51JcXAxAs2bN6NevHzfccAODBg1iwIABek8ekNDc8KWwsJAhQ4bwr3/9C8/zopBMRI5WItxUQTdUiH/vvvsuV155JStWrOCKK67gj3/8Iy1atAg6VmiZWfXZgK+++mq/Anewcrdx40bKy8sP+FopKSnViz3XlJube8iCmQhjE2h8ikeLFy8mPz+f//3f/yUzM5Px48dz/fXXc8IJJwSSx8zYsGFDddGbN28e8+fPZ8+ePQCcdNJJDBo0iEGDBjFw4EB69uxJampqvb621j08OvUdn3TmT0REpA79+/dn0aJF3Hnnndx777288cYbPPLII5x//vlBR4tpR/PGbd++fZSUlFBSUsLGjRur/6z5GKBjx45s2LCh+g1mTSkpKbRp04Y2bdrQtm1b+vTpQ5s2bWjdunX19qqP448/vvrzqkqlSCw788wz+dvf/sb7779PXl4ev/nNb3jggQf49a9/zXXXXUfz5s0b7Xs75yguLqaoqIglS5ZQVFQEQOvWrQFITU2lT58+/OIXv6gue8eyxp6KX+OK6pk/MxsFTASSgcedc7+r9XoT4BmgD7AZ+L5z7rNDfc36Hr2aPHkyP//5zykuLtaijyIxLhGOruvIemJZtGgRP/7xjykqKuKHP/whEydOjOv1uBrzyLqZsW/fPjZv3rxfeaur0FX9ebBrmQ5WzL71rW8xZsyY6jLXunVrWrRocVRTN4+k/CXC2AQanxLBwoULycvL4x//+ActW7bkhhtu4JprrjnmErhz506WL19OUVHRfmXvUNcb3nbbbUyYMOGYvq8cu5g782dmycBDgAcUA/PN7FXn3Ioau10FbHXOdTGz0cC9wPcb4vtXnfmL51/WIiISjN69ezN//nzuueceJkyYwMyZM3nooYe46KKLGqdIPfcc3HYbbs0arEMHuPtuuPTSBvvy+fn5h81cWlrKtm3bqj+2b9++3/O6PqruSpiWllZnoUpKSiI7O5tWrVqRk5NDr169yMnJqX7eqlWr/R5nZmZW3869sc7QaeFniUd9+vTh73//O/PnzycvL49bb72V+++/nxtvvJGrr76a41955ZBjiHOONWvWVJe7qo+VK1dW/z87/vjjOeOMMxg9ejQ9evSgR48edO/enYyMDJ0xj2NRO/NnZgOBPOfcSP/5LQDOuXtq7POGv888M0sB1gM57hAh63v0aty4cUycOFH/UEXiQCIcXdeR9cRVVFTElVdeycKFC7nooot46aWXGvZ3y3PPwdixsGvX/2077jh49NEDCmBFRQWlpaXs2bOHPXv2sHv37urHdT2v2nbttddy8803H7LI1TW1sqaUlBQyMzPJzMxk586drFu37oB9LrroIq655prqMteiRQuSk5OP6q8lFt5sJsLYBBqfEtG7775LXl4e06dPZ+zxx/P/SktJ3bev+vXyJk2Y/cMf8rdmzViyZAlLly7lq6++Air/b3Xu3Lm64PXs2ZMePXrQsWPHg55Nj4X/j7K/mDvzB7QFvqjxvBjof7B9nHNlZrYdyAI2Hes3rzrzJyIicix69OjBO++8w3333Vd99qxqymFycvIBf9a17VCvvfjuu7QpLd3/m+7axZdXXMHg227br8iV1t7vCPzud5VXXrRs2ZLOnTuTmZlJ+/btqwvd4T6aNWtW5yLLOkMnEn39+/fnn//8J/PmzePk4cP3K34AyaWldHnqKZ494QR69OjB5Zdfvt/ZvJrXwdaH/j/Gr2iWvwN/Q0Dt3w712QczGwuMhcq7CR1KXl7efmv4VP2iOtwdvUQkPMzsSeA8YKNzrnsdrxuV1yt/G9gFXOGcWxTdlBJLJkyYsN/vlqrrYfr06UOfPn0oLy+noqJivz/r2lbXn60OUuhOLCtj6NChNGvWjKZNm+73caTbWrVqRUVFRZ3lLVbpd7bI4Q0cOBB3kLP2Hagcqxri/73+P8avaJa/YqB9jeftgLUH2afYn/aZAWyp/YWcc48Cj0Ll1IVDfdOa12LoFLWIHMTTwP+j8oZTdTkXOMX/6A88zIEzFyREGvV3S8eOsGbNAZuTOnRgypQpDfZtGqv46YyASLCsQ4c6xxDr0AHi6ICPNI5ormg6HzjFzDqZWRowGni11j6vAmP8xxcBBYe63k9EpCE4596ijgNNNVwAPOMqvQNkmtmJ0UknoXP33ZXX+NV03HGV2xtIYxY0nREQCVgUxhCJX1Erf865MuAa4A3gA2Cac265md1pZlWLJj0BZJnZJ8B44OaGzKCjkSJylOq6ZvmANWPMbKyZLTCzBbrOODwa/HfLpZdW3tylQ4fK6x46dKjzZi/HQgVNJIFFYQyR+BXVdf4ag+5YJZJ4grijnpl1BP5xkGv+XgPucc697T+fBfy3c27hwb6exiaRxKO7fYpIrKrv+BTNaZ8iIvGqPtcsi4iIiMQ0lT8RkcN7FfiRVRoAbHfOHbiomYiIiEgMi+bdPkVEYpKZvQCcA2SbWTGQC6QCOOceAV6ncpmHT6hc6uHHwSQVEREROXoqfyISes65HxzmdQdcHaU4IiIiIo1C0z5FRERERERCQOVPREREREQkBFT+REREREREQkDlT0RERCSKzKypmb1nZkvMbLmZ5fvbO5nZu2a20sxeNLM0f3sT//kn/usda3ytW/ztH5nZyGB+IhGJFyp/IiIiItFVCgx3zvUEzgRG+cvI3As84Jw7BdgKXOXvfxWw1TnXBXjA3w8z6wqMBroBo4BJZpYc1Z9EROKKyp+IiIhIFLlKX/tPU/0PBwwHXvK3TwG+5z++wH+O//oIMzN/+1TnXKlzbjWVy9H0i8KPICJxSuVPREREJMrMLNnMFgMbgRnAp8A251yZv0sx0NZ/3Bb4AsB/fTuQVXN7HZ9T+/uNNbMFZragpKSkoX8cEYkTKn8iIiIiUeacK3fOnQm0o/Js3Tfr2s3/0w7y2sG21/X9HnXO9XXO9c3JyTmayCKSAFT+RERERALinNsGzAEGAJlmluK/1A5Y6z8uBtoD+K9nAFtqbq/jc0REDqDyJyIiIhJFZpZjZpn+42ZABPgAmA1c5O82BnjFf/yq/xz/9QLnnPO3j/bvBtoJOAV4Lzo/hYjEo5TD7yIiIiIiDehEYIp/Z84kYJpz7h9mtgKYamYTgPeBJ/z9nwD+bGafUHnGbzSAc265mU0DVgBlwNXOufIo/ywiEkes8sBR/DKzEmBNPXfPBjY1YpzGoMyNL97yQuJn7uCci+uLUjQ2xSRljo54yxyqsQk0PsWgeMsLyhwtDT4+xX35OxJmtsA51zfoHEdCmRtfvOUFZU408fh3o8zRocyNL97yRls8/v3EW+Z4ywvKHC2NkVnX/ImIiIiIiISAyp+IiIiIiEgIhK38PRp0gKOgzI0v3vKCMieaePy7UeboUObGF295oy0e/37iLXO85QVljpYGzxyqa/5ERERERETCKmxn/kREREREREIpNOXPzEaZ2Udm9omZ3Rx0nkMxs/ZmNtvMPjCz5WZ2fdCZ6svMks3sfTP7R9BZ6sPMMs3sJTP70P/7Hhh0psMxs1/5/y6WmdkLZtY06Ey1mdmTZrbRzJbV2NbSzGaY2Ur/zxZBZowV8TQ2QfyOTxqbGp/GpsSisSl6ND41Lo1N+wtF+fMXUX0IOBfoCvzAzLoGm+qQyoBfO+e+CQwAro7xvDVdD3wQdIgjMBGY7pw7HehJjGc3s7bAdUBf51x3IBl/sd8Y8zQwqta2m4FZzrlTgFn+81CLw7EJ4nd80tjUiDQ2JRaNTVGn8amRaGw6UCjKH9AP+MQ5t8o5txeYClwQcKaDcs6tc84t8h/voPI/VdtgUx2embUDvgM8HnSW+jCzE4AhwBMAzrm9zrltwaaqlxSgmZmlAMcBawPOcwDn3FvAllqbLwCm+I+nAN+LaqjYFFdjE8Tn+KSxKWo0NiUOjU1RovEpKjQ21RCW8tcW+KLG82LiYEAAMLOOQC/g3WCT1MuDwH8DFUEHqaeTgRLgKX+6xeNmlh50qENxzn0J/AH4HFgHbHfO/SvYVPXW2jm3Dip/SQOtAs4TC+J2bIK4Gp80NjUyjU0JR2NT9Gh8akQamw4UlvJndWyL+ducmtnxwF+Bcc65r4LOcyhmdh6w0Tm3MOgsRyAF6A087JzrBewkxqf7+PO9LwA6Ad8A0s3ssmBTyTGIy7EJ4md80tgUHRqbEo7GpijQ+NT4NDYdKCzlrxhoX+N5O2LwlG9NZpZK5eD1nHPu5aDz1MNg4Hwz+4zK6SHDzezZYCMdVjFQ7JyrOjL4EpUDWiyLAKudcyXOuX3Ay8CggDPV1wYzOxHA/3NjwHliQdyNTRB345PGpujQ2JRYNDZFh8anxqexqZawlL/5wClm1snM0qi80PPVgDMdlJkZlXOpP3DO/THoPPXhnLvFOdfOOdeRyr/fAudcTB9Zcc6tB74ws9P8TSOAFQFGqo/PgQFmdpz/72QEMXyhdS2vAmP8x2OAVwLMEiviamyC+BufNDZFjcamxKKxKQo0PkWFxqZaUhrii8Q651yZmV0DvEHlXX6edM4tDzjWoQwGLgeWmtlif9utzrnXA8yUqK4FnvN/ua0CfhxwnkNyzr1rZi8Bi6i8s9n7wKPBpjqQmb0AnANkm1kxkAv8DphmZldRORhfHFzC2BCHYxNofIoWjU2NQGNT/WhsksOIm/FJY1Md38u5uJjCLSIiIiIiIscgLNM+RUREREREQk3lT0REREREJARU/kREREREREJA5U9ERERERCQEVP5ERERERERCQOVPREREREQkBFT+REREREREQkDlT0REREREJAT+P+Cws1/KQD6LAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "n,p = X_train.shape\n", "\n", "R2score = []\n", "AICscore = []\n", "BICscore = []\n", "\n", "lr = LinearRegression()\n", "lrfit = lr.fit(X_train,y_train)\n", "yhat = lrfit.predict(X_train)\n", "\n", "shat = sigmaHat(X_train, y_train)\n", "\n", "for d in range(p+1):\n", " if d == 0:\n", " yhat = y_train.mean() * np.ones_like(y_train)\n", " else:\n", " l = best_model[d]\n", " lr = LinearRegression()\n", " lrfit = lr.fit(X_train[:,l],y_train)\n", " yhat = lrfit.predict(X_train[:,l])\n", " R2score.append( adjustedRSquare(y_train, yhat, d) )\n", " AICscore.append( AIC(y_train, yhat, d, shat) )\n", " BICscore.append( BIC(y_train, yhat, d, shat) )\n", "\n", "R2idx = np.argmax(R2score)\n", "AICidx = np.argmin(AICscore)\n", "BICidx = np.argmin(BICscore)\n", "\n", "xr = range(p+1)\n", "\n", "plt.rcParams['figure.figsize'] = (15,8)\n", "fig, ax = plt.subplots(1,3)\n", "\n", "ax[0].plot(xr,R2score,'k+-')\n", "ax[0].plot(xr[R2idx], R2score[R2idx], 'ro')\n", "ax[0].set_title('Number of pred.: %i (%6.4f)' % (R2idx, R2score[R2idx]))\n", "ax[0].set_ylabel('$R^2$')\n", "\n", "ax[1].plot(xr,AICscore,'k+-')\n", "ax[1].plot(xr[AICidx], AICscore[AICidx], 'ro')\n", "ax[1].set_title('Number of pred.: %i (%6.4f)' % (AICidx, AICscore[AICidx]))\n", "ax[1].set_ylabel('AIC')\n", "\n", "ax[2].plot(xr,BICscore,'k+-')\n", "ax[2].plot(xr[BICidx], BICscore[BICidx], 'ro')\n", "ax[2].set_title('Number of pred.: %i (%6.4f)' % (BICidx, BICscore[BICidx]))\n", "ax[2].set_ylabel('BIC')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Extra Task**: Implement the best subset selection algorithm as one function `bestSubsetSelection`.\n", "You might use the function `bestSubsetComputation` as well as your scoring function(s)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Solution**:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.5225204821123409\n", "(1, 2, 3, 4, 5, 7, 8)\n" ] } ], "source": [ "def bestSubsetSelection(X, y, scoring = 'aic'):\n", " \"\"\" Input: X - predictor array of size (n,p)\n", " y - array of size (n,)\n", " scoring - string variable, one of 'aic', 'bic', 'r2'\n", " Output: bps - best parameter selection\n", " \"\"\"\n", " score = []\n", " def scoring_func(y, yhat, d):\n", " if scoring == 'r2':\n", " return adjustedRSquare(y, yhat, d)\n", " else:\n", " shat = sigmaHat(X,y) \n", " if scoring == 'aic':\n", " return AIC(y, yhat, d, shat)\n", " elif scoring == 'bic':\n", " return BIC(y, yhat, d, shat)\n", " else:\n", " raise NameError('scoring not known')\n", "\n", " best_model, best_score = bestSubsetComputation(X, y)\n", " \n", " for d, l in enumerate(best_model):\n", " if d == 0:\n", " yhat = y.mean() * np.ones_like(y)\n", " else:\n", " lr = LinearRegression()\n", " lrfit = lr.fit(X[:,l],y)\n", " yhat = lrfit.predict(X[:,l])\n", "\n", " score.append( scoring_func(y, yhat, d) )\n", " \n", " if scoring == 'r2':\n", " best_idx = np.argmax(score)\n", " else:\n", " best_idx = np.argmin(score)\n", " \n", " return score[best_idx], best_model[best_idx]\n", "\n", "# Test the implementation\n", "optimal_score, optimal_model = bestSubsetSelection(X_train,y_train,'r2')\n", "print(optimal_score)\n", "print(optimal_model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Homework 7" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Implement either the **forward stepwise selection algorithm** or the **backward stepwise selection algorithm** as a function `forwardStepwiseSelection` or `backwardStepwiseSelection`, resp.\n", "\n", "Use 10-fold cross-validation as a measure for model selection (step 3 of the algorithms).\n", " \n", "Test your algorithm on the diabetes data set from above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }