04_Lab_4_Performance_evaluation_model_selection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Mining and Exploration [INFR11007]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 4: Performance evaluation and model selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab, we look at various performance metrics for classification.  We then turn our attention to cross-validation and hyper-parameter tuning. Finally, we touch on Bayesian optimisation for hyper-parameter tuning.\n",
    "\n",
    "As always, let's start by importing the basic packages and modules we will need :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Import required packages \n",
    "from __future__ import division, print_function # Imports from __future__ since we're running Python 2\n",
    "import os\n",
    "import numpy as np \n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "random_state = 10 # Ensure reproducible results\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification performance evaluation\n",
    "In this lab, we will look at the following classification metrics:\n",
    "* classification accuracy\n",
    "* logarithmic loss\n",
    "* confusion matrices\n",
    "\n",
    "As always, we will make use of the familiar [landsat satellite](https://archive.ics.uci.edu/ml/datasets/Statlog+%28Landsat+Satellite%29) dataset which is 36-dimensional and comprises 6 classes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Landsat sattelite data pre-processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1 ==========\n",
    "Load the `landsat_train.csv` dataset into a `pandas` DataFrame called  `landsat_train_full`, and the `landsat_test.csv` dataset into a DataFrame called `landsat_test`. Display the shapes of the two DataFrames. \n",
    "\n",
    "*Hint: The DataFrames should have 37 columns including the class labels, and 4435 and 2000 entries for the training and testing datasets, respectively.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2 ==========\n",
    "Load the dataset class names stored in `landsat_classes.csv'` into a dictionary. You are free to choose whatever method you wish. \n",
    "\n",
    "Replace the label numbers in both the `landsat_train_full` and `landsat_test` DataFrames with the corresponding class names.\n",
    "\n",
    "*Hint: If unsure, check out the provided solutions for Lab 3.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 3 ==========\n",
    "Store the training features, training labels, testing features and testing labels into numpy arrays `X_train_full`, `y_train_full`, `X_test`, and `y_test` respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hold-out validation\n",
    "We currently have two datasets, namely `X_train_full` and `X_test`. If we just wanted to train a simple classifier with default settings and evaluate performance on the test subset, then our current approach should be good enough. \n",
    "\n",
    "Even simple classifiers, however, have hyper-parameters that need to be carefully tuned. In order to do so, we need a separate validation subset of the data which should be different from the training set. We should never perform model (i.e. hyper-parameter) selection by using the test set, because if we do so, we won't be able to evaluate the generalisation of our model on unseen data (and yes, in case you are wondering, that means that by performing model selection we might in a way *overfit* on the validation set).\n",
    "\n",
    "The simplest approach we can follow is to split our data three-way, that is, have a training, a validation, and a test set. We are already given the test set, so all we need to do is to split our training set (which we have called `train_full`) into training and validation subsets.\n",
    "\n",
    "Thankfully, sklearn offers an implementation of this operation which is called [`train_test_split`](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html). By default, this function will shuffle our data before splitting it. The `test_size` input parameter indicates the relative size of the test set (which will be used as the validation set in our case) and the `random_state` parameter can be used to ensure we can get reproducible results if we call this function multiple times.\n",
    "\n",
    "Let's transform our full training set into two subsets called `train` and `val`. We will feed in our `X_train_full` and `y_train_full` arrays and call the new arrays `X_train`, `X_val`,  `y_train`, and `y_val`, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_val, y_train, y_val = train_test_split(X_train_full, y_train_full, \n",
    "                                                  test_size=0.33, random_state=random_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Feature standardisation\n",
    "\n",
    "We have already seen that [feature standardisation](http://scikit-learn.org/stable/modules/preprocessing.html) (i.e. transforming the data so that they have zero mean and unit variance) is important for some unsupervised dimensionality reduction methods. It turns out that it is also crucial for the training efficiency of many supervised algorithms, especially those that use some form of optimisation (e.g. logistic regression trained via gradient descent). \n",
    "\n",
    "Feature standardisation can hardly ever do any harm when it comes to training classifiers, so it is generally good practice to deploy it as a first step in our data processing pipeline.\n",
    "\n",
    "It is essential, however, to perform feature standardisation by using the training data only. If we use the whole dataset for estimating feature means and variances, we will have information leakage from the test set to the training set, and our results might be over-optimistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 4 ==========\n",
    "Create a [`StandardScaler`](http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html) instance and fit it by using the training features only (`X_train`). \n",
    "\n",
    "Then by using the object you just fit, standardise (i.e. call the `transform()` method) the training, validation and test input features and save the results into three new numpy arrays, `X_train_sc`, `X_val_sc` and `X_test_sc`.\n",
    "\n",
    "*Hint: If unsure how to perform this step, check out Lab 3.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian Naive Bayes classification\n",
    "Now we want to use a simple Gaussian Naive Bayes classifier to get a feel for the  performance baseline we can achieve on our validation set. Read about the [Naive Bayes classifier](https://en.wikipedia.org/wiki/Naive_Bayes_classifier) and the underlying assumption if you are not already familiar with it.\n",
    "\n",
    "We will make use of the `GaussianlNB` class in sklearn. **Check out the user guide [description](http://scikit-learn.org/stable/modules/naive_bayes.html#gaussian-naive-bayes) and [documentation](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html) to familiarise yourself with this class.**\n",
    "\n",
    "All classifier objects in sklearn implement a `fit()` and `predict()` method. The first learns the parameters of the model and the latter classifies inputs. For a Naive Bayes classifier, the [`fit()`](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html#sklearn.naive_bayes.GaussianNB.fit) method takes at least two input arguments `X` and `y`, where `X` are the input features and `y` are the labels associated with each example in the training dataset (i.e. targets). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's train our classifier by calling its `fit()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "gnb = GaussianNB().fit(X_train_sc,y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classification accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit-learn model objects have built in scoring methods. The default [`score()`](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html#sklearn.naive_bayes.GaussianNB.score) method for `GaussianNB` estimates classification accuracy. Alternatively, we can compute the prediction for the test data and make use of the [`accuracy_score()`](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html) function (that is in fact what the classifier's `score()` method does under the hood)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "# By using the predict() method and accuracy_score metric\n",
    "gnb_prediction = gnb.predict(X_val_sc)\n",
    "gnb_accuracy = accuracy_score(y_val, gnb_prediction) # The accuracy_score() function takes as inputs\n",
    "                                                 # the true labels and the predicted ones\n",
    "\n",
    "# By using the score() method\n",
    "gnb_accuracy_alt = gnb.score(X_val_sc, y_val) # The score() method takes as inputs \n",
    "                                              # the test input features and the associated (true) labels\n",
    "\n",
    "# Print results\n",
    "print(\"GNB classification accuracy on validation set (by using the accuracy_score() function): {:.3f}\"\n",
    "      .format(gnb_accuracy))\n",
    "print(\"GNB classification accuracy on validation set (by using the model's score() method): {:.3f}\"\n",
    "      .format(gnb_accuracy_alt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 5 [optional] ==========\n",
    "Write your own function for computing classification accuracy by taking as inputs the vector of the true labels `y_true` and the vector of predicted labels `y_pred`. Compare its outcome to the results from using the scikit-learn `accuracy_score` metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline classification [optional]\n",
    "\n",
    "How can we know if the performance of our classifier is vaguely good or bad? One simple idea is to try to compare its performance to baseline.\n",
    "\n",
    "### ========== Question 6 [optional] ==========\n",
    "\n",
    "What is the simplest classifier you can think of? (Try for a moment to forget everything you know about machine learning and think how you would classify test inputs in the simplest, and perhaps dumbest, way).\n",
    "\n",
    "Implement the baseline classifier of your choice and compute its classification accuracy on the validation set. Does the GNB model perform better than the baseline?\n",
    "\n",
    "This might seem as an unnecessary hassle at this point, but it is important to always check what the baseline performance level is for a given task. There are cases where you can spend hours optimising a classifier, only to find out later on, that its performance does not exceed the baseline (i.e. because the input features, for instance, are not informative of the labels)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 7 [optional] ==========\n",
    "\n",
    "It turns out that sklearn implements the [`DummyClassifier`](http://scikit-learn.org/stable/modules/generated/sklearn.dummy.DummyClassifier.html) class which offers sever choices for dummy classifiers (`strategy` parameter).  Is the baseline classifier you came up with in the previous question included in this class? If so, double-check that your estimate about baseline performance matches the one returned by the `DummyClassifier` in sklearn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.dummy import DummyClassifier\n",
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confusion matrix\n",
    "Scikit-learn also has a [`confusion_matrix`](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.confusion_matrix.html#sklearn.metrics.confusion_matrix) implementation which returns a numpy array (square matrix) of dimensionality `K`, where `K` is the number of classes (`6` in our case)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 8 ========== \n",
    "By using the prediction of the Gaussian Naive Bayes model, compute and display the confusion matrix on the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You are provided with the following function which uses  `seaborn's` [`heatmap`](http://seaborn.pydata.org/generated/seaborn.heatmap.html) function to visualise a confusion matrix. The parameter `normalize` can be used to  convert the rows of the confusion matrix  into normalised prediction scores/probabilities (i.e. the sum of each row will be equal to 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot confusion matrix by using seaborn heatmap function\n",
    "def plot_confusion_matrix(cm, normalize=False, classes=None, title='Confusion matrix'):\n",
    "    \"\"\"Plots a confusion matrix.\n",
    "    \n",
    "    If normalize is set to True, the rows of the confusion matrix are normalized so that they sum up to 1.\n",
    "    \n",
    "    \"\"\"\n",
    "    if normalize is True:\n",
    "        cm = cm/cm.sum(axis=1)[:, np.newaxis]\n",
    "        vmin, vmax = 0., 1.\n",
    "        fmt = '.2f'\n",
    "    else:\n",
    "        vmin, vmax = None, None\n",
    "        fmt = 'd'\n",
    "    if classes is not None:\n",
    "        sns.heatmap(cm, xticklabels=classes, yticklabels=classes, vmin=vmin, vmax=vmax, \n",
    "                    annot=True, annot_kws={\"fontsize\":9}, fmt=fmt)\n",
    "    else:\n",
    "        sns.heatmap(cm, vmin=0., vmax=1.)\n",
    "    plt.title(title)\n",
    "    plt.ylabel('True label')\n",
    "    plt.xlabel('Predicted label')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 9 ========== \n",
    "\n",
    "Use the provided `plot_confusion_matrix` to visualise the confusion matrix on the test dataset. Make two calls of the function, one with the `normalize` parameter set to `True`, and another one with the same parameter set to `False`. Use a single figure with two subplots. Inspect the outcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-entropy (logistic) loss\n",
    "Classification accuracy is not the only metric used to evaluate classification performance, and in many cases, it can be very misleading. One such case is when we deal with unbalanced datasets in binary classification tasks. In this case, we might want to use different metrics such as [precision](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.average_precision_score.html#sklearn.metrics.average_precision_score), [recall](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.recall_score.html#sklearn.metrics.recall_score), [F1-score](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.f1_score.html#sklearn.metrics.f1_score) and/or [area under the ROC curve](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.html). Unfortunately, we won't have time to cover all these in the current lab, but you are free to experiment with them if you are keen.\n",
    "\n",
    "Another particular case where the accuracy metric is not very useful, is when we care not only about the classification predictions, but also about the associated probability distributions. Sometimes, for instance, we might consider that it is not catastrophic to do a misclassification, as long the associated prediction probability is not too high. In other words, we might want to penalise wrong classifications that are made with high confidence more than misclassifications that are made with low confidence.\n",
    "\n",
    "One metric that takes into account the prediction probabilities is the [logarithmic loss or cross-entropy](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.log_loss.html).  This metric is very often used as the evaluation score in [data science competitions](https://www.kaggle.com/wiki/LogarithmicLoss). You might recognise this metric as the loss funcion used while training logistic regression models and neural networks.\n",
    "\n",
    "### ========== Question 10 ========== \n",
    "\n",
    "By using the [log-loss](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.log_loss.html) implementation in sklearn, compute the log-loss scores for the Dummy and Gaussian Naive Bayes classifers from Questions 6 and 8."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import log_loss\n",
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classifier comparison\n",
    "Now we want to compare the performance of different types of classifiers on the validation set. This is a common first step to decide which classifier might be most suitable for the task at hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 11 ========== \n",
    "\n",
    "Compare the classification accuracy of various types of classifiers of your choice (at least three different types) on the landsat dataset. Use the training set for training classifiers  and the validation set to evaluate performance. \n",
    "\n",
    "Log and print the classification accuracy and cross-entropy (aka logarithmic loss) scores for each classifier. \n",
    "\n",
    "*Reminder: Make sure you make use of the standardised versions of the data, i.e. `X_train_sc` and `X_val_sc`. *\n",
    "\n",
    "*Hint:  You might find this `sklearn` [tutorial](http://scikit-learn.org/stable/auto_examples/classification/plot_classifier_comparison.html) very useful.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 12 ========== \n",
    "\n",
    "Which classifier seems to perform best on the validation set? \n",
    "\n",
    "Are your observations as expected, or do you find them rather suprising? \n",
    "\n",
    "Would you trust the results of this  comparison?  If not, explain why."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hyper-parameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to optimise the hyper-parameters of classifiers we should never do so by using the test set. If we compare the performance of various classifiers -after we have tuned their hyper-parameters- by using a validation set, we should equally not use the same  set for parameter optimisation. \n",
    "\n",
    "Instead, we should use two independent procedures for hyper-parameter tuning and classifier comparison. For instance, we could split again our training set into two subsets and use the former to train classifiers under various parameter configurations and the latter to assess performance. Then by picking the parameter settings which gave optimal results for each classifier, we can compare the performance of the best-performing models on our original validation set. \n",
    "\n",
    "Alternatively, we can use [K-fold cross-validation](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.KFold.html) on the training set for setting hyper-parameters; this approach will yield slightly more reliable results, since hyper-parameter selection will be based on averaging across a few runs. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sklearn offers a very convenient tool for fitting model hyperparameters called [`GridSearchCV`](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html). Spend a few minutes reading the documentation of this class and make sure you understand how it works.\n",
    "\n",
    "A `GridSearchCV` classifier is constructed by defining the following parameters, among others:\n",
    "* `estimator`: this is another (base) classifier which has been constructed but not fitted yet.\n",
    "* `param_grid`: this is a dictionary, where the different hyper-parameters to be optimised are defined, as well as the search space for each hyper-parameter, which can be either discrete or continuous.\n",
    "* `scoring`: this is a string defining the objective function (i.e. scoring method) to be used for hyper-parameter optimisation. Some options are `accuracy` for classification accuracy, `neg_log_loss` for log-loss, or user-defined metrics. A list of all the available metrics can be found [here](http://scikit-learn.org/stable/modules/model_evaluation.html).\n",
    "* `cv`: cross-validation to be used for optimising the hyper-parameters. This can be either an integer `K` or a sklearn cross-validation generator. If an integer is provided, then `KFold` cross-validation will be used. By default sklearn will use 3-fold CV."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now make a first attempt to optimise the regularisation parameter `C` of an [`SVM`](http://scikit-learn.org/stable/modules/generated/sklearn.svm.SVC.html) classifier with an RBF kernel. \n",
    "\n",
    "Spend a few moments to understand what the following piece of code does. You can hopefully  see that after we have fitted a `GridSearchCV` classifier, we can access the best scores achieved and best-scoring parameter configurations by looking at the `best_params_` and `best_score_` attributes, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import KFold\n",
    "\n",
    "cv = KFold(n_splits=3, shuffle=True, random_state=random_state)\n",
    "svc = SVC(kernel='rbf')\n",
    "parameters = {'C': np.logspace(-3,3,7)}\n",
    "svc_clf = GridSearchCV(estimator=svc, cv=cv, param_grid=parameters, scoring='accuracy')\n",
    "svc_clf.fit(X_train_sc, y_train)\n",
    "print(\"Best setting of C parameter for SVC with RBF kernel: {}\".format(svc_clf.best_params_[\"C\"]))\n",
    "print(\"Best cross-validated score: {:.3f}\".\n",
    "      format(svc_clf.best_score_))\n",
    "print(\"Classification accuracy on validation set: {:.3f}\".format(svc_clf.score(X_val_sc,y_val)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### ========== Question 13 ========== \n",
    "\n",
    "By adapting the above code provided, optimise both the regularisation parameter `C` and the kernel coefficient `gamma` of an SVM classifier with RBF kernel. \n",
    "\n",
    "For `C`, you can use the previous grid, and for `gamma` you can use a logarithmic range between $10^{-4}$ to $10^{1}$.\n",
    "\n",
    "Print the best scoring parameter configuration and the classification accuracy score on the validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 14 ========== \n",
    "\n",
    "How does the performance on the validation set compare to that achieved previously with other classifers? \n",
    "\n",
    "If you did not want to use the provided `GridSearchCV` and `KFold` modules from sklearn but wanted to write your own code from scratch, how many *for-loops* would you have included in your code?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 15 [Optional] ========== \n",
    "\n",
    "Optimise the number of hidden units and regularisation constant `alpha` of an [MLP classifier](http://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.html) with one hidden layer. Use the `neg_log_loss` as the scoring function, which will perform hyper-parameter optimisation by minimising the log-loss (i.e. cross-entropy). Use default settings for the other parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 16 [Optional] ========== \n",
    "\n",
    "Why is the classification accuracy not better than the one achieved in Question 11? Did we do something wrong?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayesian optimisation is a strategy for optimising black-box functions without the need for computing derivatives. Bayesian optimisation methods can prove very useful for optimising the hyper-parameters of machine learning models in arbitrary search spaces, for both continuous and discrete-valued hyper-parameters.\n",
    "\n",
    "In this final section, we will use the [`scikit-optimize`](https://scikit-optimize.github.io) package to tune the `C` and `gamma` hyper-parameters of an SVM classifier with RBF kernel (see Question 13). We will use a Gaussian Process prior over the average cross-validated classification accuracy score, and will try to maximise this metric by tuning the SVM hyper-parameters. \n",
    "\n",
    "Spend a few moments trying to understand what the following piece of code does. If something is not clear, you can consult the skopt [documentation and examples](https://scikit-optimize.github.io)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "cv = KFold(n_splits=3, shuffle=True, random_state=random_state) # We want to use 3-fold CV to compute the loss in each iteration\n",
    "svc_clf = SVC(kernel='rbf') # Constructor without defining hyper-parameters just yet\n",
    "\n",
    "def objective_svc(params): # Here we define the metric we want to minimise\n",
    "    C, gamma = params\n",
    "    svc_clf.set_params(C=C, \n",
    "                      gamma=gamma) # Set the parameters of the clf\n",
    "    \n",
    "    return -np.mean(cross_val_score(svc_clf, X_train_sc, y_train, cv=cv, n_jobs=-1,\n",
    "                                    scoring=\"accuracy\")) # We want to maximise average accuracy, i.e. minimise minus average accuracy\n",
    "\n",
    "# Search space for the two parameters\n",
    "space  = [(10**-3, 10**3, \"uniform\"), # C\n",
    "          (10**-4, 10**1, \"uniform\")] # gamma\n",
    "\n",
    "# Initial values (optional)\n",
    "x0 = [1, 10**-2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following piece of code uses the above defined cost function (`objective_scv`) and the skopt `gp_minimise` function to select the best hyper-parameter configuration in the defined search space. The input parameters are as follows:\n",
    "* `func`: function that we wish to minimise\n",
    "* `dimensions`: the search space for the hyper-parameters\n",
    "* `x0`: inital values for the hyper-parameters\n",
    "* `n_calls`: number of times the function will be evaluated\n",
    "* `random_state`: random seed\n",
    "* `n_random_starts`: before we start modelling the optimised function with a GP Regression model, we want to try a few random choices for the hyper-parameters.\n",
    "* `kappa`: trade-off between [exploration vs. exploitation](https://en.wikipedia.org/wiki/Multi-armed_bandit)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skopt import gp_minimize\n",
    "res_gp = gp_minimize(func=objective_svc, dimensions=space, x0=x0, \n",
    "                     n_calls=25, random_state=random_state, n_random_starts=5, kappa=1.9)\n",
    "print(\"Best score with Bayesian optimisation: {:.3f}\".format(-res_gp.fun))\n",
    "print(\"Best parameters with Bayesian optimisation:\\nC: {}\\ngamma: {}\"\n",
    "      .format(res_gp.x[0],res_gp.x[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using the skopt `plot_convergence` function, we can plot the objective function (i.e. evaluated cross-validated accuracy score) against the number of iterations of the Bayesian optimisation process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skopt.plots import plot_convergence\n",
    "plt.figure(figsize=(5,2))\n",
    "plot_convergence(res_gp)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to clarify here, that the plot above shows the minimum value of the optimised function achieved at step $n$. Nevertheless, this does not mean that at a certain step the objective function should always decrease as compared to the previous step (as is the case with e.g. with gradient descent). \n",
    "\n",
    "If we wish to plot the value of the evaluated function $f \\left (x \\right)$ vs. the iteration number $n$, we can do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,2))\n",
    "plt.plot(res_gp.func_vals)\n",
    "plt.scatter(range(len(res_gp.func_vals)), res_gp.func_vals)\n",
    "plt.ylabel(r'$f(x)$')\n",
    "plt.xlabel('Number of calls $n$')\n",
    "plt.xlim([0, len(res_gp.func_vals)])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's compute the classification accuracy on the validation set by using the selected hyper-parameter configuration. Since we are not performing any further comparisons or hyper-parameter optimisation, we can finally use the test set to report the performance of our selected method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "svc_opt = SVC(kernel='rbf',\n",
    "             C=res_gp.x[0],\n",
    "             gamma=res_gp.x[1]).fit(X_train_sc,y_train)\n",
    "print(\"Classification accuracy on validation set: {:.3f}\".format(accuracy_score(y_val, svc_opt.predict(X_val_sc))))\n",
    "print(\"Classification accuracy on test set: {:.3f}\".format(accuracy_score(y_test, svc_opt.predict(X_test_sc))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a moment to reflect on what we have achieved here. By running only 25 function evaluations (i.e. iterations), we have achieved 0.916 accuracy on the validation set. This is only 0.1%  smaller than the best score achieved with grid search optimisation, where we had to evaluate our cost function 42 times (since we used a 7 X 6 grid). If we allow the GP optimiser to run for more iterations, we might even get higher performance than we can achieve with grid-search, since we are not restricting our hyper-parameters to discrete values, as opposed to grid search optimisation.\n",
    "\n",
    "Finally, it is worth noting that the classification performance scores are slightly different for the validation and test sets. In this part we have not used the validation set to perform hyper-parameter optimisation (remember we used K-fold CV within the training set to do so), hence, we can view `y_val` and `y_test` two as two separate test sets. The computed classification accuracy scores are only approximations of the performance of the classifier on unseen data. The observed difference in performance for the two sets is due to the fluctuation of the approximation when using finite size datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 17 [Optional] ========== \n",
    "\n",
    "Use the skopt package to tune the hyper-parameters of an [MLPClassifier](http://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.html) with one hidden layer (i.e. number of hidden units and `alpha`) via Bayesian optimisation. \n",
    "\n",
    "Finally, report the accuracy of the best parameter configuration on the validation and test sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:dme_env]",
   "language": "python",
   "name": "conda-env-dme_env-py"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}