04_Lab_2_Decision_trees_linear_regression_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory applied machine learning (INFR10069) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 2: Decision trees and linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first part of this lab we perform Decision trees classification on the [German credit](http://archive.ics.uci.edu/ml/datasets/Statlog+%28German+Credit+Data%29) data set. In the second part we learn how to train simple linear regression model by using the [CPU performance](https://archive.ics.uci.edu/ml/datasets/Computer+Hardware) data set. Both datasets (`credit.csv` and `cpu.csv`) are located at the `./datasets` directory which is adjacent to this file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we're going to intall a new package to visulaize trees. Run the following commands in your terminal to install (remember to exclude 'source' if you're on windows):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{bash}\n",
    "source activate iaml\n",
    "conda install python-graphviz\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this causes issue, it's not essential. Simply exclude commands involving graphviz."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's import the packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import 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",
    "import graphviz\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn import tree\n",
    "from sklearn.tree import DecisionTreeClassifier, export_graphviz\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix, r2_score\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Decision Trees\n",
    "One of the great advantages of decision trees is their interpretability. The rules learnt for classification are easy for a person to follow, unlike the opaque \"black box\" of many other methods, such as neural networks. We demonstrate the utility of this using a German credit data set. You can read a description of this dataset at the [UCI site](http://archive.ics.uci.edu/ml/datasets/Statlog+%28German+Credit+Data%29). The task is to predict whether a loan approval is good or bad credit risk based on 20 attributes. We've simplified the data set somewhat, particularly making attribute names and values more meaningful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Download the dataset and save it in a directory called `datasets` in the same folder that your notebooks live. Alternatively, you can save the dataset in any folder you wish and modify the `data_path` variable below accordingly. We will load our data into a pandas DataFrame structure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load the dataset\n",
    "data_path = os.path.join(os.getcwd(), 'datasets', 'credit.csv')\n",
    "credit = pd.read_csv(data_path, delimiter = ',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.1 ==========\n",
    "Display the number of data points and attributes in the dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "print('Number of samples: {}, number of attributes: {}'.format(credit.shape[0], credit.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.2 ==========\n",
    "Get a feeling of the data by using pandas `describe()` method. Be careful - there is a mixture of numeric and categorical data..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "credit.describe(include=[np.number])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "credit.describe(include=['O'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.3 ==========\n",
    "Display the first 10 data points of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "pd.set_option('display.max_columns', None)\n",
    "credit.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.4 ==========\n",
    "When presented with a dataset, it is usually a good idea to visualise it first. By using seaborn's [pairplot](https://seaborn.github.io/generated/seaborn.pairplot.html?highlight=pairplot#seaborn.pairplot) function, try visualising a scatter plot of the `Age` and `Duration` variables. You can use the `Approve` variable as the `hue` parameter to visualise results separately for each class. Do you notice anything unusual?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "g = sns.pairplot(data=credit, vars=['Duration', 'Age'], hue='Approve', size=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "There is a data point with negative age. This doesn't make sense whatsoever and this data point has been clearly corrupted with noise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.5 ==========\n",
    "In the previous point you should have found a data point, which seems to be corrupted, as some of its values are nonsensical. Even a single point like this can significantly affect the performance of a classifier. How do you think it would affect Decision trees? How about Naive Bayes? A good way to check this is to test the performance of each classifier before and after removing this datapoint."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "Naive Bayes fits Gaussian distributions and thus is very sensitive to outliers. Decision trees is expected to be less affected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.6 ==========\n",
    "Now we want to remove this instance from the dataset by using a filter. We want to remove all instances, where the age of an applicant is lower than 0 years, as this suggests that the instance is corrupted. Use logical indexing to get rid of these instances without creating a new dataframe. Display the number of data points after any outliers have been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "credit = credit[credit['Age']>0]\n",
    "print('Number of data points after removal of outliers: {}'.format(credit.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.7 ==========\n",
    "\n",
    "You might have noticed that most of the attributes in the dataset are in fact discrete. Now we want to know which variables exactly are discrete (both categorical and numerical, look [here](http://stats.stackexchange.com/questions/206/what-is-the-difference-between-discrete-data-and-continuous-data) if you are unsure about the difference) and which are continuous variables. In order to do so, we will inspect the number of possible values that each attribute can take. \n",
    "\n",
    "Display the number of values each attributes takes in the dataset. *Hint: As a first step, you want to loop over the columns of the DataFrame. Then you might find the numpy `unique` function quite useful.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "for column in credit:\n",
    "    unique_vals = np.unique(credit[column])\n",
    "    nr_vals = len(unique_vals)\n",
    "    if nr_vals < 11:\n",
    "        print('Number of values for attribute {}: {} -- {}'.format(column, nr_vals, unique_vals))\n",
    "    else:\n",
    "        print('Number of values for attribute {}: {}'.format(column, nr_vals)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems like the variables `Duration`, `CreditAmount` and `Age` are continuous and all the rest are discrete. The discrete variables are not in a very convenient format though. Ideally we would want the discrete attributes to take values between `0` and `n_values-1`. Scikit-learn has a handy [`LabelEncoder`](http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.LabelEncoder.html#sklearn.preprocessing.LabelEncoder) implementation which can do that for us. You are encouraged to read its documentation.\n",
    "\n",
    "Now we will create a new DataFrame called `credit_clean` and convert all the discrete variables from `credit` by using a `LabelEncoder`. Remember, we want to change the discrete variables only, so we will have to exclude the `CreditAmount`, `Age` and `Duration` attributes. Also, we don't really mind if the target variable is categorical, so we won't be transforming the `Approve` attribute either. Execute the cell below and make sure you understand what it does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "le_dict = dict() # Initialise an empty dictionary to keep all LabelEncoders\n",
    "credit_clean = credit.copy(deep=True) # Make a copy of the DataFrame\n",
    "# Loop over attributes by excluding the ones that are continuous and the target variable\n",
    "for column in credit_clean.drop(['CreditAmount', 'Age', 'Duration', 'Approve'], axis=1):  \n",
    "    le = LabelEncoder().fit(credit[column]) # Initialise the LabelEncoder and fit\n",
    "    credit_clean[column] = le.transform(credit[column]) # Transform data and save in credit_clean DataFrame\n",
    "    le_dict[column] = le # Store the LabelEncdoer in dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.8 ==========\n",
    "Display the first 10 data points of the clean data. Does it look like what you expected?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "credit_clean.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.9 ==========\n",
    "Store the input features (i.e. attributes) into a matrix `X` and the target variable (`Approve`) into a vector `y`. Remember to not include the target variable into `X`. *Hint: You can either use pandas `as_matrix()` or `values`.* \n",
    "\n",
    "Display the shapes of `X` and `y`. Confirm that you have 20 input features, one target variable and 1000 data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "X = credit_clean.drop('Approve', axis=1).as_matrix() # Input features (attributes)\n",
    "y = credit_clean['Approve'].as_matrix() # Target vector\n",
    "print('X shape: {}'.format(np.shape(X)))\n",
    "print('y shape: {}'.format(np.shape(y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hold-out validation\n",
    "In the next step we will be using a Decision Tree classifier model. To get an accurate estimate of the model's classification performance we will use hold-out validation. Familiriase yourself with the logic behind [`train_test_split CV`](http://scikit-learn.org/stable/modules/cross_validation.html#cross-validation-evaluating-estimator-performance) (also called `Hold-out` validation) and [how it is used](http://scikit-learn.org/0.16/modules/generated/sklearn.cross_validation.train_test_split.html) in `Scikit-learn`. Execute the cell below to create your training/testing sets by assigning 10% of the data to the test set (and convince yourself you understand what is going on)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size = 0.9, test_size=0.1, random_state=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.10 ==========\n",
    "Confirm that `X_train` and `X_test` matrices are subsets of `X` by displaying the number of rows in the three matrices (no need to make use of set theory)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "print('Number of instances in X: {}'.format(np.shape(X)[0]))\n",
    "print('Number of instances in X_train: {}'.format(X_train.shape[0]))\n",
    "print('Number of instances in X_test: {}'.format(X_test.shape[0]))\n",
    "print('Number of instances in X_train and X_test together: {}'.format(X_train.shape[0] + X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.11 ==========\n",
    "Now we will train a Decision Tree classifier on the training data. Read about [Decision Tree classifiers](http://scikit-learn.org/stable/modules/tree.html) in `Scikit-learn` and how they are [used](http://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html#sklearn.tree.DecisionTreeClassifier). \n",
    "Create a `DecisionTreeClassifier` instance and train it by using training data only (i.e. `X_train` and `y_tain`). Set the `criterion` attribute to `entropy` in order to measure the quality of splits by using the information gain. Use the default settings for the rest of parameters. By default, trees are grown to full depth; this means that very fine splits are made involving very few data points. Not only does this make the trees hard to visualise (they'll be deep), but also we could be overfitting the data. For now, we arbitrarily choose a depth of 2 for our tree, but this is a parameter we could tune."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "dt = DecisionTreeClassifier(criterion='entropy', max_depth=2, random_state=1337)\n",
    "dt.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have mentioned in the class that decision trees have the advantage of being interpretable by humans. Now we visualise the decision tree we have just trained. Scikit-learn can export the tree in a `.dot` format. \n",
    "\n",
    "An alternative way to visualise the tree is to open the output .dot file with an editor such as [this online .dot renderer](http://dreampuf.github.io/GraphvizOnline/). You can copy and paste the dot file in and view online (you can double click on the tree once it has been produced to view it in full screen).\n",
    "\n",
    "This code below should work if you created a decision tree classifier object called `dt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dot_data = tree.export_graphviz(dt, out_file=None, \n",
    "    feature_names=credit_clean.drop('Approve', axis=1).columns,  \n",
    "    class_names=credit_clean['Approve'].unique(),  \n",
    "    filled=True, rounded=True,  \n",
    "    special_characters=True)\n",
    "graph = graphviz.Source(dot_data)\n",
    "graph\n",
    "\n",
    "## To create a file do this:\n",
    "# with open(\"tree.dot\", 'w') as f:\n",
    "#     f = export_graphviz(dt, out_file=f,\n",
    "#                         feature_names=credit_clean.drop('Approve', axis=1).columns,  \n",
    "#                         class_names=credit_clean['Approve'].unique(),  \n",
    "#                         filled=True, rounded=True,  \n",
    "#                         special_characters=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.12 ==========\n",
    "Inspect the tree. Describe what it shows. Which attribute yields the highest information gain and what is its critical value?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "Below is a png of the tree produced. It shows that we will classify a loan as `bad` if the `CheckingAccount in (2, 3)` i.e. is '>=200' or 'none', and even worse if `OtherPlans in (1, 2)` i.e. is 'none' or 'stores'.\n",
    "\n",
    "The attribute with highest information gain is the one at the top of the tree - `CheckingAccount`. Its critical value is 1.5. To trace back the original values we have to make use of the associated LabelEncoder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(le_dict['CheckingAccount'].classes_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The critical value for the `CheckingAccount` attribute is thus <200."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(filename='tree.png') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "N.B. You should query the validity of this tree. The categories here are being treated as numeric, i.e. ordinal. Is 'none' really greater than '>=200' in the `CheckingAccount` attribute...no not really! This is a drawback of the CART method implemented by scikit learn. One way to force the tree to think of these categories as separate is to use a `OnehotEncoder` but, in practice, because the tree can separate the classes it needs by using multiple splitting conditions, this isn't usually a problem. An example of where it is a problem would be if an important condition was if `CheckingAccount` was  '<0' or 'none'. This tree **can't** make a single branch that contains this condition - if you used a `OnehotEncoder` it could."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.13 ==========\n",
    "Tree-based estimators (i.e. decision trees and random forests) can be used to compute feature importances. The importance of a feature is computed as the (normalized) total reduction of entropy (or other used `criterion`) brought by that feature. Find the relevant attribute of the classifier you just trained and display feature importances along with their names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "for i, column in enumerate(credit_clean.drop('Approve', axis=1)):\n",
    "    print('Importance of feature {}:, {:.3f}'.format(column, dt.feature_importances_[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.14 ==========\n",
    "Now we want to evaluate the performance of the classifier on unseen data. Use the trained model to predict the target variables for the test data set. Display the classification accuracy for both the training and test data sets. What do you observe? Are you surprised by the results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "prediction_train = dt.predict(X=X_train)\n",
    "prediction_test = dt.predict(X=X_test)\n",
    "print('Classification accuracy on training set: {:.3f}'.format(accuracy_score(y_train,prediction_train)))\n",
    "print('Classification accuracy on test set: {:.3f}'.format(accuracy_score(y_test,prediction_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "The decision tree classifier is not overfitting. The classification accuracy is similar for training and test. The decision tree has provided a very simple way to interpret the data - splitting it into 4 bins and applying a class to each bin. Clearly this is an oversimplification: even in the training data, the leaves of the tree contain many examples of the 'incorrect' class.\n",
    "\n",
    "**However, the very astute will make an observation:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(8,4))\n",
    "sns.countplot(y_train, ax=ax1)\n",
    "ax1.set_title('Training set', size=12)\n",
    "ax1.set_xlabel(' ')\n",
    "sns.countplot(y_test, ax=ax2)\n",
    "ax2.set_xlabel(' ')\n",
    "ax2.set_ylabel('')\n",
    "ax2.set_title('Validation set', size=12)\n",
    "fig.suptitle('Target distribution', size=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dummy classifier of predicting everything as 'good' **will beat this model (75% accuracy)!** It is worse than the baseline. **Always compare your models with very simple baselines**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.15 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit another `DecisionTreeClassifier` but this time grow it to full depth (i.e. remove the max_depth condition. Display the classification accuracy for training and test data as above. Again, what do you observe and are you surprised?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "dt_fulldepth = DecisionTreeClassifier(criterion='entropy', random_state=1337)\n",
    "dt_fulldepth.fit(X_train, y_train)\n",
    "prediction_train2 = dt_fulldepth.predict(X=X_train)\n",
    "prediction_test2 = dt_fulldepth.predict(X=X_test)\n",
    "print('Classification accuracy on training set: {:.3f}'.format(accuracy_score(y_train,prediction_train2)))\n",
    "print('Classification accuracy on test set: {:.3f}'.format(accuracy_score(y_test,prediction_test2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "This is a classic case of overfitting. The tree is allowed to grow to full depth and perfectly fit the training data. The result is that the test accuracy is much lower than the train; in fact it is less accurate than the very simple depth 2 tree! It's important to note that, when grown to full depth, decision tree classifiers can perfectly separate the training data (if no two datapoints are identical but have different classes). Performance on unseen data is likely to be much lower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.16 ==========\n",
    "By using seaborn's heatmap function, plot the normalised confusion matrices for both the training and test data sets **for the max_depth=2 decision tree from question 1.11**. Make sure you label axes appropriately. *Hint: You can make use of the `plot_confusion_matrix` function introduced in a previous lab...*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "def plot_confusion_matrix(cm, classes=None, title='Confusion matrix'):\n",
    "    \"\"\"Plots a confusion matrix.\"\"\"\n",
    "    if classes is not None:\n",
    "        sns.heatmap(cm, xticklabels=classes, yticklabels=classes, vmin=0., vmax=1., annot=True)\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_train, prediction_train)\n",
    "cm_norm = cm/cm.sum(axis=1)[:, np.newaxis]\n",
    "plt.figure()\n",
    "plot_confusion_matrix(cm_norm, classes=dt.classes_, title='Training confusion')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_test, prediction_test)\n",
    "cm_norm = cm/cm.sum(axis=1)[:, np.newaxis]\n",
    "plt.figure()\n",
    "plot_confusion_matrix(cm_norm, classes=dt.classes_, title='Test confusion')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**N.B. it will be obvious if you have plotted the full depth decision tree as the training confusion matrix will be the identity**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.17 ==========\n",
    "\n",
    "Finally we will create a [`Random decision forest`](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html) classifier and compare the performance of this classifier to that of the decision tree. The random decision forest is an ensemble classifier that consists of many decision trees and outputs the class that is the mode of the class's output by individual trees. Start with `n_estimators = 100`, use the `entropy` criterion and the same train/test split as before. Plot the classification accuracy of the random forest model on the test set and show the confusion matrix. How does the random decision forest compare performance wise to the decision tree?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "rf = RandomForestClassifier(n_estimators=100, criterion='entropy')\n",
    "rf.fit(X_train, y_train)\n",
    "prediction_test = rf.predict(X=X_test)\n",
    "print('Classification accuracy on test set: {:.3f}'.format(accuracy_score(y_test,prediction_test)))\n",
    "cm = confusion_matrix(y_test, prediction_test)\n",
    "cm_norm = cm/cm.sum(axis=1)[:, np.newaxis]\n",
    "plt.figure()\n",
    "plot_confusion_matrix(cm_norm, classes=rf.classes_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.18 ==========\n",
    "How high can you get the performance of the classifier by changing the max depth of the trees (`max_depth`), or the `max_features` parameters? Try a few values just to get a look. *Don't do a grid search or anything in-depth, just get a feel*. Try the same settings twice...do you get the same accuracy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "from itertools import product\n",
    "n_estimators = 500\n",
    "max_features = [1, 'sqrt', 'log2']\n",
    "max_depths = [None, 2, 5, 10]\n",
    "for f, d in product(max_features, max_depths): # with product we can iterate through all possible combinations\n",
    "    rf = RandomForestClassifier(n_estimators=n_estimators, \n",
    "                                criterion='entropy', \n",
    "                                max_features=f, \n",
    "                                max_depth=d, \n",
    "                                n_jobs=2,\n",
    "                                random_state=1337)\n",
    "    rf.fit(X_train, y_train)\n",
    "    prediction_test = rf.predict(X=X_test)\n",
    "    print('Classification accuracy on test set with max features = {} and max_depth = {}: {:.3f}'.format(f, d, accuracy_score(y_test,prediction_test)))\n",
    "    cm = confusion_matrix(y_test, prediction_test)\n",
    "    cm_norm = cm/cm.sum(axis=1)[:, np.newaxis]\n",
    "    plt.figure()\n",
    "    plot_confusion_matrix(cm_norm, classes=rf.classes_, title='Confusion matrix accuracy on test set with max features = {} and max_depth = {}: {:.3f}'.format(f, d, accuracy_score(y_test,prediction_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "N.B. Observing these confusion matrices you'll see something very important - for some configurations, the Random Forest **always predicts the majority class**. This highlights (again) the importance of always checking performance against a dummy classifier!!!\n",
    "\n",
    "Additionally, if you want to reproduce your results, you must set the random seed (you can do this with the `random_state` argument). Random forests are...random!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.19 ==========\n",
    "Compare the feature importances as estimated with the decision tree and random forest classifiers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "rf = RandomForestClassifier(n_estimators=500, \n",
    "                            criterion='entropy', \n",
    "                            max_features=1, \n",
    "                            max_depth=10, \n",
    "                            n_jobs=2)\n",
    "rf.fit(X_train, y_train)\n",
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "xx = np.arange(len(dt.feature_importances_))\n",
    "ax1.bar(xx, dt.feature_importances_)\n",
    "ax1.set_xticks(xx)\n",
    "ax1.set_xticklabels(credit_clean.drop('Approve', axis=1).columns, rotation='vertical')\n",
    "ax1.set_title('Decision tree depth 2 importances')\n",
    "ax2.bar(xx, rf.feature_importances_)\n",
    "ax2.set_xticks(xx)\n",
    "ax2.set_xticklabels(credit_clean.drop('Approve', axis=1).columns, rotation='vertical')\n",
    "ax2.set_title('Random forest max features=1 , max depth=10')\n",
    "\n",
    "for i, column in enumerate(credit_clean.drop('Approve', axis=1)):\n",
    "    print('Importance of feature {}, DT: {:.3f}, RF: {:.3f}'.format(column, dt.feature_importances_[i], rf.feature_importances_[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Linear regression\n",
    "In the second part of the lab we use the [CPU performance](https://archive.ics.uci.edu/ml/datasets/Computer+Hardware) dataset for a simple regression task. Famliarise yourself with the dataset before moving on to the next step. Note that the version we will be using is missing the `Model Name` and `PRP` attributes. Our task will be to use the remaining attributes to predict `ERP` values.\n",
    "\n",
    "Download the dataset and save it in a directory called `datasets` in the same folder that your notebooks live. Alternatively, you can save the dataset in any folder you wish and modify the `data_path` variable below accordingly. We will load our data into a pandas DataFrame structure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load the dataset\n",
    "data_path = os.path.join(os.getcwd(), 'datasets', 'cpu.csv')\n",
    "cpu = pd.read_csv(data_path, delimiter = ',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.1 ==========\n",
    "Display the number of data points and attributes in the dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "print('Number of samples: {}, number of attributes: {}'.format(cpu.shape[0], cpu.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.2 ==========\n",
    "Get a feeling of the data by using pandas `describe()` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "cpu.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.3 ==========\n",
    "Display the first 10 data points of the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "cpu.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.4 ========== \n",
    "You might have noticed that the `vendor` attribute is categorical. This will give problems when using a linear regression model. For now we can simply remove this attribute. Create a new DataFrame called `cpu_clean` by copying `cpu` but omit the `vendor` attribute. Display the number of samples and attributes in the clean dataset as a sanity check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "cpu_clean = cpu.copy(deep=True)\n",
    "cpu_clean=cpu_clean.drop('vendor', axis=1)\n",
    "print('Clean dataset, number of samples: {}, number of attributes: {}'.format(cpu_clean.shape[0], cpu_clean.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.5 ==========\n",
    "Now -as always- we want to perform some exploratory data analysis. Remember that our task is to predict `ERP` values, so it's a good idea to inspect individual scatter plots of the target variable (`ERP`) against our input features. For this purpose we will use once again seaborn's pairplot implementation.\n",
    "\n",
    "Create a series of [pairplots](https://seaborn.github.io/generated/seaborn.pairplot.html?highlight=pairplot#seaborn.pairplot) showing the pairwise relationship of `ERP` and the remaining attributes in the data. You can do so by specifying explicitly the `x_vars` and `y_vars` input arguments in the pairplot. *Hnt: Your final plot will consist of 6 subplots each contaning a scatter plot.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "g2 = sns.pairplot(data=cpu, x_vars=['MYCT', 'MMIN', 'MMAX', 'CACH', 'CHMIN', 'CHMAX'], y_vars = 'ERP', size=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.6 ==========\n",
    "* Do you think that ERP should be at least partially predictable from the input attributes?\n",
    "* Do any attributes exhibit significant correlations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "Yes, it seems that some input features are positively correlated with the target variable so it should be possible to partially predict it. The `MMIN` and `MMAX` features appear to exhibit the strongest correlations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.7 ==========\n",
    "Now we have a feel for the data and we will try fitting a simple linear regression model. Similarly to what we did in the first part of the lab, we want to use cross-validation to evaluate the goodness of the fit.\n",
    "\n",
    "By using the `cpu_clean` dataset extract the raw values for the input features and the target variable and store them in two matrices, called `X` and `y` respectively. \n",
    "\n",
    "Then, split the dataset into training and testing sets by using a 75%-25% split (training/testing).\n",
    "\n",
    "Display the shapes of all matrices involved and double-check that all dimensionalities appear to be as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "X = cpu_clean.drop('ERP', axis=1).as_matrix()\n",
    "y = cpu_clean['ERP'].as_matrix()\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.75, test_size=0.25, random_state=0)\n",
    "\n",
    "print('Number of instances in X: {}'.format(np.shape(X)[0]))\n",
    "print('Number of instances in X_train: {}'.format(X_train.shape[0]))\n",
    "print('Number of instances in X_test: {}'.format(X_test.shape[0]))\n",
    "print('Number of instances in X_train and X_test together: {}'.format(X_train.shape[0] + X_test.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.8 ==========\n",
    "Fit a simple linear regressor by using the [`LinearRegression`](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html) model in Scikit-learn. Report the training accuracy by using the `score` attribute. What does this represent?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "lm = LinearRegression(fit_intercept=True, normalize=True, copy_X=True)\n",
    "lm.fit(X_train, y_train)\n",
    "print('Training accuracy: {:.3f}'.format(lm.score(X_train, y_train)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*\n",
    "\n",
    "The accuracy score represents the coefficient of determination ($R^2$). This is at max 1, but can be negative. It will be 0 if you predict the mean of y for all observations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.9 ==========\n",
    "Now report the testing accuracy by using the `score` attribute of the regressor as well as the `r2_score` metric. Confirm that these two yield identical results.\n",
    "\n",
    "How does the accuracy compare to the one reported on the training dataset? Do you think that your model does well on generalising on unseen data? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "print('Testing accuracy by using score function: {:.3f}'.format(lm.score(X_test, y_test)))\n",
    "print('Testing accuracy by using r2_score meric: {:.3f}'.format(r2_score(y_test, lm.predict(X_test))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "The accuracy on the test dataset is very similar to the training accuracy. From that we can conclude that the model has not overfitted the training data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.10 ==========\n",
    "Now we want to get a feel for how good the fit is, so we wil plot the measured values against the predicted ones. Make use of the function provided below which takes as input arguments the measured (`y_true`) and predicted (`y_pred`) values of a target variable and produces a scatter plot for the two by also including a straight line going through the origin. \n",
    "\n",
    "Where would you expect the circles to be for a perfect fit?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fit_scatter(y_true, y_pred):\n",
    "    assert y_true.shape == y_pred.shape\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.scatter(y_true, y_pred)\n",
    "    ax.plot([y_true.min(), y_true.max()], [y_true.min(), y_true.max()], 'k--', lw=4)\n",
    "    ax.set_xlabel('Measured')\n",
    "    ax.set_ylabel('Predicted')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "prediction = lm.predict(X_test)\n",
    "fit_scatter(y_test, prediction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "If the fit was perfect all the points would be sitting on the straight line going through the origin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.11 ==========\n",
    "Another way of assessing the performance of the model is to inspect the distribution of the errors. Make a histogram plot by using seaborn's `displot` function. This will also show an estimate of the underlying distribution.\n",
    "\n",
    "Does it look like the errors are normally distributed? Would you trust the fit of the distribution on the graph? Explain why."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "g = sns.distplot(y_test-prediction, rug=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "We should be very careful before making any judgements since the number of samples is relatively small in this case. However, the errors seem to follow a Gaussian (normal) distribution. There is some evidence to suggest the model is very slightly over predicting more than under predicting (given the skew)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.12 ==========\n",
    "Above we deleted the vendor variable. However, we can use nominal attributes in regression by converting them to numeric, exactly in the same way that we did at the first part of this lab. \n",
    "\n",
    "Now, use the original `cpu` dataset and convert the `vendor` attribute to numeric by using a `LabelEncoder`. Then train a linear regression model to the data and compare its performance to the one we had previously. Did adding the *binazired vendor* variable help? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "vendor_le = LabelEncoder()\n",
    "vendor_le.fit(cpu[\"vendor\"])\n",
    "cpu[\"vendor\"] = vendor_le.transform(cpu[\"vendor\"])\n",
    "X = cpu.drop('ERP', axis=1).as_matrix()\n",
    "y = cpu['ERP'].as_matrix()\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.9, test_size=0.1, random_state=0)\n",
    "lm = LinearRegression(fit_intercept=True, normalize=True, copy_X=True).fit(X_train,y_train)\n",
    "prediction_new = lm.predict(X_test)\n",
    "print('New accuracy on test set: {:.3f}'.format(lm.score(X_test,y_test)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "Adding the binarized vendor variable only helped a little bit. This was somehow expected, since we made additional use of a discrete variable to predict a continuous one."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}