07_cross_validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cross-validation for parameter tuning, model selection, and feature selection ([video #7](https://www.youtube.com/watch?v=6dbrR-WymjI&list=PL5-da3qGB5ICeMbQuqbbCOQWcS6OYBr5A&index=7))\n",
    "\n",
    "Created by [Data School](http://www.dataschool.io/). Watch all 9 videos on [YouTube](https://www.youtube.com/playlist?list=PL5-da3qGB5ICeMbQuqbbCOQWcS6OYBr5A). Download the notebooks from [GitHub](https://github.com/justmarkham/scikit-learn-videos).\n",
    "\n",
    "**Note:** This notebook uses Python 3.6 and scikit-learn 0.19.1. The original notebook (shown in the video) used Python 2.7 and scikit-learn 0.16, and can be downloaded from the [archive branch](https://github.com/justmarkham/scikit-learn-videos/tree/archive)."
   ]
  },
  {
   "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.model_selection 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.9736842105263158\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.model_selection import KFold\n",
    "kf = KFold(n_splits=5, shuffle=False).split(range(25))\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], str(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.model_selection 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.9666666666666668\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.96, 0.9533333333333334, 0.9666666666666666, 0.9666666666666666, 0.9666666666666668, 0.9666666666666668, 0.9666666666666668, 0.9666666666666668, 0.9733333333333334, 0.9666666666666668, 0.9666666666666668, 0.9733333333333334, 0.9800000000000001, 0.9733333333333334, 0.9733333333333334, 0.9733333333333334, 0.9733333333333334, 0.9800000000000001, 0.9733333333333334, 0.9800000000000001, 0.9666666666666666, 0.9666666666666666, 0.9733333333333334, 0.96, 0.9666666666666666, 0.96, 0.9666666666666666, 0.9533333333333334, 0.9533333333333334, 0.9533333333333334]\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": [
       "Text(0,0.5,'Cross-Validated Accuracy')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xucm2d54P3fNUeNZ0aSD+OR7ZHjJOTk2J4xmJRzQlog9AAJAUqWtsCypdstvO0WWKD0ZWlKXlqgpe9uWbq0UMiWEiCcUhoaaEjoiUIcLNlOnBjHSayxZ+zxQZqTNTOSrv3jeTSWZUnz6DQaaa7v5zOfkZ6T7sca69J9um5RVYwxxphKtTW6AMYYY5qbBRJjjDFVsUBijDGmKhZIjDHGVMUCiTHGmKpYIDHGGFMVCyTGGGOqYoHEGGNMVSyQGGOMqUpHowuwHDZs2KDbtm1rdDGMMaapPProo6dVdWCp41ZFINm2bRt79+5tdDGMMaapiMizXo6zpi1jjDFVsUBijDGmKhZIjDHGVMUCiTHGmKpYIDHGGFOVugYSEblFRJ4UkSMi8v4C+y8TkQdFZL+IPCwiQzn7PiYij4nIIRH5HyIi7vbnicgB95qL240xxjRG3QKJiLQDnwJeDWwH7hCR7XmHfQK4W1V3AXcCH3XPfRHwYmAXsAN4PnCje86ngXcAV7k/t9TrHowxxiytnjWSG4AjqnpUVeeBe4DX5h2zHXjQffxQzn4FfEAX0A10AidFZBPgV9UfqrNG8N3ArXW8B7MKPfrsOSKxeKOLUTNzqTRf+vEx0pnaLaudXEhzz4+PkanhNU3zqmcg2QLEcp6PuttyRYHb3ce3Af0isl5Vf4gTWMbcnwdU9ZB7/ugS1wRARN4hIntFZO/ExETVN2NWj/d9bT8f/MaBRhejZu4/MMYHvn6Af/5p7f4ffHv/GO//+gH2xc7V7JqmedUzkBTqu8j/+vIe4EYR2YfTdHUcSInIc4DrgCGcQHGziLzM4zWdjaqfUdU9qrpnYGDJGf7GADCZXOCpiWmeGJ/i/Hy60cWpicgxp3ZVy1pWxA0gJ+LJml3TNK96BpJRIJzzfAg4kXuAqp5Q1dep6m7gg+62BE7t5N9VdVpVp4HvAC9wrzlU6prGVOPAaAJVSGeUx04kGl2cmoiMOvcRrWEgicaca56ctEBi6htIHgGuEpHLRaQLeBNwX+4BIrJBRLJl+ADwOffxMZyaSoeIdOLUVg6p6hgwJSIvcEdr/RrwrTreg1llcr+1t0I/yVwqzaETkwBERxM4XYvVSS6kOTTmXHMsYYHE1DGQqGoKeCfwAHAI+IqqPiYid4rIa9zDbgKeFJHDwCBwl7v9XuAp4ABOP0pUVf/O3febwF8BR9xjvlOvezCrTyQW5/INvWwJ9rREIDk0NsV8OsNLr9rA2Zl5YmfPV33Nx05MknI72cetRmKoc/ZfVb0fuD9v24dyHt+LEzTyz0sDv1HkmntxhgQbU1OqSiQW5yXP2cB8KkN0tPkDSbY56y0v3MY///Q0kdE4W9evqck1r9jQy0mrkRhsZrsxi8Ynk0xMzTE8FGA4HCB29jxnpucaXayqRGNxBvq7ufGaAXydbTXpJ4mOxgn5fewaCliNxAAWSIxZlB3dNBwOMjwUBGj6WkkkFmd4KEhnexs7Ngdq0lwXicUZDgcIBXo4NTlnc0mMBRJjsiKjcTrbhe2b/ewcCtAmEIk178itxOwCR0/PsHurExRHwkEOHk+wkM5UfM1zM/M8e2aWkfBaQv5u5tMZzs7O16rIpklZIDHGFTkWZ/smP90d7azp6uDqwf6m7nDP1qaytavhcJC5VIYnx6cqvmYke81wgFDAB8C49ZOsehZIjMGZN3LgeIKRcHBx2+6tQaKxeE2GzDZCtj9kVzgAsHhv1QTHaCyOCOzcEmDQ7wQSm0tiLJAYAxw5Nc3sfJrhnEAyPBQkcX6BZ8/MNrBklYuOxrlyoBe/rxOAobU9rO/tqqrDPRqLc9XGPvp9nRdqJBZIVj0LJMZwIeXHRYGkBt/gGyU7lDn3fkSE4XCw4vtZvKbbVDbQ102bYEOAjQUSY8DpVPf7Orh8fe/itqsH+1nT1d6UgeR4/Dynp+fZnRNIwGneOjIxzVRyoexrxs6e59zswmJw6mhvY6C/22okxgKJMeA02QyHg7S1XcgL2t4m7NhSmyGzyy1b5uG8QDIcDqIKB46XPxot29Ge248U8vssTYqxQGLM+fk0T56cWmyyyTUSDvL4iUnmU5UPmW2EaCxOV0cb14b8F20fHnI63isJjpFjcbo72rgm1L+4bdDvs852Y4HEmIMnEqQzetE37ayRcJD5dIYnxicbULLKRWMJrt/sp6vj4v/iwTVdXL6ht6IO9+honJ1bAnS2X7hmKOCz4b/GAokx0SLNQLnbmql5K5XOcOB4omANC5xaSbTMiZYL6QwHjycu+Tca9PuYTKZaZu0WUxkLJGbV2xeLsyXYw0B/9yX7Ngd8DPR3N1UgOXxymvML6cUZ7flGwkHGJ5Nl1SSeHJ9iLpW5JJBssiHABgskxhCNxQs2a4E7ZHYoWNNFoeotf0Z7vkpqWdlj80eBhdxJiWOJ6tPTm+ZlgcSsaqen5xg9d55hd/Z3ISPhAE9NzJA4X/6Q2UaIHIsTXNPJZUXSxV+3yU9nu5QdSNb1djG0tuei7YMBm91uLJCYVS5b0xgJry16THbfgdHmSOAYHXUmDTqLiF7K19nO9k3+smpZ0Vic4aHAJdfM1kjGE82dbt9UxwKJWdWisThtAju2+Ises3NxyOy55SpWxWbmUhw+OVVw4ECu4XCQA8ed0WpLmUoucGRiumCw7e3uoL+7w2okq5wFErOq7YvF3RnsxRcLDfR0csVAb1OklD9wPEFGnea4UoaHgkzPpXhqYnrpa44mUKVo89+gDQFe9SyQmFVLVYnG4kVHN+UacXNUrfRMwItDmYt0tGeNbPXe4V5oRnuuTQEfY1YjWdUskJhV65kzs0wmU0t+6ILzIXp6eo4TK/ybdyQWJ7yuh/V9lw5lznX5+l76fR2e+kmisTjb1q8huKar4P5Bv88SN65yFkjMqpXt8xjxWCMBVvwwYGcoc/GBA1ltbbJYy1pKfhbhfCG/j4npOU/9LaY1WSAxq1Y0lmBNVztXbexf8thrQ3662ttWdCA5NZnkRCK5mE9rKcNDQZ4YnyK5UHxW+ngiycnJuaLNWuD0kaQzyulpG7m1WlkgMavWvlicHVsCtLcVHiabq6ujje2b/exbwYEkEivdl5FvOBwknVEeO1F8EEGxLMK5LgwBtuat1coCiVmV5lJpDp2YvGSmdikj4SAHRhOk0iszE3B0NL6Y+t6L7CisfceKB8dILE5nu7B9U/Hh0dk0KZZOfvWyQGJWpSfGpphPX5o7qpSRcJDzC2l+emrpIbONEInFuTbUj6+z3dPxG/t9bAn2EC0x0TIai3PdJn/Ja9ra7cYCiVmVvDTZ5BtewR3umYyyP3Zpdt6lDIcDRSdapjPK/tH4kqPa1vd20dkulrhxFbNAYlalaCzOQH83m91mGS+2rV9DoKdzMSniSnL09AxTcynP/SNZI+EgsbPnOVOgo/ypiWlm5tNLXrOtTdjYb0OAVzMLJGZViiyRj6oQEWE4HCzZp9Ao5Xa0Z2VrG/sLNG+VU2sb9Nva7atZXQOJiNwiIk+KyBEReX+B/ZeJyIMisl9EHhaRIXf7y0UkkvOTFJFb3X2fF5Gnc/aN1PMeTOtJzC5wdGLG04z2fCPhIIdPTjE7n6pDySoXjcXp6+7gyoG+ss7bORSgTSg4Gi0Si9Pv6+CKDb1LXicU8FkgWcXqFkhEpB34FPBqYDtwh4hszzvsE8DdqroLuBP4KICqPqSqI6o6AtwMzALfzTnvvdn9qhqp1z2Y1rT/uLc0IoWMhANkFA4eX1lL72aXwfUylDnXmq4Orh7sL9jv42T8DdLm4Zohfw/jieSKTyFj6mPJQCIivygilQScG4AjqnpUVeeBe4DX5h2zHXjQffxQgf0Arwe+o6qzFZTBmEtkPzR3epy4lysbfFZSJuDkQppDY5Nld7RnjYSDREcvziOWXEjzxPhUyXVacoUC3czOp5maW1k1NbM8vASINwE/FZGPich1ZVx7CxDLeT7qbssVBW53H98G9IvI+gKv/6W8bXe5zWGfFJGCSYVE5B0isldE9k5MTJRRbNPqIrE4Vw70EujpLPvc9X3dhNf1lL3meT09PjbJQlrL7h/JGgkHic8u8OyZC9/VDrop5r2kW4GcIcDW4b4qLRlIVPVXgN3AU8Bfi8gP3Q/ppfJKFKoP59d73wPcKCL7gBuB48DiVxoR2QTsBB7IOecDwLXA84F1wPuKlPszqrpHVfcMDAwsUVSzWqgqkQqGyeYaHvKWo2q5RI5V1tGetTisOWc02mJHu8da2+LsdusnWZU8NVmp6iTwNZzmqU04tYefiMi7Spw2CoRzng8BJ/Kue0JVX6equ4EPuttyv+q9EfiGqi7knDOmjjngr3Ga0Izx5Hj8PKenS+eOWspIOMjx+HlOTa2MD83oaJxBfzehMoYy57pqYx89ne0XjUaLxOJsDvjY6Pd2zexrW5qU1clLH8kvicg3gO8DncANqvpqYBinRlHMI8BVInK5iHThNFHdl3ftDTn9Lx8APpd3jTvIa9ZyaymIM27zVuDgUvdgTFa2SaraQAKwf4U0bzkZfyu/n472NnYOBS6qkURH456yImcNWr6tVc1LjeQNwCdVdZeqflxVTwG4nd//sdhJqpoC3onTLHUI+IqqPiYid4rIa9zDbgKeFJHDwCBwV/Z8EdmGU6P5Qd6lvygiB4ADwAbgIx7uwRjA+YDs6mjj2lDx3FFLuX6zMzpqJTRvnZuZ55kzs1U11YETHB87Mcl8KsOZ6TliZ8+XNarN19nO2jWd1rS1ShVfX/SC/w6MZZ+ISA8wqKrPqOqDxU8DVb0fuD9v24dyHt8L3Fvk3Ge4tHMeVb3ZQ5mNKShyLM71m/10dVQ+8r2nq51rQ/0rYoZ7dInVC70aCQeZT2V4YnxyMR18ucFp0O+zfFurlJf/TV8FctOdpt1txjSVVDrDgeOJiuaP5BsOB4nG4mQavJhTNJZABHZ6zPhbTG4esUgsQVsF17RJiauXl0DS4c4DAcB9XHjNTWNWsJ+emub8wtK5o7wYGQoymUzx9JmZGpSscpHYOZ4z0Ee/r/yhzLk2B3xs6OsmEksQjcW5erCf3m4vDRYXhPw+xhO2uNVq5CWQTOT0aSAirwVO169IxtRHpfmoCsl2RDcyE7CqEh1N1OR+RJyld/fFzjkd7RVcc9Dv4/T0HPOplblei6kfL4HkPwO/JyLHRCSGM2/jN+pbLGNqLxqLE+jp5LL1a6q+1pUDffR2tTe0wz129jxnZ+ar7mjPGgkHODoxQ3x2oaJrZhe4WinDos3yWbLuqqpPAS8QkT5AVHWq/sUypvYisTjD4fIy/hbT3ibOkNkGBpJIjTras3KDRyX9SIOBCwtcDa2tPlib5uGpEVREfgG4HvBl/xOq6p11LJdpQaqKKp6SANbazFyKwyeneOX1oZpdcyS8ls/+y1HGEufpaFv+FRl+/PQZujvauCa0VJIJb3a5waOns52rB8vLIgy5a7fXv58kk1FEqMmXAlO9JQOJiPwFsAZ4OfBXOEkUf1zncpkW9K4v7UMVPvXm5y77ax88niCj3lN+eLF7a5CFtPLCj36/Ztcs157L1tLZXpsgFujp5Dkb+1jf20VHBddczjQpb/38I1y2bg1/eOuOur+WWZqXGsmLVHWXiOxX1T8QkT8Bvl7vgpnWkskoPzg8Aeo8Xu5aSS072rNuvnYjn3jDMOcX0jW7ZrlecPm6ml7vU//huXRXOMcmuKaTro62us8lSS6k+eFTp5k8X7svBaY6XgJJ9q9iVkQ2A2eAy+tXJNOKnj4zw1TSycd59PQ0z9lYm+YYr6KjccLreljfVzBZdEU629t4/fOGana9laCaZjIRYVPAx1id06Rksx3b5MeVw8tXj78TkSDwceAnwDNcmtbdmJJyO6UjDchRFY3VZiKiKW3QX/+127N/S6em5kg3eEKocZQMJG5CxQdVNa6qXwMuA67NTXNijBeRWJzernb6ujuWfVGoU1NJjsfP17RZyxQW8td/dnu2mTKdUc5M2wTIlaBkIFHVDPAnOc/n8tK8G+NJNBZn11CQXUOBZV8UqhYZf4032TQp9VxyNxqLs6arHaDuzWjGGy9NW98VkdvFxtmZCiUX0jzuLgU7HA5yaGyS5DJ2UEdjcdrbhOs3W+dsvQ36fcynMsRnF5Y+uALZbMcvv2YjYAtprRReAsnv4iRpnBORSRGZEpHJOpfLtJBDOUvBjoSDpDLKYyeW708oEotzbaifHvdbrKmf7BDgetUUstmOX7XDmQ9kHe4rg5eldvtVtU1Vu1TV7z6vfDEHs+pEc4bejoSXN0dVJqNER+M1SyNiSgvlzG6vh2y245uuGaCjTWwhrRXCy4TElxXarqr/VPvimFYUiV28FGzI71u2HFVHTzvDjkdsxNayWFxyt06BJBI7x1Ub+/D7Ohlcho59442XeSTvzXnsw1kj/VHAFpgynuRnqB0JB5dtUajF2lAZy8aaym3s70akPkvuZrMd/9x1Tv/IoL/baiQrhJemrV/K+XkFsAM4Wf+imVYQn53n6dMzFycEDAd59sws52bmS5xZG9FRZ9jxlQPl544y5etsb2N9b3ddmrZGz12c7dgW0lo5KsmFMIoTTIxZUnT00qG32ceRZaiVRNxhx+0NSBS5WoUC3XXpbN/n1i6zE0uXY/Kj8cZLH8n/BLKDwtuAESBaz0KZ1hGNxS9ZCnbnUAARZ192GGc9JBfSHBqb5O0vuaJur2EuFfL3MHputubXjcbi+DovZDsO+X3MzKeZSi5UvUKkqY6XPpK9OY9TwJdU9V/rVB7TYiKx+CVLwfZ1d3DVxr66d7hfGHZs80eWUyjQzd5nz9b8upFYnB2bA4vZjnNHiFkgaSwvgeReIKmqaQARaReRNapa+68cpqWoKtFYnJuvvbTWMRIO8r3HT6KqdVtT4kLG37V1ub4pLOT3EZ9dILmQxtdZm7k7C+kMB48n+JUXXHbR64AzZ2W5k4Cai3npI3kQ6Ml53gP8Y32KY1rJ6LnznCmyFOxwOMi52QViZ8/X7fWjecOOzfIY9Nd+LsmT41PMpTIX9bUtDjW2fpKG8xJIfKo6nX3iPrZ1NM2SSq0Bku0w3VfHBI7RUcv42wjZD/hadrgX+luqR8AylfESSGZEZHFJOxF5HlC/r5GmZURj8aJLwV4T6sfX2Va3BI7ZYcc2f2T5herwAR+NxVnX28XQ2guNI77OdoJrOm0I8ArgpY/kd4CvisgJ9/km4JfrVyTTKiKxODu2BAouBdvZ3saOzYG6pZRfHHZsNZJlV48mp0gszkg4eEl/WsjvW5Y14k1pXiYkPgJcC/wm8F+A61T10XoXzDS3hXSGgycSJVO3j4SDHDwxyUI6U/PXjxxzhx3XcI12402/r5Pervaa1RSmkgscmZgu2EzpTEq0BpJGWzKQiMhvAb2qelBVDwB9IvJf6l8008wOn5wiuZApmSxxOBxkPpXhyfGpmr9+dPTSYcdm+QwGfDVr2jpwPIEqDBcYxm01kpXBSx/Jr6vq4oB/VT0H/LqXi4vILSLypIgcEZH3F9h/mYg8KCL7ReRhERlyt79cRCI5P0kRudXdd7mI/EhEfioiXxaRLm+3apbTYudoiaalbG1lX43nk2SHHVvG38ZxPuBrE0hKDdoY9Ps4MzNXl1qt8c5LIGnLXdRKRNqBJT+83eM+Bbwa2A7cISLb8w77BHC3qu4C7gQ+CqCqD6nqiKqO4CSHnAW+657zx8AnVfUq4Bzwdg/3YJZZtnM0vK6n6DFDa3tY39tV85Ty2WHHtiJi49QykERjcbatX0NwzaUfO6GAD1Vn/XbTOF4CyQPAV0TkZ0XkZuBLwD94OO8G4IiqHlXVeeAe4LV5x2zHmacC8FCB/QCvB76jqrNuQLsZZ5IkwBeAWz2UxSyzaCzB8FCg5GRDEWE4HKx5ICn1DdYsj1DAx6mpOTKZ6pfcjcYSRWuX2RFiNpeksbwEkvfhfNj/JvBb7uP3ljzDsQWI5TwfdbfligK3u49vA/pFZH3eMW/CCV4A64G4qqZKXBMAEXmHiOwVkb0TExMeimtqZXouxeFTU56aloaHghyZmGYqWbulWUsNOzbLIxTwkcoop2eqqymMJ5KMTyaLfikYtECyIngZtZVR1b9Q1der6u3A/cC7PVy70FfR/K8n7wFuFJF9wI3AcZx8Xs4FRDYBO3FqRV6vmS33Z1R1j6ruGRgY8FBcUysHRp3OUS81gpGtQVSdc2ql1LBjszwWJwtW2RGerV0W+1Kyqc4LaRlvPP1PE5ENIvKbIvJPwMPAoIfTRoFwzvMh4ETuAap6QlVfp6q7gQ+623I/Ud4IfENVs19XTwNBEcnOf7nkmqbxInnpvksZdofn1qrDPTvs2Ga0N9Zik1OVH/CRWJzOdmH7psKrewfXdNLV0Waz2xusaCARkX4R+TUR+Qfgx8BzgCtU9UpVfY+Haz8CXOWOsurCaaK6L+81NohItgwfAD6Xd407uNCshaoqTl/K691NbwG+5aEsZhllO0fX9i49oC64povLN/TWrJ8kO+zYZrQ31oVJidXN8YjG4ly3yV80+aOI1LRj31SmVI3kFM6IqLuAK1X13YDnJe3cfox34jRLHQK+oqqPicidIvIa97CbgCdF5DBOLeeu7Pkisg2nRvODvEu/D/hdETmC02fyWa9lMssjOlre0NvhoUDNlt71MuzY1N+Gvm7a26SqGkk6oxw4vnTtMmRrtzdcqRQpv4dTi/g08Lci8uVyL66q9+P0qeRu+1DO43u5MAIr/9xnKNCRrqpHcUaEmRXo5GSSsUSyrKal4XCQb0ZOMJY4z6ZA8eHCXngZdmzqr71N2NjfXdVkwacmppmeSy3Z1zYY8LF/GVbbNMUVrZGo6idV9WeA1+B0cn8T2Cwi7xORq5ergKa5LNYIymhayn5Q1KJ5y8uwY7M8Bv3VzW5fqqM9a1PAx1giidPybRrBy6ito6p6l6ruBJ4PBIDv1L1kpilFYnE62op3jhZy3SY/ne1SdYd7OcOOTf1V2+QUicXp93VwxYbekscN+n3MpzLEZ2s3hNyUp6zxkap6QFV/T1WvrFeBTHNbqnO0EF9nO9dt8lddI8kOO7ZAsjKEAtV1gkdjcYaHgrS1la5d1mqEmKmcDbQ3NZPJKPtHS2f8LWYkHOTAaIJ0FTOhraN9ZRn0+5ieSzE9l1r64DzJhTRPjE8VTNSYLxToBiyQNJIFElMz2c7RSmoEw0NBZubTHDk1vfTBRURjcS7zOOzY1N+mKtYlOXjc+VLhZdDGhcmPFkgaxQKJqZlqclxlO+erad6KjsYtv9YKUs1SuOX8LW3sr/3SvqY8RYf/isgBiqQfAXAz9hqzKDoap7976c7RQi5f30u/r4PIaJw3Pj+89Al5Khl2bOqrmpUSo6MJNgd8bHSDUSldHW1s6Ou22e0NVGoeyS+6v3/L/f1/3N9vxknrbsxFIrE4u8KBJTtHC2lrE4aHgkSOVVYj8TpU1CyfajrBI7FzZb2XoUC39ZE0UKl5JM+q6rPAi1X1v7kjtg6o6vuBVy1fEU0zSC6keWJsqqqmpZFwkCdPTnF+Pl32uVF32PH1m70POzb11dPVjt/XUXaN5Mz0HLGz58v6W7I0KY3lpY+kV0Rekn0iIi8Cym+7MC3tsRMJUh47R4sZDgdJZ5SDJ8rPBBypYNixqb9NgZ6yawrZdDnl1EiqnfxoquMlkLwd+JSIPCMiTwP/C/iP9S2WaTaRmPPhX02NJDvUs9wO9+ywYy9DRc3yqmTt9kgsQZvAzi3e38+Q38e52QWSC+XXZk31SvWRAKCqjwLDIuIHJC/NuzGA8+HvtXO0mI39PrYEexb7O7w6ejqbk2ltxa9t6iPk7+aJscmyzonG4lw92E9v95IfTxdeJ3BhhNhl663BZLktWSMRkUER+SzwZVVNiMh2EbF10s1FIrHyMv4WMxwOlB1I9h3LDhW1GslKE/L7mJieYyGd8XS8qjrZo8tsIq1mhJipnpemrc/jpILf7D4/DPxOvQpkms/ZmXmOnZ2tyRyOkXCQ0XPnOT3tPWvshWHHfVW/vqmtwYAPVZiY8vZ+PntmlvjsQtnryVialMbyEkg2qOpXgAwsrjNiDZFmUbSGQ2+z30TL6SepZtixqa9yP+DLWV0z12Cg8smPpnpeAsmMiKzHnZwoIi8ArJ/ELIrE4mV3jhazY0uANvEeSLLDjm0i4sq02HfhsckpEovT09nO1YPl1S77uztY09Vus9sbxEtv1u/iLJF7pYj8KzAAvKGupTJNJTpafudoMb3dHVw92E9k1Nt3lcdOTJLKqKVGWaHKrZFER+Ps3BKgo7287E3ZJXetRtIYXt6tx4AbgRcBvwFcDzxRz0KZ5qGqi+m+a2UkHCQai3taqKia/F6m/tb1dtHV3uYpkMynMjx2YrLiYdzVpq03lfPyFfKHqvpcnIACgIj8BHhu3UrVwr5zYIz7oicaXYyaWUhnODe7UNPUJMPhIPc8EuM/fWEvXR2lv+s8dmKSTVUOOzb1IyJs9Hfz7egYx86Uzqw0O59mPpWp+G8p5Pfxo6fPVnQuwFf3xljf18XN1w5WfI1C1xzo7+amazbW7JorUamkjSGcNdN7RGQ3znK7AH5gzTKUrSV99l+e5tDYJFvWts6a4sNDAW66ZqBm17vpmgGGw0Fi55ZO6ebrbOOXK0jyaJbPbbu38MBj4zw1sfQSAc/dGuTFV26o6HWykx8zGS174IWq8pG/P8S2Db01CySZjPKH336c4XBw9QYSnHxabwWGgD/N2T4F/F4dy9TSxhJJXnl9iE/+8kiji7JibQr08K3fenGji2Fq5N2vvIZ3v/Kaur9OyO8jlVFOz8wtppb36tkzsyTOL3DoxCRzqTTdHdWn2nnmzAyTydSq6LcpGkhU9QvAF0TkdlX92jKWqWVPWac8AAAfnklEQVRlMsqpqeTiSBZjTO1cWOCq/ECS7WubT2c4VGXy0fxrroZ+Gy8pUr4mIr+A08nuy9l+Zz0L1orOzs6zkNbFkSzGmNpZXJFxMslOyuuwj7jZo1MZZ/BILQJJdgj7ZDLF+fk0PV2tm1DUS4qUvwB+GXgXTj/JG4DL6lyulpT9ZjJogcSYmgsFKp/dHh2N89ytaxno765qlc5cuUPYW33GvZfhvy9S1V8DzqnqHwAvBKx3swLZtlJr2jKm9jb0ddPeJmWv3Z4ddjyyNchIOFh2rrdC5lJpDp2YXJyk2+rNW14CyXn396yIbAYWgMvrV6TWlZ11a01bxtRee5sw0Ndd9uz2J8YnnWHHQ04gOXp6hsTsQlVlOTQ2xXw6wy07QgCMT55f4ozm5iWQfFtEgsDHgZ8AzwD31LNQrerkZJI2gQ19XY0uijEtqZL1Ty7kigtcyPU2Wl2tJHLsHACvut4ZSjye8J6EtBl56Wz/Q/fh10Tk24DP1iSpzHgiycZ+X9npH4wx3oT83Tw1MVPWOfticTb0dbMl2IO/pxNwgsvLrq58blR0NMHG/m6uHOijr7uj5YcAl5qQ+LoS+1DVry91cRG5Bfj/gXbgr1T1j/L2XwZ8Did/11ngV1R11N23FfgrnP4YBX5eVZ8Rkc/jpGzJBrO3qmpkqbKsBOOTycUspcaY2tsU6OHfjpwp6xxnlFYAEcHv6+TKgd6qayTZkV8iwqC/u+X7SErVSH7J/b0RJ8/W993nLwceBkoGEhFpBz4FvAIYBR4RkftU9fGcwz4B3K2qXxCRm4GPAr/q7rsbuEtVvycifbhp7F3vVdV7l7q5lebkZJLLN9jqbcbUy6Dfx9Rcipm5lKckopPJBZ6amOG23VsWt42E1/KDw6dQVUTKX5ogMbvA0dMz3P68IaCydeubTdE2FlV9m6q+Dac2sF1Vb1fV23Hmk3hxA3BEVY+q6jxOv8pr847ZDjzoPn4ou19EtgMdqvo9tyzTqrp0vowVbiyRtI52Y+ooFOgGvA+33R9zGjZy83uNhAOcnp7neLyyDvJsbSY7F2VwFWQl9tJYv01Vx3KenwSu9nDeFiCW83zU3ZYrCtzuPr4N6HfXPrkaiIvI10Vkn4h83K3hZN0lIvtF5JMi0l3oxUXkHSKyV0T2TkxMeChufc3Op5hKpqxpy5g6ys7R8tqUlP3Q35WTvTobVKKxyrqCo7E4IrBzyBn6Gwp0c2pqjnRm6WzWzcpLIHlYRB4QkbeKyFuAv8epPSylUJ0w/1/yPcCNIrIPp9/jOJDCaXJ7qbv/+cAVOHm/AD4AXOtuXwe8r9CLq+pnVHWPqu4ZGKhdQsFKjdvQX2PqLlRmINl3LM4VA70E3E52gGtDfro62ojEzlVUhkgszpUDffh9nYtlSme0rOWjm82SgURV3wn8b2AYGAE+o6rv8nDtUS6euDgEXJQ/XVVPqOrrVHU38EF3W8I9d5/bLJYCvombtl5Vx9QxB/w1ThPaijdukxGNqbtyZrerKpFYnJG8tXS6Otq4frO/ohqJqhIdvTjFSrm1pGbkaRyqqn5dVf+r+/MNj9d+BLhKRC4XkS7gTTgrLS4SkQ0iki3DB3BGcGXPXSsi2arEzcDj7jmb3N8C3Aoc9Fiehlqc1W41EmPqZk1XB36ft+G2Y4kkp6fnCq5/MjwU5MDxBKl0psCZxR2Pn+f09PxF19wUcJaMaOUO96KBRET+xf09JSKTOT9TIjK51IXdmsQ7gQeAQ8BXVPUxEblTRF7jHnYT8KSIHAYGgbvcc9M4zVoPisgBnGayv3TP+aK77QCwAfhI2XfdANkJSVYjMaa+vK6UWGp1zd1bg5xfSHP45NJrqBS8Zk4tZ9AdANDKHe6l0si/xP3dX+nFVfV+4P68bR/KeXwvUHAYrztia1eB7TdXWp5GGk+cp9/XwZqu6tc1N8YUN+j3efr2H43F6Wpv49pNl37E5c5w377Z7/m1o7E4XR0XX3NDbzcdbVJ26pZmUqpGsq7Uz3IWshWMT9rQX2OWQ8jvvUZy3WZ/wUWsLlu/huCaTiLHypuYGInF2bHZT2dO9oq2NmFjf3fZySSbSamvx4/ijLIqNvrqirqUqEWNT85Zs5YxyyAU8HF6eo5UOlM0HVE6oxw4nuCNewonMhcRhoeCZc1wT6UzHDie4I4btl6ybzDgrZbUrEpNSLxcVa9wf+f/WBAp00mbjGjMsggFfGQUJkoMt/3pqSlm59MMh4svgDUcDnL45BQzcylPr3v45DTJhUzBPpdNqzWQ5BKRtSJyg4i8LPtT74K1klQ6w8S01UiMWQ5e5pJEFzva1xY9Znc4SEbhwHFvw4BLdd4P+n0t3bTlZYXE/wT8E87oqz9wf3+4vsVqLaen50ln1FZGNGYZeJm3EYnF8fs62LZ+TdFjdrkz072umBiNxVm7ppOt6y69ZsjvY2Y+zVSyunVOViovNZLfxplF/qyqvhzYDTQ+50gTGbc5JMYsGy+TEiOxBMNudt5i1vd1E17X47mfJDoaL3rNxTK1aK3ESyBJqmoSQES6VfUJ4Jr6Fqu1LKZHsaYtY+pu3ZouOtulaCCZnU9x+OQUuws0QeUbCa/1NHJrZs655vBQ4Wsu1pJatJ/ESyAZdVdI/CbwPRH5FnmpTkxptla7McunrU1K9kkcPD5JOqMFZ7TnGx4KcCKR5NQSAeDA8QQZLdw/AuXnAGs2XlZIvM19+GEReQgIAP9Q11K1mPHJJJ3twro1tsSuMcshVGJS4oWldb3USJxjIrE4r7w+VPS4yBLXzH6JbNXZ7aUmJP69iLxZRBZXYlLVH6jqfe76Isajk+4Su21t5S+SY4wp32CJNCmRWJyhtT1s6Cu4AsVFdmwJ0N4mS/aTRGNxtq5bw7rewl8WfZ3tBNd0rsqmrc8Avwg8IyJfFpFb3eSLpkxjiaQ1axmzjLI1EtVL1wCJxOKeaiPgBIBrQ/1LZgLOLq27ZJkSrZlKvtSExG+p6h3AVpxldd8CHBORz4nIK5argK3gpKVHMWZZhfw+kgsZJs9fPJlwYmqO4/Hznjras0bCQaKxOJkiC1OdmkxyIpFcMjg5OcAqW3VxpfOyHsl5Vf2y21fySpzhv9ZH4pGqOnm2rEZizLIZLDIEuJz+kazhcJCpuRRHT88U3H9hImLxWfKwSmskWSIyKCLvEpF/xRm59V3geXUvWYuYmksxO5+2Gokxy2hTsUAyGqe9TdixufSHfq6RxaV3C/eTREfjdLQJ1y9xzVDAx5mZORbKXOOkGZTqbP91Efk+8BOcNdT/m5t7632qGlm2Eja57BBEW6vdmOVzYbjtxU1JkVicawb76em6NONvMVcO9NHX3bFY88gXicW5dlM/vs7S1wwFfKjCqanWq5WUqpG8CPgjIKyq71LVf12mMrWUMVur3Zhlt9HvjMjKbUrKZJRoGR3tWe1tws4tgYIjtzIZZX8sUXQiYq5WnktSqrP9bar6XVVdrIeJyIeXpVQtxNKjGLP8ujvaWdfbdVHT1jNnZphMppbsyyhkOBzk0NgkyYX0RduPnp5mai615IgtaO212z1l/83xmqUPMbmyTVvZb0jGmOUx6PddNAEwW6MolfG3mJFwkIW08vjYxauMR9xhwV4CiZccYM2q3EBiM+rKND6ZZF1v15Ltp8aY2tqUNykxcixOb1c7z9nYV/a1inW4R2Nx+ro7uGJg6WuuXdNJV0dbS85uLzeQ2GitMp2cTFr6eGMaIL9GEhlNsHPImalerlDAR8jvu6TDPRKLs8vjNUXE8zLAzcbL8N+PiYhfRDpxkjaeFpFfWYaytQRnrXZr1jJmuYX8Ps7MzDOXSjOXSnPoxGTZHe25hsOBi2okyYU0h8bKu2apHGDNzEuN5JWqOomTLmUUZyjwe+taqhYybulRjGmIUMD5Andqco5DY1PMpzOMeBhdVcxwOMgzZ2aJzzqpBh8fmySVUU8jtrIGA75V27TV6f7+eeBLqnq2juVpKfOpDKen561py5gGyF0DZHFp3a2VB5LcTMDA4jolu8u4ZsjfzViicA6wZuYlkPydiDwB7AEeFJEBoPVCah2cmrKhv8Y0yqZAD+C0CkRicTb2d1f1f3HnlgAiLCZwjI7GCfl9ZX1RDAV6mE9liM+21pK7XnJtvR94IbBHVReAGeC19S5YK7AFrYxpnGzQOOnWSJZaWncp/b5OnjPQRyR2DnBqJl6G/RYqU6v1k3jpbH8DkFLVtIj8PvA3wOa6l6wFZGfVWiAxZvn5ezrwdbbx5PgUR0/PlP2hX8hIOEh0NMHZmXmePTNbdud9tt9m1QUS4P9V1SkReQnwKuALwKfrW6zWMObm+bGmLWOWX3a47YNPnAK8TRpcynA4yNmZef7+wJj7vLxZ8tlmsGLLADcrL4EkmxPgF4BPq+q3AFvgyoOTk0m6O9oI9HQufbAxpuYG/T7OzswjAjuHyk+Nki8bjO7+t2cQgV1ljgLb2O8EkrFVGEiOi8j/Bt4I3C8i3R7PQ0RuEZEnReSIiLy/wP7LRORBEdkvIg+LyFDOvq0i8l0ROSQij4vINnf75SLyIxH5qbty44oNauOTc4QCvqraZY0xlcs2K1850IffV/0XumtC/XR3tPHTU9NctdHJClyOro42NvR1t9wQYC8B4Y3AA8AtqhoH1uFhHomItAOfAl4NbAfuEJHteYd9ArhbVXcBdwIfzdl3N/BxVb0OuAE45W7/Y+CTqnoVcA54u4d7aIiTCVsZ0ZhGygaScuZ6lNLZ3saOLYGqrhkKdK++PhJVnQWeAl4lIu8ENqrqdz1c+wbgiKoeVdV54B4uHe21HXjQffxQdr8bcDpU9XtuGaZVdVacr/Y3A/e653wBuNVDWSry9OkZfnLsXMXn28qIxjRW9otcNfNH8mWbtyq9ZiumSfEyauu3gS8CG92fvxGRd3m49hYglvN81N2WKwrc7j6+DegXkfU4s+fjIvJ1EdknIh93azjrgbiqpkpcM1vud4jIXhHZOzEx4aG4l/rQtw7y+984WNG5i0vsWo3EmIa5cqAPEbhh27qaXfNnLl+HCDy/wmvm5wBrBV6att4O/IyqfkhVPwS8APh1D+cV6hjIn875HuBGEdkH3AgcB1JAB/BSd//zgSuAt3q8prNR9TOqukdV9wwMDHgo7qWGh4I8eXKK8/PppQ/Oc252gflUxma1G9NAL71qAw+/5yauCfXX7Jqv2D7Iw++5iasHK7tmyO/j3OzCJWubNDMvgUS4MHIL97GX3uNRIJzzfAg4kXuAqp5Q1dep6m7gg+62hHvuPrdZLIWzVvxzgdNAUEQ6il2zlkbCQdIZ5eCJRNnnZquu1rRlTOOICJet711R18x+JrRSrcRLIPlr4Eci8mF3hcR/Bz7r4bxHgKvcUVZdwJuA+3IPEJENIpItwweAz+Wcu9ZNxwJOv8jj6iSoeQh4vbv9LcC3PJSlIrvcMeL5axB4kf0jsRqJMSbX4gJXLdRP4qWz/U+BtwFncUZJvU1V/8zDeSngnTgjvg4BX1HVx0TkThHJrrR4E/CkiBwGBoG73HPTOM1aD4rIAZwa0F+657wP+F0ROYLTZ+IlqFVkY7+PLcEe9lUQSLKjMjZZjcQYk6MV06SUHATt1hb2q+oO4CflXlxV7wfuz9v2oZzH93JhBFb+ud8DdhXYfhRnRNiyGAkHK6qRjCeSiMBAv61FYoy5YHC1NW2pagaIisjWZSrPijMcDjB67jynp+fKOm88kWRDXzed7eUuQmmMaWX93R2s6WpfzMXXCrxMy9wEPCYiP8bJ/AuAqr6m+CmtIzvpKBqL87PXDXo+z4b+GmMKyeYAa6UaiZdA8gd1L8UKtnMoQJuUH0hOTiYZWrumjiUzxjSrUMC3mNS1FRQNJCLyHGBQVX+Qt/1lOPM9VoU1XR1cPdhfdof7+GSSPdvW1qlUxphmFvL7+NHTrbPYbKkG/D8Dpgpsn3X3rRq7tzod7l6Xx0wupInPLiyu0GaMMbmya7dnMq2x5G6pQLJNVffnb1TVvcC2upVoBRoeCjKZTPHMmVlPx9scEmNMKSG/j1RGOTMz3+ii1ESpQFLqU3BVfdXOroKWXWJzKdm1Bqyz3RhTyKC/tYYAlwokj4jIJTm1ROTtwKP1K9LKc/VgP2u62onGvKVKubBWu80hMcZcKjtRuVUWuCo1aut3gG+IyJu5EDj24KyOeFu9C7aStLcJO7YEiHjscM+mPrCmLWNMIYtpUlqkRlI0kKjqSeBFIvJyYIe7+e9V9fvLUrIVZiQc5PP/+gxzqTTdHe0ljx2fTNLX3UF/DVZkM8a0ng193bS3Scus3b7kPBJVfQgnUeKqNhIOMp/O8MTY1GKfSTEnJ5MM+q1ZyxhTWHubMNDXOislWv4Ojy50uC/dvDWWsJURjTGlZYcAtwILJB5tDvgY6O/2lMDxZCJp/SPGmJI2tdCSuxZIPBIRhoeCREZLB5JMRjk1NWdDf40xJYUCFkhWpZFwgKMTMyRmF4oec3pmjlRGrWnLGFPSoN/H1FyKmblUo4tSNQskZRgJO7mz9h8vXis56aaGthqJMaaU7DyzVuhwt0BShp1DSy+9Oz5pa7UbY5a2OLu9BZq3LJCUIdDTyRUDvSVHbo27qaGtRmKMKSWb1NVqJKvQSDhIJJYomgl4fDJJe5uwvs/mkRhjimultdstkJRpJBzk9PQcx+OFF6UZT8yxsd+ZtWqMMcX0dLXj93W0xMgtCyRlGglnl94tnMDx5KRNRjTGeNMqQ4AtkJTp2pCfrvY2okXmk9ha7cYYrwZbZO12CyRl6upoY/tmP5FjRQKJzWo3xni0KeCzPpLVaiQc5MDxBKl05qLt03MppudS1rRljPEk5PcxMTV3yWdJs7FAUoGRcJDzC2kOn5y+aPu4rYxojCnDYMBHRmFieq7RRamKBZIKZDMB5/eT2FrtxphyLA4BbvIOdwskFdi2fg2Bns5LZrhn/xg2WdOWMcaDVlm73QJJBUSE4XDwkhnulh7FGFOO7JdOq5GUICK3iMiTInJERN5fYP9lIvKgiOwXkYdFZChnX1pEIu7PfTnbPy8iT+fsG6nnPRQzEg5y+OTURZk7xxNJAj2d+DpLL8VrjDEA63q76GpvY3zS+kgKEpF24FPAq4HtwB0isj3vsE8Ad6vqLuBO4KM5+86r6oj785q8896bsy9Sr3soZSQcIKNw8PiFiYk2h8QYUw4RYaO/25q2SrgBOKKqR1V1HrgHeG3eMduBB93HDxXYv2IND1269O7JySSD1qxljClDyO9jLFE45VKzqGcg2QLEcp6PuttyRYHb3ce3Af0ist597hORvSLy7yJya955d7nNYZ8UkYZkR1zf1014Xc9FI7fGE0lCfkvWaIzxzlm73Zq2iimUtTA/Ze57gBtFZB9wI3AcyHY6bFXVPcB/AP5MRK50t38AuBZ4PrAOeF/BFxd5hxuI9k5MTFR3J0UMDwUXZ7gvpDNMTM8RclNDG2OMFyF37fZiGcWbQT0DySgQznk+BJzIPUBVT6jq61R1N/BBd1siu8/9fRR4GNjtPh9Txxzw1zhNaJdQ1c+o6h5V3TMwMFDTG8saCQc5kUhyajLJxNQcqjYZ0RhTnk0BH+cX0kwmm3fJ3XoGkkeAq0TkchHpAt4E3Jd7gIhsEJFsGT4AfM7dvjbbZCUiG4AXA4+7zze5vwW4FThYx3soaTET8GgiZ+ivNW0ZY7xrhbkkdQskqpoC3gk8ABwCvqKqj4nInSKSHYV1E/CkiBwGBoG73O3XAXtFJIrTCf9Hqvq4u++LInIAOABsAD5Sr3tYyvWbA7S3CZHYucXlMm1WuzGmHNl5Z2NNPJeko54XV9X7gfvztn0o5/G9wL0Fzvs3YGeRa95c42JWrKernWtD/URjCTa4KyJa05YxphyhFli73Wa2V2k4HCQ6Gmc8kaSrvY11vV2NLpIxpolsdEd6NnM6eQskVRoZCjKVTPFvT51hMNCN03VjjDHedHe0s763ywLJajay1elwP3A8Yc1axpiKDPp91rS1ml050Edvl5NbyzrajTGVCDX5SokWSKrU3ibsctOlWI3EGFOJQXdSYrOyQFID2YWuLH28MaYSIb+PMzPzzKXSjS5KReo6/He1GAkHAGvaMsZUJrsuyav/7J9pb6vtgJ3PvuX5bF2/pqbXzGeBpAZuvHojv/7Sy3nZ1fVJxWKMaW03XjPAbbu31KVG0tVR/4YnaeZEYV7t2bNH9+7d2+hiGGNMUxGRR93kuSVZH4kxxpiqWCAxxhhTFQskxhhjqmKBxBhjTFUskBhjjKmKBRJjjDFVsUBijDGmKhZIjDHGVGVVTEgUkQng2bzNG4DTDShOvbTa/UDr3ZPdz8rXavdU7f1cpqpLpuxYFYGkEBHZ62XGZrNotfuB1rsnu5+Vr9Xuabnux5q2jDHGVMUCiTHGmKqs5kDymUYXoMZa7X6g9e7J7mfla7V7Wpb7WbV9JMYYY2pjNddIjDHG1MCqCyQicouIPCkiR0Tk/Y0uTy2IyDMickBEIiLSdAuviMjnROSUiBzM2bZORL4nIj91f69tZBnLVeSePiwix933KSIiP9/IMpZDRMIi8pCIHBKRx0Tkt93tTfk+lbifZn6PfCLyYxGJuvf0B+72y0XkR+579GUR6ar5a6+mpi0RaQcOA68ARoFHgDtU9fGGFqxKIvIMsEdVm3L8u4i8DJgG7lbVHe62jwFnVfWP3IC/VlXf18hylqPIPX0YmFbVTzSybJUQkU3AJlX9iYj0A48CtwJvpQnfpxL380aa9z0SoFdVp0WkE/gX4LeB3wW+rqr3iMhfAFFV/XQtX3u11UhuAI6o6lFVnQfuAV7b4DKteqr6T8DZvM2vBb7gPv4Czn/yplHknpqWqo6p6k/cx1PAIWALTfo+lbifpqWOafdpp/ujwM3Ave72urxHqy2QbAFiOc9HafI/HpcC3xWRR0XkHY0uTI0MquoYOP/pgY0NLk+tvFNE9rtNX03RDJRPRLYBu4Ef0QLvU979QBO/RyLSLiIR4BTwPeApIK6qKfeQunzmrbZAIgW2tULb3otV9bnAq4HfcptVzMrzaeBKYAQYA/6kscUpn4j0AV8DfkdVJxtdnmoVuJ+mfo9UNa2qI8AQTgvMdYUOq/XrrrZAMgqEc54PAScaVJaaUdUT7u9TwDdw/oCa3Um3HTvbnn2qweWpmqqedP+jZ4C/pMneJ7fd/WvAF1X16+7mpn2fCt1Ps79HWaoaBx4GXgAERaTD3VWXz7zVFkgeAa5yRzF0AW8C7mtwmaoiIr1uZyEi0gu8EjhY+qymcB/wFvfxW4BvNbAsNZH9wHXdRhO9T25H7meBQ6r6pzm7mvJ9KnY/Tf4eDYhI0H3cA/wcTt/PQ8Dr3cPq8h6tqlFbAO5wvj8D2oHPqepdDS5SVUTkCpxaCEAH8LfNdk8i8iXgJpxMpSeB/w58E/gKsBU4BrxBVZum87rIPd2E02SiwDPAb2T7F1Y6EXkJ8M/AASDjbv49nH6FpnufStzPHTTve7QLpzO9HaeS8BVVvdP9jLgHWAfsA35FVedq+tqrLZAYY4yprdXWtGWMMabGLJAYY4ypigUSY4wxVbFAYowxpioWSIwxxlTFAolpCSLysIi8Km/b74jI/1rivOlS+2tQrgE38+o+EXlp3r6HRWSP+3ibm531VQWu8XE3m+vHKyzDTSLy7ZznHxGRB0Sk2y3D3px9e0Tk4ZzzVER+KWf/t0XkpkrKYVqXBRLTKr6EM8E015vc7Y30s8ATqrpbVf+50AEiMgQ8ALxbVR8ocMhvAM9V1fd6ecGcWcyF9n0QeDFwa85cgo0i8uoip4wCH/Tyumb1skBiWsW9wC+KSDcsJuLbDPyLiPSJyIMi8hNx1m25JONzgW/tfy4ib3UfP09EfuAmxXwgb/Zz9vjL3NfY7/7eKiIjwMeAnxdnbYueAuUOAd8Ffl9VL8myICL3Ab3Aj0Tklwu9jnvc50XkT0XkIeCPC/0Dici7gZ8HfklVz+fs+jjw+4XOAaJAQkReUWS/MRZITGtQ1TPAj4Fb3E1vAr6szozbJHCbm9jy5cCfuCkyluTmY/qfwOtV9XnA54BCmQP+HGftkV3AF4H/oaoR4ENuOUbyPryz7gb+XFW/WuS+XgOcd8//cqHXyTn8auDnVPXdBS71YuA/A6/OSTWe9UNgTkReXqgMwEcoHmiMsUBiWkpu81Zus5YA/5+I7Af+ESeN9qDHa14D7AC+56bn/n2cxHf5Xgj8rfv4/wAv8Xj9fwR+VUTWeDy+1Ot8VVXTRc47gvPv8Moi+4sGi2yTXH4fjzFZFkhMK/km8LMi8lygJ7twEfBmYAB4npti+yTgyzs3xcX/H7L7BXjMrRGMqOpOVS32YZzLa+6hj+Hkq/pqqb4Nj68zU+K4kzjNWp8sVPNQ1e/j3PMLipx/F9ZXYoqwQGJahttk8zBO81NuJ3sAOKWqC+6H6GUFTn8W2O6OZArgdJIDPAkMiMgLwWnqEpHrC5z/b1yoDb0ZZ5lTr/4rMAl81kOTW8Wvo6qHgdcBf+P23+S7C/hvRc79LrAWGPb6emb1sEBiWs2XcD7s7snZ9kVgjzvM9c3AE/knqWoMJ4vtfvf4fe72eZwU3H8sIlEgAryowOv+P8Db3OazX8VZK9sTtx/nLcAmnBpKKRW/jvtajwBvA+4TkSvz9t0PTJQ4/S4KN+uZVc6y/xpjjKmK1UiMMcZUxQKJMcaYqlggMcYYUxULJMYYY6pigcQYY0xVLJAYY4ypigUSY4wxVbFAYowxpir/F19pSpGIa+FNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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.9800000000000001\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.9533333333333334\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('data/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='neg_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.6913531708051797\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.6796748419090768\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='neg_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&list=PL5-da3qGB5IA6E6ZNXu7dp89_uv8yocmf) (14 minutes), [Cross-validation the right and wrong ways](https://www.youtube.com/watch?v=S06JpVoNaA0&list=PL5-da3qGB5IA6E6ZNXu7dp89_uv8yocmf) (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",
    "\n",
    "- Email: <kevin@dataschool.io>\n",
    "- Website: http://dataschool.io\n",
    "- Twitter: [@justmarkham](https://twitter.com/justmarkham)"
   ]
  },
  {
   "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 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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}