Cours/OPT1/03-model_tuning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Model tuning and hyperparameter optimization\n",
    "\n",
    "<i>Balázs Kégl (CNRS / PS-CDS) and Alex Gramfort (Inria)</i>\n",
    "\n",
    "## Introduction\n",
    "Virtually all learning models come with a set of parameters that should be tuned. To distinguish them from the parameters of the models set during training (learned parameters like neural weigths, tree cuts, or ensemble coefficients), we call them <i>hyperparameters</i>. Typical examples are the number of iterations and the tree depth in random forests, or the <code>C</code> or <code>gamma</code> of a support vector classifier. Scikit-learn has default values for all these parameters, but there is no guarantee that any of the classifiers perform optimally (or even close to optimally) using these default parameters. Worse: the optimal parameters vary depending on the data set, so tuning them is the task of the analyst, not the task of the provider of the library.\n",
    "\n",
    "Hyperparameters come in all shapes and types. Some of them are crucial to be tuned (like <code>n_estimators</code> in ensemble methods), the default value of some others works fine most of the time (say, <code>min_samples_split</code> of a random forest). Scikit-learn also has auxiliary parameters that have no effects on the model (such as <code>verbose</code>), so knowing which parameters to tune requires understanding a little bit what the particular model is about. Most of the hyperparameters have typical ranges, although they can depend of course on the characteristics of the data set. There are even so called <i>conditional</i> hyperparameters whose <i>existence</i> depend on the value of some other hyperparameters: for example, when you add a new layer to a neural network, suddenly you have 3-5 new hyperparameters to tune.\n",
    "\n",
    "Important hyperparameters can usually be interpreted as complexity parameters of the learned function. Sometimes the link is easy to see: an ensemble with a large number of deep trees have more degrees of freedom than a small ensemble of small trees. Sometimes the link between hyperparameters and complexity is more indirect, and understanding it requires some basic notion of <i>regularization</i>. One general notion that helps to understand the behavior of the learning algorithms with respect to their hyperparameters is <i>overfitting</i> and <i>underfitting</i> (this is why it is important to understand whether complexity increases or decreases as the hyperparameter grows). Let us call the model that we could obtain if we knew the distribution (or had infinite data) the <i>ideal</i> model. The goal of learning is to get as close to this model as possible, using only a finite training set. We may have two problems\n",
    "<ol>\n",
    "<li>\n",
    "Underfitting occurs when our model is not flexible enough to pick up the details of the ideal function, so even with infinite data we could not get close to the ideal function. The difference between the error of the ideal function and the error of the best function we could pick from our model class (parametrized by the hyperparameters) is the <i>approximation</i> error. Underfitting means that the approximation error is larger than optimal. Increasing the complexity of the model class decreases the approximation error (independently of the data set).\n",
    "<li>\n",
    "Overfitting occurs when our model is too flexible so it learns the training data \"by heart\" and does not generalize well to unseen test data. The difference between the error of the function (learned on the data set) and the error of the best function we could pick from our model class (parametrized by the hyperparameters) is the <i>estimation</i> error. Overfitting means that estimation error is larger than optimal. Increasing the complexity of the model class (while keeping the data set fixed) increases the estimation error.\n",
    "</ol>\n",
    "Note the crucial difference between the two notions: the level of underfitting is independent of the data set; it only depends on the data distribution and the function set (parametrized by the hyperparameters), whereas the level of overfitting, of course, depends on the data set. Here are some general rules that could guide you in this landscape.\n",
    "<ol>\n",
    "<li> Overfitting increases as the model complexity grows.\n",
    "<li> Overfitting increases as the data complexity (for example, the number of features) grows.\n",
    "<li> Overfitting decreases as the data size grows.\n",
    "</ol>\n",
    "\n",
    "If we could compute the data distribution or had access to an infinite data set we could determine the estimation and approximation errors and set the model complexity to minimize their sum. Of course, this is impossible so these notions are not constructive (but can help you to understand the underlying phenomena). In practice, we <i>estimate</i> the <i>generalization</i> error (the sum of the approximation and estimation errors) on a held out (validation) set. The pragmatics of hyperparameter optimization are thus relatively simple: find the hyperparameters that minimize the validation error. The practical difficulties are the following:\n",
    "<ol>\n",
    "<li> Our estimate has a variance: it is based on a finite validation set. Optimizing a \"noisy\" function can be tricky. Using cross validation schemes makes this even worse: since training and test sets overlap, complex correlation structures are introduced into the error estimates.\n",
    "<li> Evaluating the validation error at a vector of hyperparameter values requires training a model which may be slow. \n",
    "<li> The number of hyperparameter combinations increases exponentially with the number of hyperparameters.\n",
    "</ol> \n",
    "\n",
    "The first step of hyperparameter tuning is usually trial and error. Try some combinations and look at how the error behaves. Once you have a \"feeling\" about what range should be explored and the sensitivity of the error to the hyperparameters, you can do a grid search (using either <code>GridSearchCV</code> from scikit learn, or by coding it). In this notebook we coded a simple grid search so you understand how it works. Once the number of hyperparameters is large (say, larger than 2) grid search will become increasingly time consuming. The main inefficency is a result of it not being intelligent: it will exhaustively look at <i>all</i> combinations, even those that are \"trivially\" bad. \n",
    "\n",
    "The main principle to improve on grid search is smoothness in the hyperparameter space: we assume that models with similar hyperparameters are not very different. This allows sophisticated automatic methods to direct the search toward promising regions while also keeping exploring unexplored regions. Most of these methods are based on <a href=\"http://en.wikipedia.org/wiki/Bayesian_optimization\">Bayesian optimization</a>. We will work with the <a href=\"https://github.com/hyperopt/hyperopt/wiki/FMin\">hyperopt</a> package. It can be installed with:\n",
    "\n",
    "<code>pip install hyperopt</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# !pip install hyperopt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import some libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fetch the data and load it in pandas, and convert it into numpy\n",
    "We can use any data set here. By the end of this section, two numpy arrays, <code>X</code> and <code>y</code> shoud be created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"train.csv\", index_col=0)\n",
    "data['real_period'] = data['period'] / data['div_period']\n",
    "data.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = data[['magnitude_b', 'magnitude_r', 'real_period']].values\n",
    "y = data['type'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape, y.shape, X.dtype, y.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're good !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train a classifier\n",
    "\n",
    "without worrying too much about hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.5, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import cross_val_score, KFold\n",
    "\n",
    "clf = RandomForestClassifier(n_estimators=100, max_depth=8)\n",
    "cv = KFold(n_splits=5, random_state=42)\n",
    "scores = cross_val_score(clf, X_train, y_train, cv=cv, scoring='accuracy')\n",
    "print(\"Accuracy: {:.4f} +/-{:.4f}\".format(np.mean(scores), np.std(scores)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train a classifier using a grid of hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_estimators_list = [2, 3, 5, 10, 20, 30, 50, 100, 200, 500]\n",
    "max_depth_list = [1, 2, 3, 5, 7, 10, 15, 20, 30]\n",
    "accuracies = {}\n",
    "accuracy_stds = {}\n",
    "\n",
    "for n_estimators in n_estimators_list:\n",
    "    print(n_estimators)\n",
    "    for max_depth in max_depth_list:\n",
    "        clf = RandomForestClassifier(n_estimators=n_estimators, max_depth=max_depth, n_jobs=1)\n",
    "        scores = cross_val_score(clf, X_train, y_train, cv=cv, scoring='accuracy', n_jobs=3)\n",
    "        accuracies[(n_estimators, max_depth)] = np.mean(scores)\n",
    "        accuracy_stds[(n_estimators, max_depth)] = np.std(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(list(accuracies.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's do it using scikit-learn <a href=\"http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html\">GridSearchCV</a> object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "tuned_parameters = {'n_estimators': n_estimators_list, 'max_depth': max_depth_list}\n",
    "clf = RandomForestClassifier(random_state=42)\n",
    "\n",
    "gs = GridSearchCV(clf, tuned_parameters, cv=cv, scoring='accuracy', n_jobs=3)\n",
    "gs.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Best parameters set found on development set:\")\n",
    "print(gs.best_params_)\n",
    "print()\n",
    "print(\"Grid scores on training set:\")\n",
    "print()\n",
    "means = gs.cv_results_['mean_test_score']\n",
    "stds = gs.cv_results_['std_test_score']\n",
    "for mean, std, params in zip(means, stds, gs.cv_results_['params']):\n",
    "    print(\"%0.3f (+/-%0.03f) for %r\" % (mean, std * 2, params))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_results = pd.DataFrame(gs.cv_results_)\n",
    "cv_results.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_results[['mean_test_score', 'param_max_depth', 'param_n_estimators']].head(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_results[['mean_test_score', 'param_max_depth', 'param_n_estimators']].sort_values(by='mean_test_score', ascending=False).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot some learning curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "this_cv_results = cv_results[cv_results['param_n_estimators'] == 10]\n",
    "plt.errorbar(this_cv_results['param_max_depth'],\n",
    "             this_cv_results['mean_test_score'],\n",
    "             yerr=this_cv_results['std_test_score'],\n",
    "             fmt='o');\n",
    "plt.xlabel('max depth')\n",
    "plt.ylabel('accuracy');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "this_cv_results = cv_results[cv_results['param_max_depth'] == 15]\n",
    "plt.errorbar(this_cv_results['param_n_estimators'],\n",
    "             this_cv_results['mean_test_score'],\n",
    "             yerr=this_cv_results['std_test_score'],\n",
    "             fmt='o');\n",
    "plt.xlabel('nb estimators')\n",
    "plt.ylabel('accuracy');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now let's automate the search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All wrapper-type hyperparameter optimization will need a function with a parameter that represents the hyperparameter, and which returns a score (usually to be minimized). The following function thus trains a random forest classifier with two hyperparameters, stored in a pair <code>x_int</code>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "i = 0\n",
    "\n",
    "def objective(params):\n",
    "    global i\n",
    "    i += 1\n",
    "    print(params)\n",
    "    clf = RandomForestClassifier(n_estimators=300,\n",
    "                                 max_depth=int(params['max_depth']),\n",
    "                                 max_features=params['max_features'])\n",
    "    scores = cross_val_score(clf, X_train, y_train, cv=cv, scoring='accuracy', n_jobs=3)\n",
    "    score = np.mean(scores)\n",
    "    print(\"SCORE: %s\" % score)\n",
    "    df_result_hyperopt.loc[i, ['score'] + list(params.keys())] = \\\n",
    "        [score] + list(params.values())\n",
    "    return {'loss': 1. - score, 'status': STATUS_OK}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hyperopt import fmin, tpe, hp, Trials, STATUS_OK\n",
    "\n",
    "space = {\n",
    "    'max_depth': hp.quniform('max_depth', 2, 30, 1),\n",
    "    'max_features': hp.quniform('max_features', 0.7, 1, 0.05)\n",
    "}\n",
    "\n",
    "df_result_hyperopt = pd.DataFrame(\n",
    "    columns=['score'] + list(space.keys()))\n",
    "\n",
    "trials = Trials()\n",
    "best = fmin(fn=objective, space=space, algo=tpe.suggest,\n",
    "            max_evals=15, trials=trials)\n",
    "\n",
    "print(\"Best: %s\" % best)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_result_hyperopt.sort_values(by='score', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "<ol>\n",
    "<li> Take any scikit-learn classifier and tune it using hyperopt.\n",
    "<li> Download any other training set from, e.g., the <a href=http://archive.ics.uci.edu/ml/>UCI Irvine data repository</a>, load it, and optimize any scikit-learn classifier and tune it using pysmac or hyperopt.\n",
    "<li> Take the <a href=\"https://drive.google.com/file/d/0B_sb8NJ9KsLUd2pzdjZ3VGxZLWM/view?usp=sharing\">notebook</a> that describes the RAMP on variable star classification, assemble a pipeline of a parametrized feature extractor and a classifier, and optimize the full pipeline.\n",
    "</ol>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}