An introduction to Machine Learning with Scikit-Learn.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Global imports and settings\n",
    "from __future__ import print_function\n",
    "from __future__ import division\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 8)\n",
    "\n",
    "import numpy as np\n",
    "np.set_printoptions(precision=3)\n",
    "\n",
    "from IPython.html.services.config import ConfigManager\n",
    "from IPython.utils.path import locate_profile\n",
    "cm = ConfigManager(profile_dir=locate_profile(get_ipython().profile))\n",
    "_ = cm.update('livereveal', {'width': 1440, 'height': 720 })"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<center>\n",
    "    <img src=\"img/scikit-learn-logo.png\" width=\"40%\" />\n",
    "    <br />\n",
    "    <h1>An introduction to Machine Learning with Scikit-Learn</h1>\n",
    "    <br /><br />\n",
    "    Gilles Louppe (<a href=\"https://twitter.com/glouppe\">@glouppe</a>), April 23, 2015\n",
    "    <br /><br />\n",
    "    CERN\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# About myself \n",
    "\n",
    "* Researcher in machine learning, expertise in tree-based methods\n",
    "* Core developer of Scikit-Learn\n",
    "* Research Fellow at GS-SIS, working on graph and text mining problems\n",
    "* Open to collaborations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Prerequisites "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "not a chance (<ipython-input-3-8ea16afe3b06>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-3-8ea16afe3b06>\"\u001b[1;36m, line \u001b[1;32m2\u001b[0m\n\u001b[1;33m    from __future__ import braces\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m not a chance\n"
     ]
    }
   ],
   "source": [
    "# This is an IPython notebook, with executable Python code inside\n",
    "from __future__ import braces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "- Materials at  <a href=\"https://github.com/glouppe/tutorial-sklearn-lhcb\">https://github.com/glouppe/tutorial-sklearn-lhcb</a>\n",
    "    \n",
    "- Python distribution with scientific packages (NumPy, SciPy, Scikit-Learn, Pandas)\n",
    "\n",
    "- ... or Anaconda http://continuum.io/downloads\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Outline\n",
    "\n",
    "* What is machine learning?\n",
    "* Scikit-Learn\n",
    "* Classifying events with decision trees\n",
    "* Model evaluation and selection\n",
    "* Ensemble methods\n",
    "* Going further\n",
    "* Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# What is machine learning?\n",
    "\n",
    "Predict, understand or identify patterns in data from observations.\n",
    "\n",
    "Applications include: _Natural language processing, Computer vision, IR and advertisement, Robotics, Bioinformatics, __High Energy Physics__, ..._\n",
    "\n",
    "<center><img src=\"img/motivation.png\" width=\"70%\" /></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Scikit-Learn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Overview\n",
    "\n",
    "* Machine learning library written in __Python__\n",
    "* __Simple and efficient__, for both experts and non-experts\n",
    "* Classical, __well-established machine learning algorithms__\n",
    "* Shipped with <a href=\"http://scikit-learn.org/dev/documentation.html\">documentation</a> and <a href=\"http://scikit-learn.org/dev/auto_examples/index.html\">examples</a>\n",
    "* __BSD 3 license__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Community driven development\n",
    "\n",
    "- 10-15~ core developers (mostly researchers)\n",
    "- 250+ occasional contributors\n",
    "- __All working together__ on [GitHub](https://github.com/scikit-learn/scikit-learn)\n",
    "- Emphasis on __keeping the project maintainable__\n",
    "    - Style consistency\n",
    "    - Unit-test coverage\n",
    "    - Documentation and examples\n",
    "    - Code review\n",
    "- Join us!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Python stack for data analysis\n",
    "\n",
    "- The __open source__ Python ecosystem provides __a standalone, versatile and powerful scientific working environment__, including: [NumPy](http://numpy.org), [SciPy](http://scipy.org), [IPython](http://ipython.org), [Matplotlib](http://matplotlib.org), [Pandas](http://pandas.pydata.org/), _and many others..._\n",
    "\n",
    "<center> \n",
    "<img src=\"img/scikit-learn-logo.png\" style=\"max-width: 120px; display: inline\" />\n",
    "<img src=\"img/numpy-logo.png\" style=\"max-width: 120px; display: inline\" />\n",
    "<img src=\"img/scipy-logo.png\" style=\"max-width: 120px; display: inline\" />\n",
    "<img src=\"img/ipython-logo.jpg\" style=\"max-width: 120px; display: inline\" />\n",
    "<img src=\"img/matplotlib-logo.png\" style=\"max-width: 120px; display: inline\"/>\n",
    "<img src=\"img/pandas-logo.png\" style=\"max-width: 120px; display: inline\" />\n",
    "</center>\n",
    "\n",
    "- Scikit-Learn builds upon NumPy and SciPy and __complements__ this scientific environment with machine learning algorithms;\n",
    "- By design, Scikit-Learn is __non-intrusive__, easy to use and easy to combine with other libraries;\n",
    "- Core algorithms are implemented in low-level languages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Algorithms\n",
    "\n",
    "__Supervised learning:__\n",
    "\n",
    "* Linear models (Ridge, Lasso, Elastic Net, ...)\n",
    "* Ensemble methods (Random Forests, Bagging, GBRT, ...)\n",
    "* Support Vector Machines\n",
    "* Nearest neighbors\n",
    "* Neural networks (soon!)\n",
    "\n",
    "<center><a href=\"http://scikit-learn.org/dev/auto_examples/classification/plot_classifier_comparison.html\"><img src=\"img/classifiers.png\" width=\"90%\" /></a>\n",
    "<em>A comparison of (some of the) classifiers in Scikit-Learn</em></center><br />\n",
    "\n",
    "__Unsupervised learning:__\n",
    "\n",
    "* Clustering (KMeans, Ward, ...)\n",
    "* Matrix decomposition (PCA, ICA, ...)\n",
    "* Density estimation\n",
    "* Outlier detection\n",
    "\n",
    "__Model selection and evaluation:__\n",
    "\n",
    "* Cross-validation\n",
    "* Grid-search\n",
    "* Lots of metrics\n",
    "\n",
    "_... and many more!_ (See our [Reference](http://scikit-learn.org/stable/modules/classes.html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Classifying events with decision trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Input data\n",
    "\n",
    "Data comes as a finite learning set ${\\cal L} = (X, y)$ where\n",
    "* Input samples are given as an array $X$ of shape `n_samples` $\\times$ `n_features` <br />\n",
    "  _E.g. samples are events and features are their kinematic properties._\n",
    "\n",
    "* Output values are given as an array $y$ <br />\n",
    "  _E.g. whether the event is background or signal._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "X = np.random.rand(4, 3)\n",
    "y = np.random.randint(0, 2, len(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# E.g. for n_samples = 4, n_features = 3\n",
    "print(\"X =\")\n",
    "print(X)\n",
    "print(\"y =\", y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Loading data with `pandas`\n",
    "\n",
    "`pandas` is a library providing easy-to-use data structures and data analysis tools. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Download simulated events from https://archive.ics.uci.edu/ml/datasets/SUSY\n",
    "\n",
    "# Load CSV file with pandas\n",
    "import pandas as pd \n",
    "df = pd.read_csv(\"SUSY.csv\", header=None)     \n",
    "df.columns = [# 0 = background, 1 = signal\n",
    "              \"target\", \n",
    "              # Kinematic properties\n",
    "              \"lepton 1 pT\", \"lepton 1 eta\", \"lepton 1 phi\", \n",
    "              \"lepton 2 pT\", \"lepton 2 eta\", \"lepton 2 phi\", \n",
    "              \"missing energy magnitude\", \"missing energy phi\", \n",
    "              # Derived features\n",
    "              \"MET_rel\",\"axial MET\", \"M_R\", \"M_TR_2\", \"R\", \"MT2\", \n",
    "              \"S_R\", \"M_Delta_R\", \"dPhi_r_b\", \"cos(theta_r1)\"]\n",
    "df.target = df.target.astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Exploration and visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.head()  # df.describe(), df.info(), ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "df_sample = df[:10000]     # Draw a sample to make things faster\n",
    "\n",
    "# Visualize samples through a scatter matrix\n",
    "colors = [\"blue\", \"red\"]   # blue = background, red = signal\n",
    "features = [\"lepton 1 pT\", \"lepton 1 eta\", \"missing energy magnitude\", \"R\"]\n",
    "_ = pd.scatter_matrix(df_sample[features], marker=\"x\", \n",
    "                      c=df_sample.target.apply(lambda x: colors[x]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Need more? See `pandas` <a href=\"http://pandas.pydata.org/pandas-docs/stable/visualization.html\">visualization</a> documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Data in Scikit-Learn\n",
    "\n",
    "- Input data = Numpy arrays or sparse matrices ;\n",
    "- Leverage efficient vector operations ;\n",
    "- Keep code short and readable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Get data as Numpy arrays from the pandas data frame\n",
    "X = df_sample.drop(\"target\", axis=1).values\n",
    "y = df_sample.target.values\n",
    "\n",
    "print(X[:5])\n",
    "print(X.shape)\n",
    "print(X.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Supervised learning \n",
    "\n",
    "The goal of supervised learning is to build an estimator $\\varphi_{\\cal L}: {\\cal X} \\mapsto {\\cal Y}$ minimizing\n",
    "\n",
    "$$\n",
    "Err(\\varphi_{\\cal L}) = \\mathbb{E}_{X,Y}\\{ L(Y, \\varphi_{\\cal L}(X)) \\}.\n",
    "$$\n",
    "\n",
    "where $L$ is a loss function (e.g., the zero-one loss for classification)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Decision trees\n",
    "\n",
    "A decision tree is a cut-based partition of the input space, where each region (= leaf) is fit with a (simple) predictive model.\n",
    "\n",
    "<center>\n",
    "    <img src=\"img/tree-partition.png\" width=\"39%\" style=\"display:inline\" />\n",
    "    <img src=\"img/tree-simple.png\" width=\"60%\" style=\"display:inline\" />\n",
    "</center>\n",
    "<small>\n",
    "<pre>\n",
    "def build(L):\n",
    "    create node t\n",
    "    if the stopping criterion is True:\n",
    "        assign a predictive model to t\n",
    "    else:\n",
    "        Find the best binary split L = L_left + L_right\n",
    "        t.left = build(L_left)\n",
    "        t.right = build(L_right)\n",
    "    return t     \n",
    "</pre>\n",
    "</small>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## A simple and unified API\n",
    "\n",
    "All learning algorithms in scikit-learn, including decision trees, share a uniform and limited API consisting of complementary interfaces:\n",
    "\n",
    "- an `estimator` interface for building and fitting models;\n",
    "- a `predictor` interface for making predictions;\n",
    "- (a `transformer` interface for converting data.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Estimators\n",
    "\n",
    "An estimator is any object that learns from data; it may be a classification, regression or clustering algorithm or a transformer that extracts/filters useful features from raw data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Estimator(object):\n",
    "    def fit(self, X, y=None):\n",
    "        \"\"\"Fits estimator to data.\"\"\"\n",
    "        # set state of ``self``\n",
    "        return self\n",
    "            \n",
    "    def predict(self, X):\n",
    "        \"\"\"Predict response of ``X``.\"\"\"\n",
    "        # compute predictions ``pred``\n",
    "        return pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Import the decision tree class\n",
    "from sklearn.tree import DecisionTreeClassifier  # Change this to try \n",
    "                                                 # something else\n",
    "\n",
    "# Set hyper-parameters, for controlling the learning algorithm\n",
    "clf = DecisionTreeClassifier(criterion=\"entropy\")\n",
    "\n",
    "# Learn a model from training data\n",
    "clf.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Estimator state is stored in instance attributes\n",
    "clf.tree_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Predictors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Make predictions \n",
    "print(clf.predict(X[:5])) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compute class probabilities\n",
    "print(clf.predict_proba(X[:5]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.tree import export_graphviz\n",
    "export_graphviz(clf, out_file=\"tree.dot\", \n",
    "                feature_names=df.columns[1:], max_depth=2)\n",
    "!dot -Tpng -o tree.png tree.dot\n",
    "from IPython.display import Image\n",
    "Image(\"tree.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Pros and cons of decision trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "- _Non-parametric_ model, proved to be consistent.\n",
    "- Support _heterogeneous data_ (continuous, ordered or categorical variables).\n",
    "- Flexibility in loss functions (but choice is limited).\n",
    "- Fast to train, fast to predict.\n",
    "- Easily _interpretable_.\n",
    "- Low bias, but usually high variance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Model evaluation and selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)\n",
    "clf = DecisionTreeClassifier(criterion=\"entropy\", \n",
    "                             random_state=1).fit(X_train, y_train)\n",
    "print(\"Training accuracy =\", clf.score(X_train, y_train))   # Biased estimate\n",
    "print(\"Test accuracy =\", clf.score(X_test, y_test))         # Unbiased estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Beware of bias when you estimate model performance:\n",
    "- Training score is often an optimistic estimate of the true performance;\n",
    "- The same data should not be used both for training and evaluation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "When ${\\cal L}$ is small, prefer K-Fold cross-validation instead of train-test split for more accurate estimates.\n",
    "\n",
    "<center><img src=\"img/kfold.jpg\" width=\"70%\" /></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import KFold\n",
    "\n",
    "for train, test in KFold(n=len(X), n_folds=5, random_state=42):\n",
    "    print(train)\n",
    "    print(test)\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "scores = []\n",
    "\n",
    "for train, test in KFold(n=len(X), n_folds=5, random_state=42):\n",
    "    X_train, y_train = X[train], y[train]\n",
    "    X_test, y_test = X[test], y[test]\n",
    "    clf = DecisionTreeClassifier(criterion=\"entropy\", \n",
    "                                 random_state=1).fit(X_train, y_train)\n",
    "    scores.append(clf.score(X_test, y_test))\n",
    "\n",
    "print(\"%f +-%f\" % (np.mean(scores), np.std(scores)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Shortcut\n",
    "from sklearn.cross_validation import cross_val_score\n",
    "scores = cross_val_score(DecisionTreeClassifier(criterion=\"entropy\", \n",
    "                                                random_state=1), \n",
    "                         X, y, cv=KFold(n=len(X), n_folds=5, random_state=42), \n",
    "                         scoring=\"accuracy\")\n",
    "print(\"%f +-%f\" % (np.mean(scores), np.std(scores)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Under- and over-fitting\n",
    "\n",
    "- Under-fitting: the model is too simple and does not capture the true relation between X and Y.\n",
    "- Over-fitting: the model is too specific to the training set and does not generalize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.learning_curve import validation_curve\n",
    "\n",
    "param_range = range(1, 16)\n",
    "train_scores, test_scores = validation_curve(\n",
    "    DecisionTreeClassifier(), X, y, \n",
    "    param_name=\"max_depth\", \n",
    "    param_range=param_range, cv=5, n_jobs=-1)\n",
    "\n",
    "train_scores_mean = np.mean(train_scores, axis=1)\n",
    "train_scores_std = np.std(train_scores, axis=1)\n",
    "test_scores_mean = np.mean(test_scores, axis=1)\n",
    "test_scores_std = np.std(test_scores, axis=1)\n",
    "\n",
    "plt.xlabel(\"max_depth\")\n",
    "plt.ylabel(\"score\")\n",
    "plt.xlim(min(param_range), max(param_range))\n",
    "plt.plot(param_range, train_scores_mean, color=\"red\", label=\"training score\")\n",
    "plt.fill_between(param_range, train_scores_mean - train_scores_std,\n",
    "                 train_scores_mean + train_scores_std, alpha=0.2, color=\"red\")\n",
    "plt.plot(param_range, test_scores_mean, color=\"blue\", label=\"test score\")\n",
    "plt.fill_between(param_range, test_scores_mean - test_scores_std,\n",
    "                 test_scores_mean + test_scores_std, alpha=0.2, color=\"blue\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Best trade-off\n",
    "print(\"max_depth = %d, accuracy = %f\" % (param_range[np.argmax(test_scores_mean)],\n",
    "                                         np.max(test_scores_mean)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Hyper-parameter search\n",
    "\n",
    "- Learning algorithms are not black boxes. \n",
    "- Finding good hyper-parameters is crucial to control under- and over-fitting, hence achieving better performance.\n",
    "- This is automatized with the `GridSearchCV` estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.grid_search import GridSearchCV\n",
    "grid = GridSearchCV(DecisionTreeClassifier(),\n",
    "                    param_grid={\"max_depth\": range(1, 16),\n",
    "                                \"criterion\": [\"gini\", \"entropy\"],\n",
    "                                \"min_samples_leaf\": [1, 5, 10, 50]},\n",
    "                    scoring=\"accuracy\",\n",
    "                    cv=5, n_jobs=-1)\n",
    "grid.fit(X, y)\n",
    "\n",
    "# Warning: don't report these numbers in your experiments!\n",
    "print(\"Best score = %f, Best parameters = %s\" % (grid.best_score_, \n",
    "                                                 grid.best_params_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Selection and evaluation, _simultaneously_\n",
    "\n",
    "- The resulting `grid.best_estimator_` model is not independent from `grid.best_score_` since its construction was guided by the maximization of this quantity. \n",
    "\n",
    "- As a result, the optimized `grid.best_score_` accuracy _may_ in fact be a biased, optimistic, estimate of the true performance of the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import cross_val_score\n",
    "from sklearn.grid_search import GridSearchCV\n",
    "\n",
    "scores = cross_val_score(\n",
    "            GridSearchCV(DecisionTreeClassifier(),\n",
    "                         param_grid={\"max_depth\": range(1, 16),\n",
    "                                     \"criterion\": [\"gini\", \"entropy\"],\n",
    "                                     \"min_samples_leaf\": [1, 5, 10, 50]},\n",
    "                         scoring=\"accuracy\",\n",
    "                         cv=5, n_jobs=-1), \n",
    "            X, y, cv=5, scoring=\"accuracy\")\n",
    "\n",
    "# Unbiased estimate of the accuracy\n",
    "print(\"%f +-%f\" % (np.mean(scores), np.std(scores)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Ensemble methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Random Forests\n",
    "\n",
    "- Single decision trees suffer from high variance. \n",
    "- This can be fixed by combining several randomized trees into a single model. \n",
    "\n",
    "<center>\n",
    "    <img src=\"img/forest.png\" />\n",
    "</center>\n",
    "\n",
    "- Randomization strategies:\n",
    "    * Bootstrap samples\n",
    "    * Random selection of $K$ split features\n",
    "    * Random selection of cut-off threshold\n",
    "- From the bias-variance decomposition, it can be shown that the resulting generalization performance is better.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)\n",
    "\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "rf = RandomForestClassifier(n_estimators=1000).fit(X_train, y_train)\n",
    "print(\"Training accuracy =\", rf.score(X_train, y_train)) \n",
    "print(\"Test accuracy =\", rf.score(X_test, y_test))   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Variable importances\n",
    "\n",
    "- Scores derived from a tree-based model to assess the individual importance of input features.\n",
    "- Can be accessed through  `estimator.feature_importances_`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Feature importances from totally randomized trees (TRT)\n",
    "from sklearn.ensemble import ExtraTreesClassifier\n",
    "trt = ExtraTreesClassifier(n_estimators=1000, \n",
    "                           max_features=1, max_depth=3).fit(X_train, y_train)      \n",
    "\n",
    "# Plot importances\n",
    "feature_names = df.columns[1:]\n",
    "\n",
    "importances = pd.DataFrame()\n",
    "importances[\"DT\"] = pd.Series(grid.best_estimator_.feature_importances_,\n",
    "                              index=feature_names)\n",
    "importances[\"RF\"] = pd.Series(rf.feature_importances_, index=feature_names)\n",
    "importances[\"TRT\"] = pd.Series(trt.feature_importances_, index=feature_names)\n",
    "importances.plot(kind=\"barh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Disclaimer: Importances are measured only through the eyes of the model. _They may not tell the entire nor the same story!_\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Partial dependence plots\n",
    "\n",
    "Relation between the response Y and a subset of features, marginalized over all other features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn.ensemble.partial_dependence import plot_partial_dependence\n",
    "\n",
    "gbrt = GradientBoostingClassifier(n_estimators=100, \n",
    "                                  learning_rate=0.1, \n",
    "                                  max_leaf_nodes=5)\n",
    "gbrt.fit(X, y)\n",
    "\n",
    "plot_partial_dependence(gbrt, X, feature_names=feature_names, features=[0, 6])\n",
    "plot_partial_dependence(gbrt, X, feature_names=feature_names, features=[(1, 4)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More detailed information can be found <a href=\"http://scikit-learn.org/dev/modules/ensemble.html#partial-dependence\">here</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Sample and class weights\n",
    "\n",
    "- Most learning algorithms assume that input samples are _i.i.d._\n",
    "- When it is not the case,\n",
    "    * sample probabilities can be specified through `sample_weight`\n",
    "    * class probabilities can be specified through `class_weight`\n",
    "- These are supported by _some_ estimators (but not all)\n",
    "- Beware that _you are changing the objective function_!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# Make classes equiprobable (shortcut: class_weight=\"auto\"), \n",
    "# even if classes are unbalanced\n",
    "clf = RandomForestClassifier(class_weight={0: 0.5, 1: 0.5}) \n",
    "\n",
    "# Specify sample weights\n",
    "sample_weight = np.random.rand(X_train.shape[0]) \n",
    "clf.fit(X_train, y_train, sample_weight=sample_weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Going further"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Model persistence\n",
    "\n",
    "As any Python objects, estimators can be saved to disk for future reuse using `pickle`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "clf = DecisionTreeClassifier().fit(X_train, y_train)\n",
    "\n",
    "# Save to disk\n",
    "import pickle\n",
    "pickle.dump(clf, open(\"my-model.dat\", \"w\"))\n",
    "\n",
    "# Load from disk\n",
    "pickle.load(open(\"my-model.dat\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## TMVA in Python \n",
    "\n",
    "* <a href=\"https://github.com/yandex/rep\">Reproducible Experiment Platform</a> (REP, from Yandex) provides Python wrappers for TMVA.\n",
    "* It exposes TMVA through a scikit-learn compatible API, hence making its usage possible within the Python ecosystem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)\n",
    "\n",
    "from rep.estimators.tmva import TMVAClassifier\n",
    "clf = TMVAClassifier(method=\"kBDT\", BoostType=\"Bagging\", NTrees=100, \n",
    "                     UseRandomisedTrees=True, UseNVars=4, nCuts=-1, \n",
    "                     PruneMethod=\"NoPruning\", MaxDepth=10000, MinNodeSize=10e-7)\n",
    "clf.fit(X_train, y_train)\n",
    "print(\"Training accuracy =\", clf.score(X_train, y_train)) \n",
    "print(\"Test accuracy =\", clf.score(X_test, y_test)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Advanced topics\n",
    "\n",
    "Scikit-Learn is more than training classifiers. It also covers:\n",
    "\n",
    "- Clustering\n",
    "- Matrix decomposition\n",
    "- Kernel Density Estimation\n",
    "- Outlier detection\n",
    "- Out-of-core learning\n",
    "- ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Further readings\n",
    "\n",
    "\n",
    "- Scikit-Learn <a href=\"http://scikit-learn.org/stable/documentation.html\">documentation</a>, <a href=\"http://scikit-learn.org/stable/auto_examples/index.html\">example gallery</a>\n",
    "- PyCon 2015 tutorial: Parts <a href=\"https://www.youtube.com/watch?v=L7R4HUQ-eQ0\">1</a> and <a href=\"https://www.youtube.com/watch?v=oGqGxvqA9-k\">2</a>\n",
    "- Complementary packages: <a href=\"http://statsmodels.sourceforge.net/\">statsmodel</a>, <a href=\"http://dan.iel.fm/emcee/current/\">emcee</a>, <a href=\"http://www.nltk.org/\">NLTK</a>, <a href=\"http://deeplearning.net/software/theano/\">Theano</a>, <a href=\"https://github.com/lisa-lab/pylearn2\">Pylearn2</a>, <a href=\"https://github.com/JasperSnoek/spearmint\">spearmint</a>, ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Summary\n",
    "\n",
    "- Scikit-Learn provides all essential tools for machine learning.\n",
    "- It integrates within a larger Python scientific ecosystem.\n",
    "- Try it for yourself!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "questions?"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}