07_cross_validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cross-validation for parameter tuning, model selection, and feature selection\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Agenda\n",
    "\n",
    "- What is the drawback of using the **train/test split** procedure for model evaluation?\n",
    "- How does **K-fold cross-validation** overcome this limitation?\n",
    "- How can cross-validation be used for selecting **tuning parameters**, choosing between **models**, and selecting **features**?\n",
    "- What are some possible **improvements** to cross-validation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review of model evaluation procedures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Motivation:** Need a way to choose between machine learning models\n",
    "\n",
    "- Goal is to estimate likely performance of a model on **out-of-sample data**\n",
    "\n",
    "**Initial idea:** Train and test on the same data\n",
    "\n",
    "- But, maximizing **training accuracy** rewards overly complex models which **overfit** the training data\n",
    "\n",
    "**Alternative idea:** Train/test split\n",
    "\n",
    "- Split the dataset into two pieces, so that the model can be trained and tested on **different data**\n",
    "- **Testing accuracy** is a better estimate than training accuracy of out-of-sample performance\n",
    "- But, it provides a **high variance** estimate since changing which observations happen to be in the testing set can significantly change testing accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn import metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the iris data\n",
    "iris = load_iris()\n",
    "\n",
    "# create X (features) and y (response)\n",
    "X = iris.data\n",
    "y = iris.target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.973684210526\n"
     ]
    }
   ],
   "source": [
    "# use train/test split with different random_state values\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=4)\n",
    "\n",
    "# check classification accuracy of KNN with K=5\n",
    "knn = KNeighborsClassifier(n_neighbors=5)\n",
    "knn.fit(X_train, y_train)\n",
    "y_pred = knn.predict(X_test)\n",
    "print(metrics.accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** What if we created a bunch of train/test splits, calculated the testing accuracy for each, and averaged the results together?\n",
    "\n",
    "**Answer:** That's the essense of cross-validation!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Steps for K-fold cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Split the dataset into K **equal** partitions (or \"folds\").\n",
    "2. Use fold 1 as the **testing set** and the union of the other folds as the **training set**.\n",
    "3. Calculate **testing accuracy**.\n",
    "4. Repeat steps 2 and 3 K times, using a **different fold** as the testing set each time.\n",
    "5. Use the **average testing accuracy** as the estimate of out-of-sample accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diagram of **5-fold cross-validation:**\n",
    "\n",
    "![5-fold cross-validation](images/07_cross_validation_diagram.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration                   Training set observations                   Testing set observations\n",
      "    1     [ 5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]        [0 1 2 3 4]       \n",
      "    2     [ 0  1  2  3  4 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]        [5 6 7 8 9]       \n",
      "    3     [ 0  1  2  3  4  5  6  7  8  9 15 16 17 18 19 20 21 22 23 24]     [10 11 12 13 14]     \n",
      "    4     [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 20 21 22 23 24]     [15 16 17 18 19]     \n",
      "    5     [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19]     [20 21 22 23 24]     \n"
     ]
    }
   ],
   "source": [
    "# simulate splitting a dataset of 25 observations into 5 folds\n",
    "from sklearn.cross_validation import KFold\n",
    "kf = KFold(25, n_folds=5, shuffle=False)\n",
    "\n",
    "# print the contents of each training and testing set\n",
    "print('{} {:^61} {}'.format('Iteration', 'Training set observations', 'Testing set observations'))\n",
    "for iteration, data in enumerate(kf, start=1):\n",
    "    print('{:^9} {} {:^25}'.format(iteration, data[0], data[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Dataset contains **25 observations** (numbered 0 through 24)\n",
    "- 5-fold cross-validation, thus it runs for **5 iterations**\n",
    "- For each iteration, every observation is either in the training set or the testing set, **but not both**\n",
    "- Every observation is in the testing set **exactly once**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing cross-validation to train/test split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Advantages of **cross-validation:**\n",
    "\n",
    "- More accurate estimate of out-of-sample accuracy\n",
    "- More \"efficient\" use of data (every observation is used for both training and testing)\n",
    "\n",
    "Advantages of **train/test split:**\n",
    "\n",
    "- Runs K times faster than K-fold cross-validation\n",
    "- Simpler to examine the detailed results of the testing process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation recommendations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. K can be any number, but **K=10** is generally recommended\n",
    "2. For classification problems, **stratified sampling** is recommended for creating the folds\n",
    "    - Each response class should be represented with equal proportions in each of the K folds\n",
    "    - scikit-learn's `cross_val_score` function does this by default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation example: parameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Goal:** Select the best tuning parameters (aka \"hyperparameters\") for KNN on the iris dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import cross_val_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.          0.93333333  1.          1.          0.86666667  0.93333333\n",
      "  0.93333333  1.          1.          1.        ]\n"
     ]
    }
   ],
   "source": [
    "# 10-fold cross-validation with K=5 for KNN (the n_neighbors parameter)\n",
    "knn = KNeighborsClassifier(n_neighbors=5)\n",
    "scores = cross_val_score(knn, X, y, cv=10, scoring='accuracy')\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.966666666667\n"
     ]
    }
   ],
   "source": [
    "# use average accuracy as an estimate of out-of-sample accuracy\n",
    "print(scores.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.95999999999999996, 0.95333333333333337, 0.96666666666666656, 0.96666666666666656, 0.96666666666666679, 0.96666666666666679, 0.96666666666666679, 0.96666666666666679, 0.97333333333333338, 0.96666666666666679, 0.96666666666666679, 0.97333333333333338, 0.98000000000000009, 0.97333333333333338, 0.97333333333333338, 0.97333333333333338, 0.97333333333333338, 0.98000000000000009, 0.97333333333333338, 0.98000000000000009, 0.96666666666666656, 0.96666666666666656, 0.97333333333333338, 0.95999999999999996, 0.96666666666666656, 0.95999999999999996, 0.96666666666666656, 0.95333333333333337, 0.95333333333333337, 0.95333333333333337]\n"
     ]
    }
   ],
   "source": [
    "# search for an optimal value of K for KNN\n",
    "k_range = list(range(1, 31))\n",
    "k_scores = []\n",
    "for k in k_range:\n",
    "    knn = KNeighborsClassifier(n_neighbors=k)\n",
    "    scores = cross_val_score(knn, X, y, cv=10, scoring='accuracy')\n",
    "    k_scores.append(scores.mean())\n",
    "print(k_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x8b69ef0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEPCAYAAACHuClZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8VOWd5/HP9yIggiyKiLKJMaKAsukVSKZzE2O7jRpN\npxM7GbfE2N0aEzvTo3E6I2bSHU26HbNMOtoxGcdOYmJiopmk3b3ZuBdBLqhsalAEBFwIgqDI8ps/\nnnOgKKruPbWcqjp1f+/X676oe5annkPB+dXz/M7zPDIznHPOuUq01LsCzjnnss+DiXPOuYp5MHHO\nOVcxDybOOecq5sHEOedcxTyYOOecq1jqwUTSGZKWS3pW0rUF9g+VdK+kxZI6JU3M2XeNpGckPSXp\nB5L6RduHSXpI0gpJD0oakvZ1OOecKy7VYCKpBfgWcDowCbhQ0nF5h10PdJnZFOBi4BvRuUcCnwGm\nm9mJwAHAx6JzrgMeMbMJwGPAF9K8Duecc91Lu2XSCjxnZqvMbAdwN3Be3jETCQEBM1sBHCXpsGhf\nH2CgpAOAg4C10fbzgDuj13cCH0rvEpxzzvUk7WAyClid8/uaaFuuxcAFAJJagbHAaDN7GfgX4CVC\nENlkZo9G54wwsw0AZrYeGJHaFTjnnOtRIyTgbwKGSVoIXAl0AbskDSW0QMYBRwKDJP1VkTJ8Thjn\nnKujA1Iufy2hpREbzd6uKgDMbAtwWfy7pJXASuAMYKWZbYy23wvMBn4IbJB0uJltkDQSeKXQm0vy\nIOOcc2UwM5VyfNotk/nAMZLGRU9ifQy4P/cASUMk9Y1eXw781szeJHRvzZR0oCQBpwLLotPuBy6J\nXl8M3FesAmbWtD833HBD3evQrNe3e7dx2GHGwIHGjh3NdW1mxksvGWDMnl39sk86yRg6tL7Xl/ZP\nvT+/tH/KkWowMbNdwFXAQ8AS4G4zWybpCkmfjg47HnhG0jLCU1+fjc59AvgpodtrMSDg9uicm4HT\nJK0gBJmb0rwO1/v88Y/Qrx+MHQtPPVXv2lTf3LnwgQ/AokWwfXv1yt22DZYsgc2bYffu6pXrGl/a\n3VyY2QPAhLxtt+W87szfn7PvRuDGAts3Ah+sbk2d26ujA2bPhiFDwuvp0+tdo+rq6IDTT4eNG6Gr\nC2bOrE65CxbAiSeGAPzKKzByZHXKdY2vERLwrkxtbW31rkKq6nl9c+fCrFnhZ+7c6pdf788ureuL\nyx01qo01a6pXbqOp9+fXiFRu/1gWSLJmvj6XnqlT4bbbQsvkrLNg5cp616h63noLhg+H116Dn/0M\n7rsP7rmnOmWfey78l/8Cd90Fl14K559fnXJdbUnCGiwB71zmbN4Mzz0H06bBscfCpk2wbl29a1U9\nCxbApEkwYEDoyps7F6rxncssdJ/NmgVjxsDq1T2f45qHBxPn8jzxRAgk/fpBS0u4OXZ01LtW1TN3\nbggiAOPHw86d8NJLlZf7/PMhQI0e7cGkN/Jg4lyeOPkemz27uYJJ7vVJ1bu+3HI9mPQ+HkycyxMn\nkWNpJeHrwSy968std8wYmjoB7/bnwcS5HLt3Q2fnvjfb1lZYvLi64zHqZeVK6N8/3Oxj1WqZ5Haf\njR7tLZPexoOJczmWL4dhw/YdHzFoELz73WE8Rtblt0oAZsyApUth69byy928OQSqKVPC76NGhYcW\ndu0qv0yXLR5MnMuR++06V/zUU9YVur4BA2Dy5PCUV7nmzQsDO/v1C7/37w+HHALr15dfpssWDybO\n5chPvseaJQmf1vUVKteT8L2LBxPnchTqBoK9Seosj4HdsiU8vjt16v77Km15Ffp78yR87+LBxLnI\nxo3h5nfCCfvvGz8+9P9n+Zt27viZfPFYmnKCZaGHFsCT8L2NBxPnIp2dcPLJcECB6U+l7D8iXKzV\nBeHGP2BAaLmUatmyMD3LiLz1Tr2bq3fxYOJcpFjyPZb1JHxP11fuSP9i5Xow6V08mDgXKZacjmU5\nCV+sKypXucGy2N+bB5PexYOJc4T5qZ54ovt1PeLxGNu21a5e1bJ8eXhU9/DDix9TbjAp1n3mCfje\nxYOJc8Azz4S8wSGHFD/mwANDcr6S8Rj10lOrC8KAw5UrwwDEpF5/HV5+OYxTyXfEEWGBrJ07S6ur\nyyYPJs7RfXI6V1aT8Emur1+/MPBw3rzk5XZ2hulm+vTZf1/fvnDYYSHYuObnwcQ5kn1zh+wm4XtK\nvsdKzQv1VK7nTXoPDybOkfxmW8l4jHrZuBHWri3cFZWv1JZXT0HYg0nv4cHE9XobNoQb7nHH9Xxs\nJeMx6qW78TP5Zs0Kx+/e3fOxO3fC/PlwyinFj/EkfO/hwcT1eh0d4SmuloT/G7L2iHDSLjwIAw+H\nDw8DEXvy1FMwdmyYZbkYHwXfe3gwcb1e0i6uWNbyJmldX7zee3e8m6v38GDier0kN8VcWXqiK+6K\n6m78TL6kLa8kQcqDSe/hwcT1au+8AwsXhsdbk5o6tfTxGPUSj5/prisqX9JgmaT7zINJ7+HBxPVq\nXV1wzDEweHDyc/r2LX08Rr0kHT+Ta/LksEri668XP2bdOti0CY49tvuyRo4MDze8805pdXDZ48HE\n9WqlJKdzZSUJX8719ekTWmqdnd2XO2tWzw8t9OkTAsrataXVwWWPBxPXq5WanI5lJQlf7vX11NVV\nSp7Ju7p6h9SDiaQzJC2X9KykawvsHyrpXkmLJXVKmhhtP1ZSl6SF0Z9vSLo62neDpDXRvoWSzkj7\nOlxzKjX5Hps5M/l4jHrZsAH+9CeYMKH0c3tqeZUSpDyY9A6pBhNJLcC3gNOBScCFkvKHhl0PdJnZ\nFOBi4BsAZvasmU0zs+nADGArcG/OebeY2fTo54E0r8M1p9WrYft2eNe7Sj93xIgw71SS8Rj1Uur4\nmVynnBKeAis0SeP27bBoUfKHFjyY9A5pt0xagefMbJWZ7QDuBs7LO2Yi8BiAma0AjpJ0WN4xHwT+\naGa5Y2mVUp1dLxEnp1Xmv6RGf0S4nOR7bNiwMCDx6af339fVFRLvgwYlK8tHwfcOaQeTUUDud5I1\n0bZci4ELACS1AmOB0XnHfBT4Ud62qyQtkvRdSUOqV2XXW5SbfI81ehK+GtdXKFiWmofxUfC9QyMk\n4G8ChklaCFwJdAG74p2S+gLnAvfknPNt4GgzmwqsB26pXXVdsyg3OR1r5CT8O++EFkQp42fyFWt5\nlZpn8m6u3iHB1G8VWUtoacRGR9v2MLMtwGXx75JeAFbmHHIm8KSZvZpzzqs5+/8N+GWxCsyZM2fP\n67a2Ntra2kqpv2tSb70FS5bASSeVX8akSXvHYxx6aPXqVg3x+JmDDy6/jNmz4ctf3nebWQgwN9+c\nvBwPJo2vvb2d9vb2isqQpTiXtqQ+wArgVGAd8ARwoZktyzlmCLDNzHZIuhx4j5ldkrP/R8ADZnZn\nzraRZrY+en0NcLKZ/VWB97c0r89l1+9+B5//fFiqtxKnnQaf+xycfXZ16lUtt94Kzz4L3/52+WXs\n3h0mfVy6NIwVAVi1KiTn161LnmvavTvMtLxpU/jTNT5JmFlJ2cRUu7nMbBdwFfAQsAS428yWSbpC\n0qejw44HnpG0jPDU12fj8yUdREi+37tvyXxV0lOSFgHvA65J8zpc8yn3keB8jZqEryT5Hmtp2bt+\nSyz+eyvloYWWFhg1ygcuNru0u7mIHtudkLfttpzXnfn7c/ZtA/Kf7MLMLqpyNV0vM3cuXHhh5eXM\nng1f/Wrl5VRbRwd85SuVlxPnhc4/P/xebp4pTsIfc0zldXKNqRES8M7VVNzvX0nyPTZzZvHxGPWy\nenVIwB99dOVlFWqZlPP35nmT5ufBxPU6K1dCv37hBlepoUPDeIynnqq8rGqpdPxMrtbWMEBx+3bY\nujXkT2bMKL0cDybNz4OJ63Wq1SqJNdp4k0rHl+QaNCgMUOzqggUL4IQT4MADSy/Hg0nz82Diep1q\nJd9jjZaEr0byPVfc1VXJ35uPgm9+Hkxcr9PMLZNqjJ/JFyfhK/l781Hwzc+DietVtmyB55+HadOq\nV+axx8Ibb4SxF/W2YEFY3Kqa4znillelLRMPJs3Ng4nrVZ54Iiy7269f9cpsaQlPdTVC66TaXVwA\n48fDrl0hQI3OnzUvoeHDQ6tp69bq1s01jtTHmTgH4XHc5cvh+OPrW49qd3HFZs+GH/84LOlbT7/6\nFVx1VXXLlML1VRKApb1dXcflL0KRkuXLw1ou1XiqzfWsx2Ai6RzgV2bWwMsAuUa3aBG8972hO+iA\nOn6F6eiAyy+vfrnnnx/WhL/99uqXXYpDD4UPfKD65X7qU2EJ3krESfhaBJNdu8K0L52d9f8C01sk\n+W/9UeBWST8Dvmdmy1Ouk2tCc+fCtm1hPMb06fWpw+7dIZh873vVL3viRLj//uqX2yjOOqvyMmqZ\nhF+yBDZvhpde8mBSKz3mTMzsE8A04I/A/5HUIenTkiqYj9T1NnPnwuDB9X2EdvlyOOSQvZMWutqq\nZRI+/nfmSf/aSZSAN7PNwE8JKyUeAZwPLJT0mRTr5ppIRwdccUV9k9TVHl/iSlPLYNLRAUcc4cGk\nlnoMJpLOlfRzoB3oC7Sa2ZnAFODz6VbPNYP168P045dcUt+WSVrJd5dMrVsmH/mID5SspSQtkw8D\n/8vMTjCzr5nZK7BnRt9Pplo71xQ6OsKjs8cdV9/xGN4yqa9ajYJ/9dXwc/rp3jKppSTBZA5hUSsA\nJA2QdBSAmT2aSq1cU4lbBIXWx6iVjRvDjeWEE2r/3i6oVQK+oyM8yTVunAeTWkoSTO4Bch8L3sW+\n67E7163c7qV6TT3S2RlmwK3nY8m93bBhYar+zZvTfZ+4BRp3q/liq7WRJJgcYGbvxL9Er6s4ftg1\ns+3bwxiT1tbwe70mRfQurvqTapM3ib+8DB4cWsObNqX7fi5IEkxelXRu/Iuk84DX0quSayZdXWHu\nqkGDwu+562PUkiffG0PawWTHDnjyydDNFb+fJ+FrI0kw+WvgekkvSVoNXAtckW61XLPIbxHkro9R\nKzt3htUQZ86s3Xu6wtK+uS9eHOYSGzJk7/t53qQ2euxBNrM/AjMlDYp+fzP1WrmmMXcunHfevtvi\nKc1rdXN/5hkYNSoMWHT1lXYSPr8F6sGkdhKlIyWdDUwCDlQ0a5qZfSnFerkmEK+1fvPN+26fPRt+\n8Yva1cO7uBrHmDHpPoDR0REeCc59Pw8mtZFk0OJ3CPNzfQYQ8BFgXMr1ck1g9erQxTR+/L7b4yR8\nrZ6y8eR740j75p7/xcEX5aqdJDmT2WZ2EfAnM7sRmAUcm261XDOI/2PnTwEer4/x0ku1rYervzSD\nydq1Yb2Ud7973/fzBHxtJAkmb0d/bpN0JLCDMD+Xc90q1iKQajd4ccOGMGCxVmtouO7FN/c0WqXx\nv7fcLy/ezVU7SYLJLyUNBb4GLAReBH6YZqVcc+iuRRAn4dMWT+XS4muKNoQ0x34U+veWZvBy++r2\nv5ikFuBRM9tkZj8j5EqOM7P/UZPauczatg2WLoUZMwrvr9VIeO/iajxptRYKtYQHDoQDD4TXX6/+\n+7l9dRtMotUV/3fO79vN7I3Ua+Uyb8ECmDw5rBteyIwZIdikvSa4J98bTxrB5O23w8JrJ5+8/z5P\nwtdGksb/o5I+LPlKyi65nloEBx4YJl1csCC9OrzzDixcuHcqF9cY0ggmTz4ZVlQcOLDw+3kSPn1J\ngskVhIkdt0vaLGmLpJSnanNZN3duzy2CtJPwXV1wzDGhn941jjRu7t21QD0JXxtJlu092MxazKyf\nmQ2Ofk/831PSGZKWS3pW0rUF9g+VdK+kxZI6JU2Mth8rqUvSwujPNyRdHe0bJukhSSskPShpSCkX\n7dJlFv5z95SrSDsJn6QOrvbS6HbqriXswaQ2kgxa/LNCP0kKjxL43wJOJ4ygv1BS/kOa1wNdZjYF\nuBj4BoCZPWtm08xsOjAD2ArcG51zHfCImU0AHgO+kKQ+rjaefz7kSkaP7v64uGWS1pM2nnxvTNW+\nucdfXrxlUl9Jurn+Pufni8AvCQtmJdEKPGdmq8xsB2EN+byZmphICAiY2QrgKEmH5R3zQeCPZhY3\njs8D7oxe3wl8KGF9XA0kTXqPHh2CzvPP17cerraqfXN/8cUwtmRckXk5PAFfG0m6uc7J+TkNmAz8\nKWH5o4Dcj3FNtC3XYuACAEmtwFgg/zvtR4Ef5fw+wsw2RPVbD4xIWB9XA6W0CNJ6RHj16jDN/bve\nVf2yXWWqPfYj7s4s9oiQJ+Bro5x159YAx1exDjcBX5e0EHga6CKs5giApL7AuYSurWKK/rOcM2fO\nntdtbW20tbVVVlvXo7lz4bLLkh0bz9N10UXVr0P+aGjXGAYODC3S11+H4cMrL6+nhz1Gjw5Treze\n7YNXi2lvb6e9vb2iMnoMJpK+yd6bdQswlTASPom1hJZGbHS0bQ8z2wLsufVIegFYmXPImcCTZvZq\nzrYNkg43sw2SRgKvFKtAbjBx6du8GVauhKlTkx0/ezbccUf16+HJ98YWdz1VK5h8/OPF9w8YAAcf\nDK++CocfXvn7NaP8L9o33nhjyWUkidMLgCejnw7gWjP7RMLy5wPHSBonqR/wMeD+3AMkDYlaH0i6\nHPhN3popF7JvFxdRGZdEry8G7ktYH5eyefNg+nTol3Bh56lTQ/Cp9rrgSR5NdvVTrbzJm2/CihXh\n31wt3s8Vl6Sb66fA22a2C0BSH0kHmdm2nk40s12SrgIeIgSuO8xsmaQrwm67ndBldqek3cAS4JPx\n+ZIOIiTfP51X9M3ATyRdBqwC/jLBdbgaKDXp3bdvuBHMmwennVadOrz1FixZAiedVJ3yXPVV6+Y+\nfz5MmQL9+3d/XNwS8n8T6UkSTB4l3NDj1sIAQnBI1IlgZg8AE/K23ZbzujN/f86+bUD+k12Y2cao\nTq7BzJ0Lf/M3pZ0TJ+GrFUwWLIBJk+Cgg6pTnqu+agWTpN2ZnoRPX5JurgNzu52i1/7f1O1n927o\n7Cy9eylOwleLd3E1vmrd3JN+1t7Nlb4kwWSrpD09kpJmAG+lVyWXVcuWhYTqiBIf1J41KwSh3bur\nUw9Pvje+aoz96GmwYi4PJulLEkw+B9wj6XeSfg/8GLgq3Wq5LCp3xPmIEXDYYSEYVSped95bJo2t\nGjf3Z58NT2kdeWRt3s91r8eciZnNj6ZAifMaK6LR7M7to5IR53FX16RJldVh5crwJNmYMZWV49JV\njbEfpXx58VHw6UsyN9eVwEAze8bMngEGSfrb9KvmsqaSubCqNRK+2LrzrrEMGBBmc36l6AixnpXS\nnTlqFKxbB7t29XysK0+S7wSXm9meRTbN7E/A5elVyWXR66/Dyy+HBbHKUa0kvHdxZUelSfhSPuv+\n/eGQQ2DDhvLfz3UvSTDpk7swlqQ+QMIhaa636OwMi1D16VPe+ZMnh2BU6fKqnnzPjkq6njZtglWr\n4MQTk5/jeZN0JQkmDwA/lnSqpFMJo9EfSLdaLmsqne69Tx845ZQQlMq1ZQs89xxMm1Z+Ga52Krm5\nz5sXln7u27c27+d6liSYXEuYIv5vop9HCdPRO7dHNaZ7r7Sr64knQiBJOpWLq69Kbu7lfHnxJHy6\nkkxBv9vMvmNmf2FmfwH8Gvh8+lVzWbFzZ5jWYubMysqpNAnvi2FlSyXBpJzuTB8Fn65ED+VJOkzS\n30r6HdAO+Nybbo+nngr/UYcNq6ycU04JQWnnzvLO9+R7tpR7c9+1K3Rzlfrlxbu50lU0mEg6WNLF\nkh4EngDeBYw3s3eZ2X+tWQ1dw6tW0nvYMBg7Fp5+uvRzy53KxdVPud1OS5fCyJGlT1/vwSRd3bVM\nXiGsM/Jl4Ggz+zzwTk1q5TKlmt1Ls2eXlzdZvjwEo5Ejq1MPl75yx36U2wL1YJKu7oLJF4D+wLeB\nL0jyBVBdQdXsXio3Ce/5kuyJx36sX1/aeeV+1kccEQZJltuN6rpXNJiY2a1mNhM4L9r0C+BISddK\nOrYmtXMNb906eOMNmFBwEYHSlZuE9/El2VROa6Hcz7pv3zAH3Lp1pZ/repbkaa6VZvZPZnYCcBIw\nmPBEl3N0dIREaLXW1j722DAgrZxvq54vyZ5Sk/CvvhpaFxMnlv9+3tWVjpJuAdH8XP/dzI5Jq0Iu\nW6rdImhpCUGhlNbJxo3hhnTCCdWrh6uNUpPwnZ3hqb9yv7x4MElPlb5Put4qjRZBqXmTzk44+WQ4\nIMm6oa6hlHpzr/TfmweT9HgwcWXbvh0WLQpzclVTqU90efI9u8oJJpV81j4KPj3+XS5F27fDRRfB\ntm31rkk6tm4NOY6DD65uua2t0NUF55yT7Pgnn4Tvfre6dXC1MXYsPP548s96/vzQzVWuMWPgD38o\n/3yAn/88PIbe1lZZOfnWrIEf/xg+n9H5RYoGE0lPA1Zsv5mVMF9n7/TCC+Gb1Le/Xe+apOeYFLJn\ngwbBI48kn0G4Tx847bTq18Olr7UV7ror+eO6110HQ4aU/37V6Ob6znfCGJlqB5Nf/Qpuu60Jgwnw\nn6M/r4z+vCv68+PpVae5rF4dvrkn/dbl9vJuq97hgAPgzDNr936VBpN4poUjjqhenWJz54a6mWVz\ncbeiwcTMVgFIOs3Mcif1vk7SQuC6tCuXdWvWhD5a51xjGDkytHjfeae82aWXLg3TuKxbF8o59NDq\n1a2jA95+OzydWM1yayVJAl6S3pPzy+yE5/V6q1f7WuTONZI+fUJAWbu2vPPnzoX3vjc8PVjJ2jv5\n4vEzkyZl9wGBJEHhk8C3Jb0o6UXC9CqXpVqrJuHBxLnGU8lU9PHTZOXOIVdMR0d4sGDcuCYOJmb2\npJlNAaYAU8xsqpktTL9q2efBxLnGU+k6KrNmlT6wNmm5WR4H02MwkXS4pDuAu83sDUkTJX2yBnXL\nPA8mzjWecm/Yr70WpvmZNClMIVTJ2jv54hZPlsfBJOnm+j/Ag8CR0e/PAp9Lq0LNxBPwzjWecoNJ\nPJVLnz6Vrb2Tb8cOWLgwlJ3l1SCTBJPhZvYTYDeAme0EEq9AIOkMScslPSvp2gL7h0q6V9JiSZ2S\nJubsGyLpHknLJC2RdEq0/QZJayQtjH7OSFqfWtm8OazTMHRovWvinMtV7rf//Klcyl0uId/ixXDU\nUWH8TFN3cwFbJR1KNIBR0kzgjSSFS2oBvgWcDkwCLpR0XN5h1wNdUV7mYuAbOfu+DvzazI4n5GyW\n5ey7xcymRz8PJKlPLcVdXFl8Xty5Zlbut//8qVyqlYTPLbfZg8nfAfcD75L0B+D/AlcnLL8VeM7M\nVpnZDuBu9q6PEpsIPAZgZiuAo6I15wcD/8nMvh/t22lmm3POa+jbtOdLnGtM5dywd+wI0/bkTuVS\nrSR8nHyH0GpauzYMjsyaJMFkCfA+YDZwBaGFsTxh+aOA3I9tTbQt12LgAgBJrcBYYDQwHnhN0vej\nrqzbJQ3IOe8qSYskfVdSBRMspMODiXONacSI0A399tvJz3nqqfDYbm639YQJ5a29ky+3ZTJgQJhO\n6LXXKiuzHpJM9NhhZtMJQQWAaAT89CrV4Sbg61GZTwNdhJxM3+g9rjSzBZJuJYy6v4Ew1uVLZmaS\nvgzcQhgPs585c+bsed3W1kZbtSfUKcKT7841ppYWOPLI8H806dxyhWYrzl175/zzy6vL2rVhwtR3\nv3vvtrjlNGJEeWWWo729nfb29orK6G6ix5GEVsQASdPY2600GDgoYflrCS2N2Oho2x5mtoWcQZCS\nXgBWAgOB1Wa2INr1U+Da6JxXc4r4N+CXxSqQG0xqafVqn1/KuUYVJ+GTBpOOjsKTicZJ+HKDSdzF\nlZtbjYPJjBnllVmO/C/aN954Y8lldNfNdTrwz4QAcAvwL9HP3xGS5knMB46RNE5SP+BjhPzLHtET\nW32j15cDvzGzN81sA7A6Z735U4Gl0XEjc4q4AHgmYX1qxru5nGtcpSbhi62jUmkSvlC5WU3CdzfR\n453AnZI+bGY/K6dwM9sl6SrgIULgusPMlkm6Iuy224Hjo/fZTehKy+2uuhr4QRRsVgKXRtu/Kmkq\n4XHlFwm5nIbiwcS5xlXKDfvll2HLljADeL7W1rBA3Pbt0L9/6fXo6ICvfKX8ujWSHnMmZvYzSWcT\nEu8H5mz/UpI3iB7bnZC37bac1535+3P2LQZOLrD9oiTvXS9mHkyca2RjxsAzCfszCnVFxQYNCkGm\nqyuMii/F22+HxP7JeXe40aPD9qxJMp3Kd4CPAp8h5E0+AoxLuV6ZtmlTWKeh2isQOueqo5Rv/z0t\nFTx7dnmPCC9cCMcfDwMHll+3RpLk0eDZUUvgT2Z2IzALKNDgczFvlTjX2Eq5YeeOAymk3JHw+SPq\ny6lbI0kSTN6K/twm6UhgB5DCOmPNw4OJc41t9OhkCfjt28N0J/ldUbniJLwVXeS8sGItnlGjwuJb\nuxJPWtUYkgST/ydpKPA1YCEh4f2jNCuVdR5MnGtsw4fDtm3hpzsLF4bBiYMGFT9m/Pgwe3AprQmz\n4i2e/v3DRJIbNiQvrxEkWc/kf5rZpuiJrnHAcWb2xfSrll0eTJxrbFKyCR97ypfEZZX6iPCLL4bz\nxhXJPidtOTWSosFE0gX5P8DZwKnRa1eEj353rvElyU0kCSZQehK+oyOcU2wi2CzmTbp7NPic6M8R\nhHm5Hot+fz8wF7g3xXplmrdMnGt8Pd2wzUIw+ed/7rmsWbPgmmuSv3ex5HvSujWi7gYtXgog6SFg\nopmti34/grBglivCg4lzja+nrqRVq8KfRx3Vc1kzZsDSpSEHc1CCyabmzoWPf7z4/iwGkyQJ+DFx\nIIlsYN/5tlwOs/AP1IOJc42tpxt2d4MV8w0YAJMnw4IFPR/75puwYgVM72aq3GYNJo9KelDSJZIu\nAX4FPJJutbLrtdfCN5Mk306cc/XT0w07ab4kljQJP38+TJnS/fQrTZWAj5nZVcBthJUOpwC3m9ln\n0q5YVnny3blsqFcwSVJuFlsmSdYzwczuxRPuiXi+xLls6O6GvXUrLF/efVdUvlmz4G//NnR1d9c1\n1tEBl16eP7BeAAAWWUlEQVRafD+E9VY2bAjjVw5IdJeuv+4eDf599OcWSZtzfrZI2lzsvN7Og4lz\n2TB0aLhZby5wN5s/H048EQ48cP99xYweHXInzz9f/JjuBivm6tsXDjssjITPiqLBxMzeG/15sJkN\nzvk52MwG166K2eLBxLlskIqva5Lkhl9IT+vCP/tsmAD2yCN7LitrXV3dtUwO6e6nlpXMEg8mzmVH\nsRt2qfmSWE95k1LKzVoSvrveuCcBY+9yvbkMODqVGmWcJ+Cdy45CwSTuirrttsLndGf2bLjjjuL7\nSwkmWWuZdDdocXwtK9IsvGXiXHYUumE/91xYYyRJV1S+KVNg5cqQhxlcIBnQ0QF//dfJ6/bSS6XX\noV6SjDNB0jBJrZL+LP5Ju2JZtHs3rF3rLRPnsqJQV1K5XVwA/fqFJ8Dmzdt/36ZNYVT9iScmKytr\nLZMkKy1+Cvgt8CBwY/TnnHSrlU2vvAJDhpT2BIhzrn4K3bDLTb7HiiXh580L06707Vt+3RpZkpbJ\nZwnrsK8ys/cD04BNqdYqo7yLy7lsKXTDrqRlAsWT8KWWm7UEfJJg8raZvQ0gqb+ZLQcmpFutbPLk\nu3PZEgeTeJXEN96AF14IuY9yzZoFnZ2h2ztXqcHkiCPC9EzvvFN+XWopSTBZE620+AvgYUn3AavS\nrVY2ecvEuWwZPBj69An5DCi9K6qQESPCSo7Llu3dtmsXPPEEzJyZvJw+fWDkSHj55fLrUktJ5uY6\nP1ppcQ7wReAO4ENpVyyLPJg4lz253UmVdnHF8ru6liwJgWH48NLKyVLepLtBi7+W9AlJe1Y/NrPf\nmNn9ZpaRhldteTBxLntyb9iVJt9j+Un4csttimBCmCn4bOAFST+RdL6kfjWqVyZ5MHEue+Ib9q5d\nIddRjWCS3zIpt8WTpSR8d3Nz3WdmFwLjgJ8BFwEvSfq+pNNqVcEs8QS8c9kTB5OlS+Hww8MEi5Wa\nPDnkOl5/PfxebjBplpYJAGa2zcx+bGbnA38OTAUeSL1mGbNrF6xfD6NG1bsmzrlSxDfsanVxQUie\nt7aGls6rr4afiRPLr1sW9DhTvqTDgb8EPgYcAfwEuCTdamXP+vVwyCFhBKxzLjvirqRqJd9jcVfX\nrl1wyinQkmi+kX1lKZh0l4C/XNJjwELg3cDfm9nRZnadmS1O+gaSzpC0XNKzkq4tsH+opHslLZbU\nKWlizr4hku6RtEzSEkmnRNuHSXpI0opoSeEhJV11Cjxf4lw2pdEygb1J+ErKbYpgAswCvgKMMbOr\nzSzBgpT7ktQCfAs4HZgEXCjpuLzDrge6zGwKcDHwjZx9Xwd+bWbHE5YMjp/cvg54xMwmAI8BXyi1\nbtXmwcS5bBozBl58MfQuTJpUvXJnzgyLbP32t+W3eEaMCAMp3367evVKS3cJ+MvM7GEz2zOOU9Kc\nEstvBZ4zs1VmtgO4Gzgv75iJhICAma0AjpJ0mKTBwH8ys+9H+3aaWbwm2nnAndHrO2mAcS+efHcu\nmwYOhEGDQldUnz7VK3fYMBg7NuRNTjmlvDJaWsLsxWvXVq9eaSm1F+/cEo8fBeQ20tZE23ItBi4A\nkNQKjAVGA+OB16KnxxZKul3SgOicEWa2AcDM1gMjSqzXPnbsCNMWVMJbJs5l15gx1e3iis2aFRLv\nQyroiM9KV1epS9UXWiirUjcBX5e0EHga6AJ2AX2B6cCVZrZA0q2E7q0bCtTDihU+Z86cPa/b2tpo\na2vb75i77oJHH4Uf/KD8i1i9uvxvH865+nr/++HMM6tf7jnnhNZJJWoRTNrb22lvb6+oDJkVvQ/v\nf7DUktvtleD4mcAcMzsj+v06wMzs5m7OeQE4ARgIdJjZ0dH29wLXmtk5kpYBbWa2QdJI4PEor5Jf\nliW5vmXL4KyzwgRv5Zo5E265pbpPgzjn3HXXhTnErr++du8pCTMrqfGQZD2Tr0oaLKkvYaLHVyV9\nImH584FjJI2LRs9/DLg/r/whUdlIuhz4jZm9GXVjrZZ0bHToqcDS6PX97H08+WLgvoT1KWjChJDk\nWreu/DK8m8s5l4asjIJPkjP58yjx/Z+BF4FjgL9PUriZ7QKuAh4ClgB3m9kySVdI+nR02PHAM1Fr\n43TC+imxq4EfSFpEeJrrn6LtNwOnSVpBCDI3JalPMS0txRe0SWLHjjAo6YgjKqmFc87tr5lyJvEx\nZwP3mNkbUvLWj5k9QN76J2Z2W87rzvz9OfsWExbmyt++Efhg4kokEAeTCy4o/dx168IjfAeUmoFy\nzrkeZCWYJGmZ/D9Jy4EZwKOSDgMy8NRzaYqtjpaEd3E559LSNMHEzK4DZgMnRWNFtrL/WJHMa22F\nRYtg+/bSz/Vg4pxLy/DhsG1b+GlkSRLwHwF2mNkuSf8A/DtwZOo1q7FBg+DYY6Grq/RzPZg459Ii\nhQlkGz0Jn6Sb64tmtiV6NPeDhJUW/zXdatVHuV1dPvrdOZemLHR1JQkmu6I/zwZuN7NfAU05N265\nT3R5y8Q5l6ZmCSZrJd0GfBT4taT+Cc/LnLhlUsI4TsCDiXMuXc0STP4SeBA43cw2AYeQcJxJ1owf\nDzt3wksvlXaeBxPnXJqaIpiY2Tbgj8Dpkq4iTLL4UOo1qwMptE5K6eravh02bgzLfTrnXBqyMAo+\nydNcnwV+QJiZdwTw75I+k3bF6qXUJPzLL4eR79Wcuto553JloWWSZMz2J4FTzGwrgKSbgQ7gm2lW\nrF5mzYKf/CT58d7F5ZxLWxaCSZKcidj7RBfR6zSmom8IM2bA0qXJBwh5MHHOpW3YsDAH4JYt9a5J\ncUmCyfeBeZLmRCstdhLGmjSlAQNg8mRYsCDZ8R5MnHNpkxq/dZIkAX8LcCmwMfq51MxuTbti9VRK\n3sSDiXOuFho9Cd9tzkRSH2CJmR0HLKxNlepv9uyw+mISa9bAB6s6f7Fzzu0v0y2TaD2SFZIqXHgy\nW+KR8EkGL3rLxDlXC40eTJI8zTUMWCLpCcKMwQCY2bmp1arORo8OuZPnn4d3v7v7Yz2YOOdqYcwY\n6Oysdy2KSxJMvph6LRpQ3DrpLpi89VZ4uuKww2pXL+dc7zRmDNxzT71rUVzRbi5Jx0h6j5n9JveH\n8GhwA6eBqiNJEn7NmjA1dEtTzlTmnGskjZ6A7+42eCuwucD2N6J9TS1pMPGp551ztRDnTEqdiLZW\nugsmh5vZ0/kbo21HpVajBjFlCqxcCZsLhdOI50ucc7UyZEgYb/LGG/WuSWHdBZOh3ewbUO2KNJp+\n/WD6dJg3r/gxHkycc7XUyE90dRdMFki6PH+jpE8BT6ZXpcbR02JZHkycc7XUyMGku6e5Pgf8XNLH\n2Rs8TiKssnh+2hVrBLNnw792s0Dx6tVw9tm1q49zrndr5CR80WBiZhuA2ZLeD0yONv/KzB6rSc0a\nwKxZcPHFsHt34Se2PAHvnKulrLZMADCzx4HHa1CXhjNiBAwfDsuWwaRJ++/3bi7nXC2NGQO//W29\na1GYj5DoQbFHhLduhbffhkMPrX2dnHO9UyO3TDyY9KBYEn716tDFpaZd2cU512g8mGRYsZaJd3E5\n52otTsA34sDF1IOJpDMkLZf0rKRrC+wfKuleSYsldUqamLPvxWh7VzTRZLz9BklrJC2Mfs5Iq/6T\nJ4d13l9/fd/tnnx3ztXaoEHQvz9s3Fjvmuwv1WAiqQX4FnA6MAm4UNJxeYddD3SZ2RTgYuAbOft2\nA21mNs3MWvPOu8XMpkc/D6R0CfTpA62t+8/W6S0T51w9NGpXV9otk1bgOTNbZWY7gLuB8/KOmQg8\nBmBmK4CjJMXz8KqbOtYsW1Goq8uDiXOuHnprMBkF5F72mmhbrsXABQCSWoGxQNyBZMDDkuYXGI1/\nlaRFkr4raUj1q75XoSS8BxPnXD301mCSxE3AMEkLgSuBLsI09wDvMbPpwFnAlZLeG23/NnC0mU0F\n1gO3pFnBmTNh/nzYuXPvNg8mzrl6GDOmMUfBJ1kcqxJrCS2N2Oho2x5mtgW4LP5d0gvAymjfuujP\nVyX9nNBt9nszezWniH8DflmsAnPmzNnzuq2tjba2tpIvYtgwGDsWnnoqTP4InoB3ztXH6NHwyCPV\nLbO9vZ329vaKypCl+IyZpD7ACuBUYB3wBHChmS3LOWYIsM3MdkRdWe8xs0skHQS0mNmbkgYCDwE3\nmtlDkkaa2fro/GuAk83srwq8v1Xr+j71KZg2Da68MkxLf+SRYZVFH2finKulxx6DL30JKrz3d0sS\nZlbS3S3Vbi4z2wVcRQgES4C7zWyZpCskfTo67HjgGUnLCE99fTbafjjwe0ldQCfwSzN7KNr3VUlP\nSVoEvA+4Js3rgH2T8HEXlwcS51ytNWrOJNWWSb1Vs2WybFmYIXjlSnjgAbjlFnjooZ7Pc865anrr\nrdD1vm1bekuGN1zLpJlMmACbNsG6dZ58d87Vz4ABcPDB8Npr9a7JvjyYJNTSsvcRYU++O+fqafTo\nxuvq8mBSgjiYeMvEOVdPjZg38WBSgjgJ78HEOVdPjRhM0h5n0lRaW2HRIjj8cA8mzrn6acRg4i2T\nEgwaBMceCy+84MHEOVc/jTgK3oNJiWbPhsGDw9MUzjlXD42YgPdurhLNmgW/+U29a+Gc683GjIGl\nSyFntqi682BSonPPDQOGnHOuXsaNg+uug61b612TvXwEvHPOuX34CHjnnHN14cHEOedcxTyYOOec\nq5gHE+eccxXzYOKcc65iHkycc85VzIOJc865inkwcc45VzEPJs455yrmwcQ551zFPJg455yrmAcT\n55xzFfNg4pxzrmIeTJxzzlXMg4lzzrmKeTBxzjlXMQ8mzjnnKubBxDnnXMVSDyaSzpC0XNKzkq4t\nsH+opHslLZbUKWlizr4Xo+1dkp7I2T5M0kOSVkh6UNKQtK/DOedccakGE0ktwLeA04FJwIWSjss7\n7Hqgy8ymABcD38jZtxtoM7NpZtaas/064BEzmwA8BnwhrWtoZO3t7fWuQqqa+fqa+drAr683Srtl\n0go8Z2arzGwHcDdwXt4xEwkBATNbARwl6bBon4rU8Tzgzuj1ncCHql3xLGj2f9DNfH3NfG3g19cb\npR1MRgGrc35fE23LtRi4AEBSKzAWGB3tM+BhSfMlXZ5zzggz2wBgZuuBESnU3TnnXEIH1LsCwE3A\n1yUtBJ4GuoBd0b73mNm6qKXysKRlZvb7AmVYjerqnHOuAJmldx+WNBOYY2ZnRL9fB5iZ3dzNOS8A\nJ5jZm3nbbwC2mNktkpYRcikbJI0EHjez4wuU5UHGOefKYGYq5fi0WybzgWMkjQPWAR8DLsw9IHoS\na5uZ7Yi6sn5jZm9KOghoiV4PBP4cuDE67X7gEuBmQtL+vkJvXupfhnPOufKkGkzMbJekq4CHCPmZ\nO8xsmaQrwm67HTgeuFPSbmAJ8Mno9MOBn0etiwOAH5jZQ9G+m4GfSLoMWAX8ZZrX4ZxzrnupdnM5\n55zrHZpyBHxPAyWzrthgzqySdIekDZKeytnWNANTi1zfDZLWSFoY/ZxRzzpWQtJoSY9JWiLpaUlX\nR9sz/xkWuLbPRNub4vOT1F/SvOhe8nSUmy7rs2u6lkk0UPJZ4FTgZULe5mNmtryuFasiSSuBGWb2\np3rXpRokvRd4E/i/ZnZitO1m4HUz+2r0hWCYmV1Xz3qWq8j17XmgpK6Vq4LoIZiRZrZI0iDgScJY\nsEvJ+GfYzbV9lOb5/A4ys22S+gB/AK4GPkyJn10ztkySDJTMumKDOTMpetw7PzA2zcDUItcH4XPM\nPDNbb2aLotdvAssIY8Uy/xkWubZ4rFyzfH7bopf9Cflpo4zPrmluSDmSDJTMumKDOZtJbxiYepWk\nRZK+m8UuoEIkHQVMBTqBw5vpM8y5tnnRpqb4/CS1SOoC1gMPm9l8yvjsmjGY9AbvMbPpwFnAlVE3\nSrNrrv5Y+DZwtJlNJfwnbobukkHAT4HPRt/i8z+zzH6GBa6taT4/M9ttZtMIrclWSZMo47NrxmCy\nljAlS2x0tK1pmNm66M9XgZ8TuvaazQZJh8OefutX6lyfqjKzV21vwvLfgJPrWZ9KSTqAcLO9y8zi\ncV9N8RkWurZm+/wAzGwz0A6cQRmfXTMGkz0DJSX1IwyUvL/OdaoaSQdF35LIGcz5TH1rVRVi3z7o\neGAqdDMwNUP2ub7oP2jsArL/GX4PWGpmX8/Z1iyf4X7X1iyfn6ThcRedpAHAaYS8UMmfXdM9zQXh\n0WDg6+wdKHlTnatUNZLGE1ojuYM5M319kn4ItAGHAhuAG4BfAPcAY4gGpprZpnrVsRJFru/9hP73\n3cCLwBVxH3XWSHoP8FvC3HoW/VwPPAH8hAx/ht1c21/RBJ+fpBMICfaW6OfHZvaPkg6hxM+uKYOJ\nc8652mrGbi7nnHM15sHEOedcxTyYOOecq5gHE+eccxXzYOKcc65iHkycc85VzIOJy7RoevDT8rZ9\nVtL/7uG8LSnXa7ikTklPRmMVcvc9Lml69Hp8tFTCaQXK+Fo0LXjRZa57qMP7JP0y5/cvS/q1pL6S\n2iXNz9k3Q9LjOeftlnR2zv5fSvqzcurhegcPJi7rfkjeUtCEWQ9+2MN5aQ+w+iDwlJnNMLM/FDpA\n0mjgP4BrzOzhAodcDpxoZonW5ImmEM9n0b5/AGYBH4pm0zbgMEmn5x8bWQP89yTv6xx4MHHZ9zPg\nrGj+JCSNA44wsz9IGijpEUkLFBYTOzf/5ALf3r8p6aLo9fT4G7yk/4jnKso7f5ykR6PyH1ZYTGkK\nYWnp86KFk/oXqPeRwIPAF8zsVwXKvQ8YBDwp6SM577Mofp/ouO9L+ldJndF7FihKfwecDpxjZu/k\n7Psa8A8F/1ZhMfCGpFOL7HduHx5MXKZFC4Q9AZwZbfoYYRoIgLcJ38RPAj4A/EuxYvI3RMHpm8CH\nzexk4PvAPxU495vA981sCqE19E0zWwz8D8LUFNPNbHuB8+6Mjv15kes6D9gWnX9PzvtMjd8n5/BR\nZjbTzP5rgaLeA1wBnJmzbkV8zR3AdknvK1QF4B+BLxaqn3P5PJi4ZnA3IYgQ/fmj6LWAr0haDDwC\nHCkp6ZoaE4DJhHVjughdPkcWOG5WzvvdRbh5J/Ew8AlJB3ZzTO7El929zz3dlPF8VM6fFym7aMCI\nFvWy/JyPc4V4MHHN4D7gVEnTgAFm1hVt/zgwHJgWrdfwCpB/897Jvv8P4v0CnolaBtPMbIqZncn+\nys29fJUww/VPFZaaLsSKvM63tZt96wnr3twqqW2/NzB7nHDNM4uc/0+ErjCfxM91y4OJyzwz20pY\nh+F77P32DjAEeMXMdkt6PzAuZ1/8zXwVMDF6wmkoEOcIVhAS1DMhdHtJmljg7eey9wGATwC/K6He\n1wBvRPUuJLdlUsn7PE+YJv3fJZ1Y4JB/BP5bkXMfBoYBhc5zbg8PJq5Z/Ihww8sNJj8ATo66uT5B\nWKchZgBmtoaQY3mG0F22MNq+A/gL4GZJi4AuQldTvquBS6NjPg58NkFdc7/lXwKMLPL4b+5xxd4n\nUYvBzBYAlwL3R8sYWM6+/yC02oqV9Y+EqcidK8qnoHfOOVcxb5k455yrmAcT55xzFfNg4pxzrmIe\nTJxzzlXMg4lzzrmKeTBxzjlXMQ8mzjnnKubBxDnnXMX+PygS5di7aan1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x87e1a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# plot the value of K for KNN (x-axis) versus the cross-validated accuracy (y-axis)\n",
    "plt.plot(k_range, k_scores)\n",
    "plt.xlabel('Value of K for KNN')\n",
    "plt.ylabel('Cross-Validated Accuracy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation example: model selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Goal:** Compare the best KNN model with logistic regression on the iris dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.98\n"
     ]
    }
   ],
   "source": [
    "# 10-fold cross-validation with the best KNN model\n",
    "knn = KNeighborsClassifier(n_neighbors=20)\n",
    "print(cross_val_score(knn, X, y, cv=10, scoring='accuracy').mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.953333333333\n"
     ]
    }
   ],
   "source": [
    "# 10-fold cross-validation with logistic regression\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "logreg = LogisticRegression()\n",
    "print(cross_val_score(logreg, X, y, cv=10, scoring='accuracy').mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation example: feature selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Goal**: Select whether the Newspaper feature should be included in the linear regression model on the advertising dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the advertising dataset\n",
    "data = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv', index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a Python list of three feature names\n",
    "feature_cols = ['TV', 'Radio', 'Newspaper']\n",
    "\n",
    "# use the list to select a subset of the DataFrame (X)\n",
    "X = data[feature_cols]\n",
    "\n",
    "# select the Sales column as the response (y)\n",
    "y = data.Sales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-3.56038438 -3.29767522 -2.08943356 -2.82474283 -1.3027754  -1.74163618\n",
      " -8.17338214 -2.11409746 -3.04273109 -2.45281793]\n"
     ]
    }
   ],
   "source": [
    "# 10-fold cross-validation with all three features\n",
    "lm = LinearRegression()\n",
    "scores = cross_val_score(lm, X, y, cv=10, scoring='mean_squared_error')\n",
    "print(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.56038438  3.29767522  2.08943356  2.82474283  1.3027754   1.74163618\n",
      "  8.17338214  2.11409746  3.04273109  2.45281793]\n"
     ]
    }
   ],
   "source": [
    "# fix the sign of MSE scores\n",
    "mse_scores = -scores\n",
    "print(mse_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.88689808  1.81595022  1.44548731  1.68069713  1.14139187  1.31971064\n",
      "  2.85891276  1.45399362  1.7443426   1.56614748]\n"
     ]
    }
   ],
   "source": [
    "# convert from MSE to RMSE\n",
    "rmse_scores = np.sqrt(mse_scores)\n",
    "print(rmse_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.69135317081\n"
     ]
    }
   ],
   "source": [
    "# calculate the average RMSE\n",
    "print(rmse_scores.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.67967484191\n"
     ]
    }
   ],
   "source": [
    "# 10-fold cross-validation with two features (excluding Newspaper)\n",
    "feature_cols = ['TV', 'Radio']\n",
    "X = data[feature_cols]\n",
    "print(np.sqrt(-cross_val_score(lm, X, y, cv=10, scoring='mean_squared_error')).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improvements to cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Repeated cross-validation**\n",
    "\n",
    "- Repeat cross-validation multiple times (with **different random splits** of the data) and average the results\n",
    "- More reliable estimate of out-of-sample performance by **reducing the variance** associated with a single trial of cross-validation\n",
    "\n",
    "**Creating a hold-out set**\n",
    "\n",
    "- \"Hold out\" a portion of the data **before** beginning the model building process\n",
    "- Locate the best model using cross-validation on the remaining data, and test it **using the hold-out set**\n",
    "- More reliable estimate of out-of-sample performance since hold-out set is **truly out-of-sample**\n",
    "\n",
    "**Feature engineering and selection within cross-validation iterations**\n",
    "\n",
    "- Normally, feature engineering and selection occurs **before** cross-validation\n",
    "- Instead, perform all feature engineering and selection **within each cross-validation iteration**\n",
    "- More reliable estimate of out-of-sample performance since it **better mimics** the application of the model to out-of-sample data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources\n",
    "\n",
    "- scikit-learn documentation: [Cross-validation](http://scikit-learn.org/stable/modules/cross_validation.html), [Model evaluation](http://scikit-learn.org/stable/modules/model_evaluation.html)\n",
    "- scikit-learn issue on GitHub: [MSE is negative when returned by cross_val_score](https://github.com/scikit-learn/scikit-learn/issues/2439)\n",
    "- Section 5.1 of [An Introduction to Statistical Learning](http://www-bcf.usc.edu/~gareth/ISL/) (11 pages) and related videos: [K-fold and leave-one-out cross-validation](https://www.youtube.com/watch?v=nZAM5OXrktY) (14 minutes), [Cross-validation the right and wrong ways](https://www.youtube.com/watch?v=S06JpVoNaA0) (10 minutes)\n",
    "- Scott Fortmann-Roe: [Accurately Measuring Model Prediction Error](http://scott.fortmann-roe.com/docs/MeasuringError.html)\n",
    "- Machine Learning Mastery: [An Introduction to Feature Selection](http://machinelearningmastery.com/an-introduction-to-feature-selection/)\n",
    "- Harvard CS109: [Cross-Validation: The Right and Wrong Way](https://github.com/cs109/content/blob/master/lec_10_cross_val.ipynb)\n",
    "- Journal of Cheminformatics: [Cross-validation pitfalls when selecting and assessing regression and classification models](http://www.jcheminf.com/content/pdf/1758-2946-6-10.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comments or Questions?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        width: 90%;\n",
       "/*        margin-left:auto;*/\n",
       "/*        margin-right:auto;*/\n",
       "    }\n",
       "    ul {\n",
       "        line-height: 145%;\n",
       "        font-size: 90%;\n",
       "    }\n",
       "    li {\n",
       "        margin-bottom: 1em;\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: Helvetica, serif;\n",
       "    }\n",
       "    h4{\n",
       "        margin-top: 12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "    div.text_cell_render{\n",
       "        font-family: Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 145%;\n",
       "        font-size: 130%;\n",
       "        width: 90%;\n",
       "        margin-left:auto;\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\", source-code-pro,Consolas, monospace;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #4057A1;\n",
       "        font-style: italic;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}