08_Lab_4_Clustering_PCA_and_Evaluation_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory applied machine learning (INFR10069) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 4: Clustering, PCA, and Evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab we consider unsupervised learning in the form of clustering methods and principal component analysis (PCA), as well as more thorough performance evaluation of classifiers.\n",
    "\n",
    "All the datasets that you will need for this lab are located at the `./datasets` directory which is adjacent to this file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": 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\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.metrics import adjusted_rand_score\n",
    "from sklearn.decomposition import PCA\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Clustering the landsat dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first consider clustering of the Landsat data. For information about the Landsat data you can read [this description](http://www.inf.ed.ac.uk/teaching/courses/dme/html/landsat.html). Since there are 6 classes in the data, it would be interesting to try clustering with k=6 centres..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.1 ==========\n",
    "1. With pandas, read the csv dataset located at './datasets/landsat.csv'\n",
    "1. Split the data into the features `X` (pandas dataframe), and the labels `y` (easier to make it a numpy array)\n",
    "1. Give it a once over\n",
    "    * Get a feel for the size\n",
    "    * Check it looks ok\n",
    "    * Understand what the features are\n",
    "    * Plot the class distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "data_path = os.path.join(os.getcwd(), 'datasets', 'landsat.csv')\n",
    "landsat = pd.read_csv(data_path, delimiter = ',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "X = landsat.drop('class', axis=1)\n",
    "y = landsat['class'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "print('Number of instances: {}, number of attributes: {}'.format(X.shape[0], X.shape[1]))\n",
    "print(landsat.describe())\n",
    "# Notice that class is categorical (not numeric) and count\n",
    "labeldict = {1:'red soil', 2:'cotton crop', 3:'grey soil',\n",
    "             4:'damp grey soil', 5:'soil with vegetation stubble',\n",
    "             6:'mixture class (all types present)',\n",
    "             7:'very damp grey soil'}\n",
    "fig, ax = plt.subplots()\n",
    "landsat['class'].astype('category').value_counts().plot(kind='bar', ax=ax)\n",
    "labels = [int(ticklabel.get_text()) for ticklabel in ax.get_xticklabels()]\n",
    "ax.set_xticklabels([labeldict[l] for l in labels])\n",
    "plt.xlabel('Classes')\n",
    "plt.ylabel('Count')\n",
    "plt.title('Class distribution')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Advanced (not expecting you to have done this) - describe pixel values by band\n",
    "\n",
    "landsat_melt = pd.melt(landsat, id_vars=['class'], var_name='pixel_label')\n",
    "landsat_melt['band'] = landsat_melt['pixel_label'].str[-1]\n",
    "landsat_melt.groupby('band')['value'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.2 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Tip - don't get stuck on this, move on after 10 mins or so (it's not critical)*\n",
    "\n",
    "Plot a few datapoints. You'll need to understand and reshape the datapoints to do this. *Hint: try reading the [detailed description](http://www.inf.ed.ac.uk/teaching/courses/dme/html/satdoc.txt), it'll take you 2 minutes...`plt.imshow()` or `sns.heatmap()` (with an `ax.invert_yaxis()`) may also be handy ;)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "def get_images(row, bands=None):    \n",
    "    # Get each of the 4 3x3 images contained in the row\n",
    "    # Pixels are labeled 1 to 9 from topleft to bottom right\n",
    "    # They are measured at 4 spectral bands\n",
    "    #\n",
    "    # row - a numpy array\n",
    "    if bands is None:\n",
    "        img = [[]] * 4\n",
    "        for ii in range(4):\n",
    "            img[ii] = row[[4*p + ii for p in range(9)]].values.reshape((3,3)).astype(int)\n",
    "    elif isinstance(bands, int):\n",
    "        img = row[[4*p + bands for p in range(9)]].values.reshape((3,3)).astype(int)\n",
    "    elif isinstance(bands, list):\n",
    "        img = [[]] * len(bands)\n",
    "        for ii, band in enumerate(bands):\n",
    "            img[ii] = row[[4*p + band for p in range(9)]].values.reshape((3,3)).astype(int)\n",
    "    return img\n",
    "        \n",
    "for ii in range(4):\n",
    "    fig, ax = plt.subplots(1,4)\n",
    "    plt.suptitle('Row {}'.format(ii), fontsize=16)\n",
    "    for jj, img in enumerate(get_images(landsat.iloc[ii,:])):\n",
    "        ax[jj] = sns.heatmap(img, annot=True, fmt=\"d\", ax=ax[jj], vmin=0, vmax=255, \n",
    "                    cbar=False, square=True, cmap=plt.cm.gray)\n",
    "        ax[jj].invert_yaxis()\n",
    "        ax[jj].set_title('Band {}'.format(jj))\n",
    "    plt.tight_layout()\n",
    "    plt.subplots_adjust(top=1.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### ========== Question 1.3 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Read [this demonstration of k-means clustering assumptions](http://scikit-learn.org/0.17/auto_examples/cluster/plot_kmeans_assumptions.html#sphx-glr-auto-examples-cluster-plot-kmeans-assumptions-py) from the sklearn documentation. Get a feel for how to create and fit a k-means object and what the different arguments do.\n",
    "\n",
    "Initialise a [k-means clustering](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html#sklearn.cluster.KMeans) object with 6 clusters, and one other parameter that ensures you can reproduce your results (other arguments kept as default). Call the object `kmeans`. Use the `fit()` method to fit to the training data (`X` - the features of `landsat` created above)\n",
    "\n",
    "**Be careful to fit `X` - only the features - not the class labels!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "# It's really important that you set a random_state such that you can\n",
    "# reproduce your results\n",
    "kmeans = KMeans(n_clusters=6, random_state=1337)  \n",
    "kmeans.fit(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.4 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, how well did that work? Are the classes well separated and form 6 nice clusters? The sklearn documentation gives a great introduction to k-means [here](http://scikit-learn.org/stable/modules/clustering.html#k-means). It describes what the algorithm is trying to minimise - the squared difference between datapoints and their closest cluster centre - a.k.a. the **inertia**. Lower inertia implies a better fit.\n",
    "\n",
    "Since we have the true class labels in this case, we can use another metric: the [adjusted rand index](http://scikit-learn.org/stable/modules/clustering.html#clustering-performance-evaluation). Understand what it is and roughly how it is calculated (try the [mathematical formulation on sklearn](http://scikit-learn.org/stable/modules/clustering.html#mathematical-formulation) or [on wikipedia](https://en.wikipedia.org/wiki/Rand_index)).\n",
    "\n",
    "Print the `inertia` and the `adjusted_rand_score` of the kmeans object. The inertia is contained within the `kmeans` object you just fitted as a property. You need to use `y` and the cluster labels (another property of the `kmeans` object you just made), and are welcome to use the sklearn metrics function [adjusted_rand_score](http://scikit-learn.org/0.17/modules/generated/sklearn.metrics.adjusted_rand_score.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "kmeans.inertia_, adjusted_rand_score(y, kmeans.labels_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 1.5 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the counts of the labels within each cluster. If the clustering has worked well, and the labels are inticative of genuine difference in the data, we should expect each cluster to have one dominant label.\n",
    "\n",
    "Use `sns.countplot` on `kmeans.labels_` with a hue of `y` to get a plot that counts the number of instances within each cluster, and breaks them down by the class labels.\n",
    "\n",
    "Below the plot, comment on:\n",
    "1. How successful the clustering has been at separating data with different labels\n",
    "1. Focussing on clusters, which are the best?\n",
    "1. Focussing on labels, which are well identified by the clustering?\n",
    "1. Which labels are the hardest to determine by the data point cluster assignment?\n",
    "\n",
    "**Extension**: create a vector `y_labels` from `y` which has the names for the classes for a more interpretable plot with respect to the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax = sns.countplot(x=kmeans.labels_, hue=y, ax=ax)\n",
    "plt.xlabel('cluster id')\n",
    "plt.show()\n",
    "\n",
    "# --- Extension ---\n",
    "cats = ['red soil','cotton crop','grey soil', 'damp grey soil', \n",
    "        'soil with vegetation stubble', 'mixture class (all types present)',\n",
    "        'very damp grey soil']\n",
    "labeldict = dict(zip(range(1,8), cats))\n",
    "y_labels = pd.Categorical([labeldict[y_i] for y_i in y],\n",
    "                          categories=[c for c in cats if c != 'mixture class (all types present)']) # ensures colours same\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax = sns.countplot(x=kmeans.labels_, hue=y_labels, ax=ax)\n",
    "ax.legend(loc='center left', bbox_to_anchor=[1.1, 0.5], title='Class')\n",
    "plt.xlabel('cluster id')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "1. How successful the clustering has been at separating data with different labels\n",
    "    * Hmmm...so so! There are 3 clusters that contain one main type of datapoint: label 3 (grey soil), label 2 (cotton crop), and label 1 (red soil). The other 3 clusters contain a mix.\n",
    "1. Focussing on clusters, which are the best?\n",
    "    * If a datapoint in the training set is in the majority label 2 cluster (cotton crop), you can be certain it has label 2, the majority label 1 (red soil), and majority label 3 (grey soil) clusters are similar but not quite as pure\n",
    "1. Focussing on labels, which are well identified by the clustering?\n",
    "    * Label 3 (grey soil) data seems well separated in general - most of the datapoints with this label are in one cluster. It's important to note that, though the majority label 2 (cotton crop) cluster contains only label 2 datapoints, there are label 2 data in other clusters. From this we can say that it appears label 2 is mixed up with label 1 (red soil) and label 5 (soil with vegetation stubble).\n",
    "1. Which labels are the hardest to determine by the data point cluster assignment?\n",
    "    * label 4 (damp grey soil) is confused with about 1/3 of the label 7 (very damp grey soil) data...this doesn't seem particularly surpising. This label 7 data (very damp grey soil) is often confused with label 5 data (soil with vegetation stubble)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Dimensionality reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The landsat data is 36 dimensional, so we cannot visualise it, with respect to class, on a nice two dimensional plot. Additionally, as dimensionality increases, euclidean distance [becomes less meaningful](https://en.wikipedia.org/wiki/Curse_of_dimensionality#Distance_functions)...\n",
    "\n",
    "Perhaps if we found a lower dimensional subspace the data lies upon, we could more easily distinguish the datapoints..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.1 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have a look at the [PCA example](http://scikit-learn.org/0.17/auto_examples/decomposition/plot_pca_iris.html) in the sklearn documentation. For more information about PCA and decomposition in general check out the sklearn [user guide on decomposition](http://scikit-learn.org/stable/modules/decomposition.html#pca).\n",
    "\n",
    "We are going to project the data down to 2 dimensions and visualise it using PCA. \n",
    "\n",
    "1. Create `pca`, an instance of an [sklearn PCA object](http://scikit-learn.org/0.17/modules/generated/sklearn.decomposition.PCA.html#sklearn.decomposition.PCA), setting n_components to 2. \n",
    "1. Create `X_2d` by using the pca method `fit_transform()` and supplying the features `X` to fit and transform to 2d."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "pca = PCA(n_components=2)\n",
    "X_2d = pca.fit_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.2 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualise the data! Use a scatterplot and colour the datapoints by their class. You'll find [this example](http://scikit-learn.org/stable/auto_examples/decomposition/plot_pca_vs_lda.html) very helpful to stea-I mean, adapt. \n",
    "\n",
    "Below the plot, comment on whether the data looks more or less seperable now. Do the data look like they will be confused by a k-means clustering in the same way now?\n",
    "\n",
    "**Extension**: Unless the data is somehow magically perfectly seperable, you may want to try and describe the space a little better than a scatterplot (bacause points are plotted on top of one another). Try and make a plot that clarifies the location of the classes. We're actually interested in where the **density** is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "colours = matplotlib.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "lw = 2\n",
    "\n",
    "sub_cats = [c for c in cats if c != 'mixture class (all types present)']\n",
    "for color, target_name in zip(colours, sub_cats):\n",
    "    plt.scatter(X_2d[y_labels == target_name, 0], X_2d[y_labels == target_name, 1], color=color, alpha=.5, lw=lw,\n",
    "                label=target_name)\n",
    "plt.axis('equal')\n",
    "plt.legend(loc='center left', scatterpoints=3, bbox_to_anchor=[1.01, 0.5])\n",
    "plt.title('Labelled data in PCA space')\n",
    "plt.xlabel('PC1')\n",
    "plt.ylabel('PC2')\n",
    "top_plot = plt.gca()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# --- extension ---\n",
    "import matplotlib.lines as mlines\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "colours = matplotlib.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "lines = []\n",
    "for ii, cat in enumerate(sub_cats):\n",
    "    idx = y_labels==cat\n",
    "    g = sns.kdeplot(X_2d[idx, 0], X_2d[idx, 1], ax=ax, n_levels=10, \n",
    "                    cmap=sns.light_palette(colours[ii], as_cmap=True))\n",
    "ax.set_xlim(top_plot.get_xlim())\n",
    "ax.set_ylim(top_plot.get_ylim())\n",
    "patches = [mlines.Line2D([], [], color=colours[ii], label=cat) for ii, cat in enumerate(sub_cats)]\n",
    "plt.legend(patches, sub_cats, loc='center left', scatterpoints=3, bbox_to_anchor=[1.01, 0.5])\n",
    "plt.title('Kernel density plots for each class')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "Without performing a clustering, it's hard to say whether the data is more or less seperable. However, we can see 6 definite clusters with some overlap in the top plot, but it's not all that clear. The kernel density plot allows us to see better where the density lies. These clusters are not spherical and there is overlap, so k-means will always fail. See in particular the 'red soil' data points vs. the 'soil with vegetation stubble points. As expected, we see the 'cotton crop' data points are very distinct from the rest. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.3 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find out whether it's easier to model the transformed data. Fit k-means to the transformed data and report the inertia and the adjusted rand index. Below, comment on whether it is better or worse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "kmeans_pca = KMeans(n_clusters=6, random_state=1337)  \n",
    "kmeans_pca.fit(X_2d)\n",
    "kmeans_pca.inertia_, adjusted_rand_score(y, kmeans_pca.labels_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer goes here*\n",
    "\n",
    "The fit is not as good as before, the adjusted rand index is lower. Note that we cannot compare the inertia with before. You should think about why."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ========== Question 2.4 =========="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The blobs in the 2 dimensional plot do look quite gaussian...try another classifier on the 2d data and see if it can perform better. What about using 3 principal component? \n",
    "\n",
    "Maybe there are subclasses within each class? Maybe increasing the number of clusters will increase your `adjusted_rand_score`.\n",
    "\n",
    "Use the adjusted rand score for fair comparison. Why do you think it works better or worse? Discuss with your colleagues and lab tutors why you think you got better/worse/the same results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here\n",
    "\n",
    "# Here are just some examples you could try\n",
    "# Anything goes for this question\n",
    "\n",
    "\n",
    "# In my solution I find that there seems to be a sweet spot around 10 clusters\n",
    "# This provides some evidence that the classes are multimodal. You can see that\n",
    "# in the plots.\n",
    "#\n",
    "# increasing to 3 principal component dimensions\n",
    "# provides a big benefit for kmeans and a Gaussian Naive Bayes classifier.\n",
    "# In fact, it's better than using the full data. \n",
    "# Using a random forest, I get the best results using the full data - N.B.\n",
    "# I'm using the Out-of-bag error here (the training error is 0% as discussed before)\n",
    "# so we may find that this is overfitting if given a test set.\n",
    "\n",
    "\n",
    "\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "\n",
    "for ii in range(1,20,2):\n",
    "    km = KMeans(n_clusters=ii, random_state=1337)  \n",
    "    km.fit(X_2d)\n",
    "    y_pred = km.labels_\n",
    "    print('\\nk-means with {} clusters results\\nARI: {}, Inertia: {}'.\\\n",
    "          format(ii, adjusted_rand_score(y, y_pred), km.inertia_))\n",
    "\n",
    "\n",
    "for ii in range(2, 6):\n",
    "    X_pca = PCA(n_components = ii).fit_transform(X)\n",
    "    km = KMeans(n_clusters=6, random_state=1337)  \n",
    "    km.fit(X_pca)\n",
    "    y_pred = km.labels_\n",
    "    print('\\n{}D PCA data, kmeans results\\nARI: {}, Inertia: {}'.\\\n",
    "          format(ii, adjusted_rand_score(y, y_pred), km.inertia_))\n",
    "\n",
    "for ii in range(2, 6):\n",
    "    X_pca = PCA(n_components = ii).fit_transform(X)\n",
    "    gnb = GaussianNB()\n",
    "    gnb.fit(X_pca, y)\n",
    "    y_pred = gnb.predict(X_pca)\n",
    "    print('\\n{}D PCA data, GNB results\\nARI: {}, Mean Accuracy: {}'.\\\n",
    "          format(ii, adjusted_rand_score(y, y_pred), gnb.score(X_pca, y)))\n",
    "\n",
    "gnb = GaussianNB()\n",
    "gnb.fit(X, y)\n",
    "y_pred = gnb.predict(X)\n",
    "print('\\nFull data, GNB results\\nARI: {}, Mean Accuracy: {}'.format(adjusted_rand_score(y, y_pred), gnb.score(X, y)))\n",
    "\n",
    "for ii in range(2, 6):\n",
    "    X_pca = PCA(n_components = ii).fit_transform(X)\n",
    "    rf = RandomForestClassifier(n_estimators=100, random_state=1337, oob_score=True)\n",
    "    rf.fit(X_pca, y)\n",
    "    y_pred = np.argmax(rf.oob_decision_function_, axis=1)\n",
    "    print('\\n{}D PCA data, RF results\\nARI: {}, Mean Accuracy: {}'.\\\n",
    "          format(ii, adjusted_rand_score(y, y_pred), rf.oob_score_))\n",
    "\n",
    "rf = RandomForestClassifier(n_estimators=100, random_state=1337, oob_score=True)\n",
    "rf.fit(X, y)\n",
    "y_pred = np.argmax(rf.oob_decision_function_, axis=1)\n",
    "print('\\nFull data, RF results\\nARI: {}, Mean Accuracy: {}'.\\\n",
    "      format(adjusted_rand_score(y, y_pred), rf.oob_score_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "fig = plt.figure(1, figsize=(8, 6))\n",
    "ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=30, azim=300)\n",
    "\n",
    "pca = PCA(n_components=3)\n",
    "pca.fit(X)\n",
    "X_3d = pca.transform(X)\n",
    "for ii, cat in enumerate(sub_cats):\n",
    "    idx = y_labels==cat\n",
    "    ax.scatter(X_3d[idx, 0], X_3d[idx, 1], X_3d[idx, 2], c=colours[ii],\n",
    "               label = cat)\n",
    "plt.legend(loc='center left', bbox_to_anchor=[1.1, .5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Credits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lab prepared by Lawrence Murray and Chris Williams, November 2008; revised Athina Spiliopoulou Nov 2009; revised Sean Moran Nov 2011; revised Boris Mitrovic Oct 2013; revised and converted python by James Owers and Agamemnon Krasoulis Oct 2016"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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
}