{ "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": [ "**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": [ "**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": [ "**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": [ "### 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": [ "**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": [ "**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": [ "# 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 }