machine_learning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Machine Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook provides an introduction to machine learning and walks through how to develop workflows for training machine learning models.\n",
    "\n",
    "This lesson is prepared by:\n",
    "- [Kevin Chai](http://computation.curtin.edu.au/about/computational-specialists/health-sciences/)\n",
    "- [Rebecca Lange](http://computation.curtin.edu.au/about/computational-specialists/humanities/)\n",
    "\n",
    "from the [Curtin Institute for Computation](http://computation.curtin.edu.au) at Curtin University in Perth, Australia. \n",
    "\n",
    "Some of the materials in this notebook have been referenced and adapted from:\n",
    "- [Randal Olsen's Data Science Notebook](https://github.com/rhiever/Data-Analysis-and-Machine-Learning-Projects/tree/master/example-data-science-notebook)\n",
    "- [Sebastian Raschka's Python Machine Learning Notebooks](https://github.com/rasbt/python-machine-learning-book)\n",
    "- [Kevin Markham's Scikit Learn Notebooks](https://github.com/justmarkham/scikit-learn-videos)\n",
    "\n",
    "Make sure to open this notebook in the root directory of the code repository.\n",
    "\n",
    "This work is made available under the [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. [Introduction](#1.-Introduction)\n",
    "2. [Problem definition](#2.-Problem definition)\n",
    "3. [Data cleaning](#3.-Data-cleaning)\n",
    "4. [Exploratory analysis](#4.-Exploratory-analysis)\n",
    "5. [Data preparation](#5.-Data-preparation)\n",
    "6. [Classification](#6.-Classification)\n",
    "7. [Regression](#7.-Regression)\n",
    "8. [Clustering](#8.-Clustering)\n",
    "9. [Dimensionality reduction](#9.-Dimensionality-reduction)\n",
    "10. [Reproducibility](#10.-Reproducibility)\n",
    "11. [Conclusions](#11.-Conclusions)\n",
    "12. [Further reading](#12.-Further-reading)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Required libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "This notebook uses several Python packages that come standard with the [Anaconda Python distribution](http://continuum.io/downloads). The primary libraries that we'll be using are:\n",
    "\n",
    "* **NumPy**: a fast numerical array structure and helper functions.\n",
    "* **pandas**: a DataFrame structure to store data in memory and work with it easily and efficiently.\n",
    "* **scikit-learn**: a machine learning package.\n",
    "* **matplotlib**: a basic plotting library; most other plotting libraries are built on top of it.\n",
    "* **seaborn**: a advanced statistical plotting library.\n",
    "\n",
    "To make sure you have all of the packages you need, install them with `conda`:\n",
    "\n",
    "    conda install numpy pandas scikit-learn matplotlib seaborn\n",
    "\n",
    "`conda` may ask you to update some of the packages if you don't have the most recent version. Allow it to do so.\n",
    "\n",
    "Alternatively, if you can install the packages with [pip](https://pip.pypa.io/en/stable/installing/) (a Python package manager):\n",
    "\n",
    "    pip install numpy pandas scikit-learn matplotlib seaborn\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Machine learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Machine learning refers to algorithms that learn from examples and experiences (data). These algorithms often perform better than explicitly hard coded rules for complex tasks. i.e. where it is difficult to describe precise rules.\n",
    "\n",
    "We will explore two types of machine learning algorithms in this notebook:\n",
    "\n",
    "**1. Supervised Learning** \n",
    "* make predictions using data\n",
    "* There is an outcome we are trying to predict\n",
    "* Example: Is an e-mail spam or ham?\n",
    "\n",
    "![Spam filter](images/spam_filter.png)\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Classification example</div>\n",
    "\n",
    "The workflow for developing and deploying a supervised machine learning model is shown below.\n",
    "\n",
    "![Supervised Learning Workflow](images/supervised_learning_workflow.jpg)\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Supervised learning workflow</div>\n",
    "\n",
    "**2. Unsupervised Learning**\n",
    "* Extract structure from data\n",
    "* There is no \"right answer\"\n",
    "* Example: Segment grocery store shoppers into clusters that exhibit similar behaviours\n",
    "\n",
    "![Clustering](images/clustering.png)\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Clustering example</div>\n",
    "\n",
    "The unsupervised learning workflow is illustrated in  the figure below. Observe that no labels are used for training. \n",
    "\n",
    "![Unsupervised Learning Workflow](images/unsupervised_learning_workflow.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Problem definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "We have been tasked to develop a machine learning model to classify galaxies by their morphology (appearance) from a dataset containing measurements such as:\n",
    "\n",
    "- redshift\n",
    "- magnitude in _ugriz_ bands\n",
    "- exponential scale radius and ellipticity\n",
    "- de Vaucouleurs scale radius and ellipticity\n",
    "- stellar mass\n",
    "- ...\n",
    "- etc.\n",
    "\n",
    "Galaxies in the dataset have labels taken from [Galaxy Zoo](https://www.galaxyzoo.org/) DR1 - Table 2. Galaxy Zoo is described in Lintott et al. 2008, MNRAS, 389, 1179 and the data release is described in Lintott et al. 2011, 410, 166. \n",
    "\n",
    "We use the final debiased labels to categorise a galaxy as:\n",
    "\n",
    "- spiral\n",
    "- elliptical\n",
    "\n",
    "![Galaxies](images/spiral_ellipse_galaxies.jpg)\n",
    "\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: A spiral galaxy (left) and elliptical galaxy (right)</div>\n",
    "\n",
    "Our goal is to train a model that can accurately classify galaxies. We want our model to generalise well. i.e. it can correctly classify unseen galaxies (i.e. galaxies not in our training dataset)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "The first step of developing a machine learning model is to clearly understand and define the problem we want to solve. Here are some useful questions to ask. \n",
    "\n",
    "#### Does the dataset contain labels?\n",
    "\n",
    "Yes, the dataset contains galaxy morphology labels (spiral and elliptical). Therefore, we can use supervised learning methods to utilise the labelled data.\n",
    "\n",
    "#### What is the type of problem?\n",
    "\n",
    "We want to classify galaxies by their morphology so this is a classification problem.\n",
    "\n",
    "#### What metric can be used to evaluate the model?\n",
    "\n",
    "Since we are performing classification, we can use a classification metric such as [accuracy](https://en.wikipedia.org/wiki/Accuracy_and_precision) to evaluate and quantify the performance of our model. Accuracy is the percentage of correctly classified galaxies. As a challenge, we have been asked to develop a model that achieves ≥ 70% accuracy. Is this feasible? How can we find out?\n",
    "\n",
    "#### Can the question actually be answered with the available data?\n",
    "\n",
    "Yes, the dataset contains labels of the galaxy types we want to classify. i.e. `{spiral, elliptical}`. If we were asked to develop a model to detect other types of galaxies then we would need to collect data and labels for these other types.\n",
    "\n",
    "#### Are there problem and application requirements we need to consider?\n",
    "\n",
    "How will the model be used after it has been developed? Are there requirements for the model to be used on specific devices? e.g. personal computers, smart phones, embedded systems and/or servers. For the purpose of this lesson, we can assume that it is fine to train and run the model on our personal computers. However, this may not always be the case and it is something to consider if we are required to deploy our models to systems with limited compute and memory resources (e.g. embedded systems). In these instances, we may be restricted to build small, less complex and fast models at the cost of reduced performance (e.g. accuracy).\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Thinking about and documenting the problem we're working on is an important step to performing effective machine learning (and data analysis in general) that sometimes gets overlooked.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## 3. Data cleaning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "#### Dataset\n",
    "\n",
    "For this lesson, we will be using data from [Galaxy Zoo DR1](https://www.galaxyzoo.org/) and the Sloan Digital Sky Survey (SDSS) ([using the DR9 SQL search](http://skyserver.sdss.org/dr9/en/tools/search/sql.asp)). \n",
    "\n",
    "The data dictionary for this dataset is presented in Table 1. This dataset is limited to first 5,000 Galaxy Zoo classified galaxies which have spectra in the SDSS database. The debiased fraction of the votes in elliptical and spiral categories is given, along with columns identifying systems classified as spiral, elliptical or uncertain.\n",
    "\n",
    "<p style=\"text-align:center;font-weight:bold\">Table 1: Data dictionary</p>\n",
    "\n",
    "| Column           | Description                                                            |\n",
    "|:-----------------|:-----------------------------------------------------------------------|\n",
    "| id               | Unique SDSS ID composed of [skyVersion, rerun, run, camcol, field, obj]|\n",
    "| ra               | right Ascension  (HMS)                                                 |\n",
    "| dec              | declination (DMS)                                                      |\n",
    "| redshift         | redshift                                                               |\n",
    "| mag_u            | magnitude _u_ band                                                     |\n",
    "| mag_g            | magnitude _g_ band                                                     |\n",
    "| mag_r            | magnitude _r_ band                                                     |\n",
    "| mag_i            | magnitude _i_ band                                                     |\n",
    "| mag_z            | magnitude _z_ band                                                     |\n",
    "| deVRad_r         | de Vaucouleurs scale radius fit in _r_ band                            |\n",
    "| deVAB_r          | ellipticity from de Vaucouleurs fit in _r_ band                        |\n",
    "| expRad_r         | exponential scale radius fit in _r_ band                               |\n",
    "| expAB_r          | ellipticity from exponential fit in _r_ band                           |\n",
    "| stellar_mass     | log galaxy mass (in units of solar mass)                               |\n",
    "| votes            | number of Galaxy Zoo annotators                                        |\n",
    "| p_el_debiased    | debiased labelling probability the galaxy is elliptical                |\n",
    "| p_cs_debiased    | debiased labelling probability the galaxy is spiral                    |\n",
    "| spiral           | label = spiral galaxy {0=False, 1=True}                                |\n",
    "| elliptical       | label = elliptical galaxy {0=False, 1=True}                            |\n",
    "| uncertain        | label = uncertain {0=False, 1=True}                                    |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Understand and clean the data\n",
    "\n",
    "The next step is to explore and familiarise ourselves with the data. Datasets can have errors, and it's important to identify and resolve these errors before rushing in to develop our models. Errors in the dataset that aren't resolved will propogate through the machine learning pipeline.\n",
    "\n",
    "Generally, we're looking to answer the following questions:\n",
    "\n",
    "* Is there anything wrong with the data?\n",
    "* Are there any quirks with the data?\n",
    "* Do I need to clean or remove some of the data?\n",
    "\n",
    "Let's start by reading the data into a `pandas` DataFrame and inspecting the first few rows.\n",
    "\n",
    "If you want to learn more about pandas you can find an introductory notebook from a session previously run by [ADACS](https://github.com/ADACS-Australia/Face-to-Face)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Import pandas\n",
    "import pandas as pd\n",
    "\n",
    "# Read the data file\n",
    "data = pd.read_csv('data/galaxies.csv')\n",
    "\n",
    "# Display the first few rows\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "One of the first things we should do is to look for is missing data. We can use `pandas` to inspect rows with missing data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Get records containing any missing values\n",
    "data[data.isnull().any(axis=1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Here we see a number of rows that contain missing values (`NaN`). There are many [methods for handling and imputing missing data](http://machinelearningmastery.com/handle-missing-data-python/). However, from our results we see that the galaxies with missing values have been labelled as `uncertain` so we will simply remove these records. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Drop the rows containing any missing values\n",
    "data.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's double check to make sure the records have been removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Should return 0 rows\n",
    "data[data.isnull().any(axis=1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "It is important to construct a good training dataset for developing a robust model. In machine learning, we also refer to the training dataset as the ground truth. Any uncertainities or errors in the ground truth can confuse and reduce the performance of our models.\n",
    "\n",
    "In order to build a good training dataset, we only want to keep records where we are confident there is agreement between the Galaxy Zoo annotators in labelling galaxies. From the dataset, it appears we can achieved this by:\n",
    "\n",
    "- Keeping records where the uncertain label column is equal to False (0)\n",
    "\n",
    "Let's try this method and perform some sanity checking to see if the resulting data subset has debiased probabilities for ellipital and spiral labels ≥ 0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# 1. Filter records using the uncertain column \n",
    "f1 = data[data['uncertain'] == 0]\n",
    "print('# records after removing uncertain rows: %d' % (len(f1)))\n",
    "\n",
    "# 2. Filter records using the debiased probabilities ≥ 0.5\n",
    "f2 = f1[(f1['p_el_debiased'] >= 0.5) | (f1['p_cs_debiased'] >= 0.5)]\n",
    "print('# records with ≥ 0.5 debiased probability: %d' % (len(f2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "There appears to be a discrepency of 5 records between the two data subsets. Let's investigate further by inspecting these records."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Display filter #1 rows that don't appear in the filter #2 dataset\n",
    "f1.loc[f1.index.difference(f2.index)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "These records have been labelled as spiral galaxies but the debiased probabilities (confidence) for these labels are below the 0.5 confidence threshold. i.e. `p_cs_debiased < 0.5`. Additionally, for the last 2 records, we see that debiased probabilities for labelling the galaxy as elliptical is actually larger than the spiral probabilities which suggests these might be due to labelling errors. \n",
    "\n",
    "It appears we can't rely on filtering the dataset on the `uncertain` column alone so we will apply both filters to clean the dataset. Sanity checking your dataset is a step that is regularly overlooked but dedicating time to explore and resolve these errors can (sometimes dramatically) improve the performance of our machine learning models.\n",
    "\n",
    "We will make a copy of the filtered dataset in the variable `df` (abbreviation for DataFrame) for our clean dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "df = f2.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Let's clean and prepare our dataset by:\n",
    "\n",
    "- creating a class column to represent the galaxy type label `{0=elliptical, 1=spiral}`\n",
    "- removing unwanted columns that won't be used in our analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Create the class column - do this by renaming the 'spiral' column to 'class'\n",
    "df.rename(columns={'spiral': 'class'}, inplace=True)\n",
    "\n",
    "# Remove unwanted columns\n",
    "unwanted_columns = ['ra', 'dec', 'votes', 'p_el_debiased', \n",
    "                    'p_cs_debiased', 'elliptical', 'uncertain']\n",
    "df.drop(unwanted_columns, axis=1, inplace=True)\n",
    "\n",
    "# Inspect the first few rows of our prepared dataset\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Next, it's always a good idea to look to have a closer look at our data — especially for outliers. Let's start by printing out some summary statistics about the data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Display dataset summary statistics\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The summary table provides some useful information:\n",
    "\n",
    "- there are 1,801 galaxies in our dataset\n",
    "- the mean of the `class` column indicates that ~75.3% of the dataset are spiral galaxies\n",
    "\n",
    "There also appears to be negative values for the magnitude (`mag`) columns, radius and ellipticity columns (`deVRad` and `expRad`) and the `stellar_mass` which look like errors.\n",
    "\n",
    "Let's examine all records that have ≤ 0 values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Get records with ≤ 0 values \n",
    "# Exclude the first (id) and last (class) column in the search\n",
    "df[(df.iloc[:, 1:-1] <= 0).any(axis=1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Using our domain knowledge, we know that -9999 is often used to represent a missing value in (optical) astronomy datasets. It also doesn't make sense to have galaxies with a zero or negative magnitude and stellar mass values. We can treat these as missing values / errors and remove them from our dataset. Let's display the summary statistics on the reduced dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Only retain records with ≥ 0 values in our dataset\n",
    "df = df[(df.iloc[:, 1:-1] > 0).all(axis=1)]\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "The summary of the reduced dataset looks better:\n",
    "\n",
    "- min values for the `mag` columns range from 11.729 to 14.966\n",
    "- min values for the `deVRad` and `expRad` columns range from 0.078 to 0.388\n",
    "- min value for the `stellar_mass` = 8.096\n",
    "\n",
    "Our cleaned dataset contains 1,797 galaxies and 75.3% of the records are spiral galaxies.\n",
    "\n",
    "Tables like this are useful when we know that our data should fall in a particular range. e.g. It did not make sense to have zero or negative values for some of our columns. However, it is usually better to visualize the data in some way. Visualization makes outliers and errors immediately stand out, whereas they might go unnoticed in a large table of numbers.\n",
    "\n",
    "Since we know we're going to be plotting in this section, let's set up the notebook so we can plot inside of it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Command to show plots inside of the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Next, let's create a **scatterplot matrix**. Scatterplot matrices plot the distribution of each column along the diagonal and a scatterplot matrix for the combination of each variable. They make for an efficient tool to look for errors in our data.\n",
    "\n",
    "Since we have quite a few columns, let's generate two scatterplots for data containing:\n",
    "\n",
    "1. measurements\n",
    "2. magnitudes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Measurements scatterplot\n",
    "measurements = ['redshift', 'stellar_mass', 'deVRad_r', 'deVAB_r', 'expRad_r', 'expAB_r']\n",
    "selected_columns = measurements + ['class']\n",
    "sb.pairplot(df[selected_columns], hue='class', plot_kws={'alpha': 0.3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "From the scatterplot matrix, we see potential outliers: \n",
    "\n",
    "- the `stellar_mass` value of a `spiral` galaxy falls outside its normal range (a value of ~8).\n",
    "- two spiral galaxies with `expRad_r` ≥ 17\n",
    "\n",
    "Let's inspect the outlier records to determine whether if we should keep it in our training dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Inspect possible outliers\n",
    "df[(df['stellar_mass'] < 8.1) | (df['expRad_r'] >= 17)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "There is nothing obvious that leads us to believe the records are erroneous. Let's keep them in our dataset.\n",
    "\n",
    "Fixing outliers can be difficult. It's often unclear whether the outlier was caused by measurement error, recording the data in improper units, or if the outlier is a real anomaly. For that reason, we should be careful when working with outliers. If we decide to exclude any data, we need to make sure to document what data was excluded and provide reasons for excluding that data. i.e. The data didn't fit my hypothesis will not stand peer review.\n",
    "\n",
    "Let's plot the magnitudes scatterplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Magnitudes scatterplot\n",
    "magnitudes = ['mag_u', 'mag_g', 'mag_r', 'mag_i', 'mag_z']\n",
    "selected_columns = magnitudes + ['class']\n",
    "sb.pairplot(df[selected_columns], hue='class', plot_kws={'alpha': 0.3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "From the magnitude scatterplot matrix, we see a potential outlier: the `mag_u` value of a `elliptical` galaxy falls outside its normal range (a value of ~24)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Inspect possible outliers\n",
    "df[df['mag_u'] > 24]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "There is nothing obvious that leads us to believe the record is erroneous. Let's keep it in our dataset just like before.\n",
    "\n",
    "After all this hard work, we don't want to repeat the data cleaning process every time we work with the data set. Let's save the cleaned dataset as a separate file and work directly with that data file from now on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# Save the clean dataset to a file \n",
    "df.to_csv('data/galaxies-clean.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Exploratory analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Exploratory analysis is the step where we start delving deeper into the data set beyond the outliers and errors. We'll be looking to answer questions such as:\n",
    "\n",
    "- How is my data distributed?\n",
    "- Are there any correlations in my data?\n",
    "- Are there any confounding factors that explain these correlations?\n",
    "\n",
    "This is the stage where we plot the data in as many ways as possible. Create many charts, but don't bother making them pretty as these charts are for internal use.\n",
    "\n",
    "For completeness, let's generate a scatterplot for all features in one figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot scatterplot matrix with all features\n",
    "selected_columns = measurements + magnitudes + ['class']\n",
    "sb.pairplot(df[selected_columns], hue='class', plot_kws={'alpha': 0.3})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a number of observations from this plot:\n",
    "\n",
    "1. `reshift`, `stellar_mass`, `mag_u` and `mag_g` follow a normal distribution. This is something to keep in mind if we use modelling methods that assume the data is normally distributed.\n",
    "\n",
    "2. There is a strong positive correlation between the `expAB_r` and `deVAB_r` which indicates we might not need to use both of these features in our model.\n",
    "\n",
    "3. There is a strong positive correlation between the `mag_r`, `mag_i` and `mag_z`.\n",
    "\n",
    "4. There is a positive correlation between `mag_u` and `mag_g` with the other magnitude bands but not as strong as in observation #3. This suggests these features might provide added information not present in the `mag_r`, `mag_i` and `mag_z` bands.\n",
    "\n",
    "5. There are no clear pair-wise combinations of features that can easily seperate the two classes. \n",
    "\n",
    "Distnguishing between `spiral` and `elliptical` galaxies might be difficult given how much these features interrelate.\n",
    "\n",
    "We can also make **violin plots** of the data to compare the measurement distributions of the classes. Violin plots contain the same information as [box plots](https://en.wikipedia.org/wiki/Box_plot), but also scales the box according to the density of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# General violin plots of all features\n",
    "selected_columns = measurements + magnitudes + ['class']\n",
    "plt.figure(figsize=(13, 13))\n",
    "for column_index, column in enumerate(selected_columns):\n",
    "    if column == 'class':\n",
    "        continue\n",
    "    plt.subplot(3, 4, column_index + 1)\n",
    "    sb.violinplot(x='class', y=column, data=df[selected_columns])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are numerous errors and possible scenarios you may face while exploring and cleaning your data.\n",
    "\n",
    "The general takeaways here should be:\n",
    "\n",
    "- Handle missing data: replace it if you can or drop it\n",
    "- Ensure your data is encoded properly\n",
    "- Check if your data falls within the expected range and use domain knowledge whenever possible to define that expected range\n",
    "- Avoid tidying your data manually because that is not easily reproducible\n",
    "- Use code as a record of how you tidied your data\n",
    "- Plot everything you can about the data at this stage of the analysis so you can visually confirm everything looks correct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Cleaning and exploring the data is a important component to any machine learning project. If we had jumped straight in to modelling, we would have trained a model with errors in the dataset. Bad data leads to bad models.\n",
    "\n",
    "Now it's a good time to introduce some machine learning terminology."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Terminology\n",
    "\n",
    "- Each row is an `example` (also known as (aka): observation, sample, instance, record)\n",
    "- Each column is a `feature` (aka: predictor, attribute, dimension, independent variable)\n",
    "- The value we are predicting is the `label` (aka: target, outcome, dependent variable) \n",
    "\n",
    "In classification, the label is also referred to as the class or category."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Notation\n",
    "\n",
    "Let:\n",
    "- `m` = the number of examples in our dataset\n",
    "- `n` = the number of features in our dataset\n",
    "\n",
    "The dataset can be represented by two variables as shown in Table 2:\n",
    "\n",
    "1. A matrix `X` containing the examples and features of size `m x n`\n",
    "2. A vector `y` containing the labels of size `m`\n",
    "\n",
    "<p style=\"text-align:center;font-weight:bold\">Table 2: Dataset notations</p>\n",
    "\n",
    "<table>\n",
    "    <thead>\n",
    "        <tr>\n",
    "            <th style=\"text-align:center\" colspan=3>X</th>\n",
    "            <th style=\"text-align:center\">y</th>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <th style=\"text-align:center\">feature 1 ($x_1$)</th>\n",
    "            <th style=\"text-align:center\">feature 2 ($x_2$)</th>\n",
    "            <th style=\"text-align:center\">feature 3 ($x_3$)</th>\n",
    "            <th style=\"text-align:center\">label</th>\n",
    "        </tr>\n",
    "    </thead>\n",
    "    <tbody>\n",
    "        <tr>\n",
    "            <td style=\"text-align:center\">5</td>\n",
    "            <td style=\"text-align:center\">12</td>\n",
    "            <td style=\"text-align:center\">3</td>\n",
    "            <td style=\"text-align:center\">0</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"text-align:center\">1</td>\n",
    "            <td style=\"text-align:center\">3</td>\n",
    "            <td style=\"text-align:center\">8</td>\n",
    "            <td style=\"text-align:center\">1</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"text-align:center\">9</td>\n",
    "            <td style=\"text-align:center\">3</td>\n",
    "            <td style=\"text-align:center\">2</td>\n",
    "            <td style=\"text-align:center\">1</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"text-align:center\">5</td>\n",
    "            <td style=\"text-align:center\">8</td>\n",
    "            <td style=\"text-align:center\">7</td>\n",
    "            <td style=\"text-align:center\">0</td>\n",
    "        </tr>\n",
    "        <tr>\n",
    "            <td style=\"text-align:center\">2</td>\n",
    "            <td style=\"text-align:center\">7</td>\n",
    "            <td style=\"text-align:center\">8</td>\n",
    "            <td style=\"text-align:center\">0</td>\n",
    "        </tr>\n",
    "    </tbody>\n",
    "</table>\n",
    "\n",
    "This is the standard notation used in machine learning. Following this convention makes it easier for other machine learning practitioners to understand your scripts. This notation also allows us to concisely define supervised machine learning and classification described in the [Introduction](#1.-Introduction) section.\n",
    "\n",
    "#### Supervised learning\n",
    "Learn a function `f` that maps features to labels\n",
    "<p style=\"text-align:center;font-weight:bold\">$f(X) → y$</p>\n",
    "\n",
    "\\begin{align}\n",
    "y &= \\begin{bmatrix}\n",
    "    y_1 \\\\\n",
    "    \\vdots \\\\\n",
    "    y_i \\\\\n",
    "    \\vdots \\\\\n",
    "    y_{m}\n",
    "    \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "With classification, $y_i = \\{c_1, ..., c_k\\},where\\ k = number\\ of\\ classes$\n",
    "\n",
    "We only have two classes in our dataset so $y_i \\in \\{0=elliptical, 1=spiral\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scikit Learn\n",
    "We will be using the [scikit-learn](http://scikit-learn.org/) machine learning library for developing our models. As such, we need to prepare our dataset into a format that scikit-learn expects:\n",
    "\n",
    "1. Features and labels are **separate objects**\n",
    "2. Features and labels should be **numeric**\n",
    "3. Features and labels should be **`numpy` arrays**\n",
    "4. Features and labels should have **specific shapes**\n",
    "\n",
    "**Note**: `pandas` is built on top of `numpy` so we can create the features matrix as a pandas `DataFrame` and labels vector as a pandas `Series`.\n",
    "\n",
    "Let's prepare the dataset for scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Set a random seed number to repoduce our results\n",
    "seed = 19\n",
    "\n",
    "# Load the clean dataset\n",
    "df = pd.read_csv('data/galaxies-clean.csv')\n",
    "\n",
    "# Selected columns for modelling\n",
    "measurements = ['redshift', 'stellar_mass', 'deVRad_r', 'deVAB_r', 'expRad_r', 'expAB_r']\n",
    "magnitudes = ['mag_u', 'mag_g', 'mag_r', 'mag_i', 'mag_z']\n",
    "features = measurements + magnitudes\n",
    "\n",
    "# features matrix\n",
    "X = df[features]\n",
    "# labels vector\n",
    "y = df['class']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's verify our features matrix and labels vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Features matrix\n",
    "# expected output: pandas DataFrame, (number of rows, number of features)\n",
    "print('X: %s, %s' % (type(X), X.shape))\n",
    "\n",
    "# Labels vector \n",
    "# expected output: pandas Series, (number of rows,)\n",
    "print('y: %s, %s' % (type(y), y.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks great. Now we can start modelling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Classification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "### Splitting the data\n",
    "\n",
    "It's time to make the next big step in our analysis: splitting the data into training and test sets.\n",
    "\n",
    "- A **training set** is a random subset of the data that we use to train our models.\n",
    "- A **test set** is a random subset of the data (mutually exclusive from the training set) that we use to test our models.\n",
    "\n",
    "In machine learning we are always concerned that our models will **overfit** the data. i.e. The model learns the training set so well that it won't be able to handle examples it's never seen before. This is why it's important for us to build the model with the training set, but score it with a separate unseen testing set.\n",
    "\n",
    "Once we split the data into a training and test set, we should treat the test set like it no longer exists. We cannot use any information from the testing set to build our model or else we're cheating. The training dataset can also be split again to create a validation set to tune our models. It is important that the validation set is separate from the holdout test dataset.\n",
    "\n",
    "This process is illustrated in the figure below.\n",
    "\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Holdout set method to split the dataset</div>\n",
    "<img src=\"images/holdout.png\" />\n",
    "\n",
    "\n",
    "\n",
    "Now let's split the dataset using the `train_test_split` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Split the dataset into 80% training, 20% testing\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.8, random_state=seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the dataset split, we can start fitting models to our data. We have heard from our colleagues that they have had success using decision tree classifiers in their projects so let's start with those.\n",
    "\n",
    "Decision tree classifiers are incredibly simple in theory. In their simplest form, decision tree classifiers ask a series of Yes/No questions about the data — each time getting closer to finding out the class of each example — until they either classify the data set perfectly or simply can't differentiate a set of examples.\n",
    "\n",
    "An example decision tree classifier for approving loan applications is shown in the following figure.\n",
    "\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Decision tree example</div>\n",
    "<img src=\"images/decision_tree_example.png\" />\n",
    "\n",
    "Notice how the classifier asks yes/no questions about the data. e.g. whether the applicant owns a house so it can differentiate the records. \n",
    "\n",
    "Decision tree classifiers are *scale-invariant*, i.e. the scale of the features does not affect their performance unlike many machine learning models. In other words, it doesn't matter if our features range from 0 to 1 or 0 to 1,000; decision tree classifiers will work with them just the same.\n",
    "\n",
    "There are several [parameters](http://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html) that we can tune for decision tree classifiers and different [metrics](https://en.wikipedia.org/wiki/Evaluation_of_binary_classifiers) we can use to evaluate its performance. For now let's use a basic decision tree and the accuracy performance metric.\n",
    "\n",
    "\\begin{align}\n",
    "accuracy &= \\dfrac{correct\\ classifications}{total\\ number\\ of\\ classifications}\n",
    "\\end{align}\n",
    "\n",
    "Scikit learn provides a 4 step modelling pattern which makes it easy to switch in different models / algorithms for your dataset. This pattern is described in the code and comments below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Step 1: Import the model\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "# Step 2: Instantiate the model\n",
    "decision_tree_classifier = DecisionTreeClassifier(random_state=seed)\n",
    "\n",
    "# Step 3: Fit the model on data (i.e. train the model)\n",
    "decision_tree_classifier.fit(X_train, y_train)\n",
    "\n",
    "# Step 4: Generate predictions / scores\n",
    "decision_tree_classifier.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's great! Our model achieves 91% accuracy without much effort. We have already beaten the challenge to build a model with ≥ 70% accuracy.\n",
    "\n",
    "However, there's a catch. Depending on how our training and testing set is sampled, our model can achieve anywhere from 87% to 96% accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model_accuracies = []\n",
    "\n",
    "# Split the dataset differently and fit a model on this split, 1,000 times\n",
    "for i in range(1000):\n",
    "    X2_train, X2_test, y2_train, y2_test = train_test_split(X, y, train_size=0.8, random_state=i)\n",
    "    \n",
    "    decision_tree_classifier = DecisionTreeClassifier(random_state=seed)\n",
    "    decision_tree_classifier.fit(X2_train, y2_train)\n",
    "    classifier_accuracy = decision_tree_classifier.score(X2_test, y2_test)\n",
    "    model_accuracies.append(classifier_accuracy)\n",
    "\n",
    "sb.distplot(model_accuracies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a problem! The model performance varies a lot depending on the subset of the data it's trained on. This means the model is **overfitting**: the model learns to classify the training set so well that it doesn't generalize and perform well on data it hasn't seen before (i.e. the test set)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross-validation\n",
    "\n",
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Cross validation is a method of splitting the dataset to help address this problem. For this lesson, we will use the ***k*-fold cross-validation** method. This involves splitting the original data set into *k* subsets, use one of the subsets as the testing set and, the rest of the subsets are used as the training set. This process is repeated *k* times such that each subset is used as the testing set exactly once. This process is illustrated in Figure X.\n",
    "\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: k-fold cross-validation</div>\n",
    "<img src=\"images/k-fold.png\" />\n",
    "\n",
    "10-fold cross-validation is the most common choice so let's use that here. Performing 10-fold cross-validation on our data set looks something like this on a subset of 100 examples in our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.cross_validation import StratifiedKFold\n",
    "\n",
    "def plot_cv(cv, n_samples):\n",
    "    masks = []\n",
    "    for train, test in cv:\n",
    "        mask = np.zeros(n_samples, dtype=bool)\n",
    "        mask[test] = 1\n",
    "        masks.append(mask)\n",
    "        \n",
    "    plt.figure(figsize=(15, 15))\n",
    "    plt.imshow(masks, interpolation='none')\n",
    "    plt.ylabel('Fold')\n",
    "    plt.xlabel('Row #')\n",
    "\n",
    "# Use a subset of 100 examples\n",
    "subset = y[0:100]\n",
    "    \n",
    "plot_cv(StratifiedKFold(subset, n_folds=10), len(subset))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each square in the plot represents an example in our dataset.\n",
    "\n",
    "You'll notice that we used the `StratifiedKFold` function in the code above. Stratified means we keep the class percentage the same across all of the folds (~75% spiral), which is important for maintaining a representative subset of our dataset. i.e. we don't want to end up having 100% spiral galaxies in one of the folds.\n",
    "\n",
    "We can fit a decision tree classifier using 10-fold cross-validation with the `cross_val_score` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "decision_tree_classifier = DecisionTreeClassifier(random_state=seed)\n",
    "\n",
    "# cross_val_score returns a list of the scores, which we can visualize\n",
    "# to get a reasonable estimate of our classifier's performance\n",
    "cv_scores = cross_val_score(decision_tree_classifier, X_train, y_train, cv=10)\n",
    "sb.distplot(cv_scores)\n",
    "plt.title('Average score: {}'.format(np.mean(cv_scores)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks much better. We have a more consistent score (less variance) of our classifier's classification accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Learning curves allow us to evaluate the performance of our models as they are provided more training examples. i.e. They show us how or if our model learns. Plotting these curves also allows us to diagnose our models in order to gain insights on how we can improve them. More specifically, learning curves allows us to identify:\n",
    "\n",
    "- underfitting (high bias): model performs poorly on training and validation\n",
    "- overfitting (high variance): model performance is good on training but much poorer on validation\n",
    "\n",
    "Ideally, we want to develop a model that achieves similar performance on both the training and validation sets (good bias and variance trade off) as shown in the figure below.\n",
    "\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Learning curves</div>\n",
    "<img src=\"images/learning_curves.png\" />\n",
    "\n",
    "Let's plot the learning curve for our decision tree classifier using the `learning_curve` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate the learning curve scores\n",
    "train_sizes, train_scores, val_scores = learning_curve(estimator=decision_tree_classifier,\n",
    "                                                       X=X_train,\n",
    "                                                       y=y_train,\n",
    "                                                       train_sizes=np.linspace(0.1, 1.0, 20),\n",
    "                                                       cv=10,\n",
    "                                                       n_jobs=1)\n",
    "\n",
    "# Calculate the result averages and standard deviation\n",
    "train_mean = np.mean(train_scores, axis=1)\n",
    "train_std = np.std(train_scores, axis=1)\n",
    "val_mean = np.mean(val_scores, axis=1)\n",
    "val_std = np.std(val_scores, axis=1)\n",
    "\n",
    "# Plot the training learning curve\n",
    "plt.plot(train_sizes, train_mean,\n",
    "         color='green', marker='o',\n",
    "         markersize=5, label='training accuracy')\n",
    "\n",
    "plt.fill_between(train_sizes,\n",
    "                 train_mean + train_std,\n",
    "                 train_mean - train_std,\n",
    "                 alpha=0.15, color='green')\n",
    "\n",
    "# Plot the validation learning curve\n",
    "plt.plot(train_sizes, val_mean,\n",
    "         color='blue', linestyle='--',\n",
    "         marker='s', markersize=5,\n",
    "         label='validation accuracy')\n",
    "\n",
    "plt.fill_between(train_sizes,\n",
    "                 val_mean + val_std,\n",
    "                 val_mean - val_std,\n",
    "                 alpha=0.15, color='blue')\n",
    "\n",
    "# Plot settings\n",
    "plt.xlabel('Number of training samples')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.legend(loc='lower right')\n",
    "plt.ylim([0.8, 1])\n",
    "plt.tight_layout()\n",
    "# plt.savefig('images/learning_curve.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the learning curve we can see that the decision tree classifier overfits the training dataset. The training accuracy remains constant at 100% (fits the data perfectly) while the average validation scores steadily increase from 87% to 93% as the number of training samples increases.\n",
    "\n",
    "Now that we know our decision tree classifier overfits the dataset, let's do something about it!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameter tuning\n",
    "\n",
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Every machine learning model comes with a variety of parameters to tune, and these parameters can be vitally important to the performance of our classifier. For example, using the default parameters for decision trees has resulted in our classifier overfitting the training dataset. \n",
    "\n",
    "As an extreme example, let's limit the depth of our decision tree classifier to a `max_depth = 1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set the max_depth of the decision tree to 1\n",
    "decision_tree_classifier = DecisionTreeClassifier(max_depth=1, random_state=seed)\n",
    "\n",
    "# Calculate the cross validation score\n",
    "cv_scores = cross_val_score(decision_tree_classifier, X_train, y_train, cv=10)\n",
    "\n",
    "sb.distplot(cv_scores, kde=False)\n",
    "plt.title('Average score: {}'.format(np.mean(cv_scores)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The classification accuracy falls tremendously to 86%.\n",
    "\n",
    "We need to use a systematic method to discover the best parameters for our model and dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Grid search\n",
    "\n",
    "The most common method for model parameter (hyperparameter) tuning is grid search. The idea behind grid search is to explore a range of parameters and find the best-performing parameter combination. Focus your search on the best range of parameters, then repeat this process several times until the best parameters are discovered.\n",
    "\n",
    "Let's tune our decision tree classifier using the `GridSearchCV` function. We'll stick to only two parameters for now (`max_depth` and `max_features`), but it's possible to simultaneously explore dozens of parameters if we want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "decision_tree_classifier = DecisionTreeClassifier(random_state=seed)\n",
    "\n",
    "parameter_grid = {'max_depth': [1, 2, 3, 4, 5, 6],\n",
    "                  'max_features': [7, 8, 9, 10, 11]}\n",
    "\n",
    "cross_validation = StratifiedKFold(y_train, n_folds=10, random_state=seed)\n",
    "\n",
    "grid_search = GridSearchCV(decision_tree_classifier,\n",
    "                           param_grid=parameter_grid,\n",
    "                           cv=cross_validation)\n",
    "\n",
    "grid_search.fit(X_train, y_train)\n",
    "print('Best score: {}'.format(grid_search.best_score_))\n",
    "print('Best parameters: {}'.format(grid_search.best_params_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualize the grid search to see how the parameters interact."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid_visualization = []\n",
    "\n",
    "for grid_pair in grid_search.grid_scores_:\n",
    "    grid_visualization.append(grid_pair.mean_validation_score)\n",
    "\n",
    "grid_visualization = np.array(grid_visualization)\n",
    "grid_visualization.shape = (6, 5)\n",
    "sb.heatmap(grid_visualization, cmap='Blues')\n",
    "plt.xticks(np.arange(5) + 0.5, grid_search.param_grid['max_features'])\n",
    "plt.yticks(np.arange(6) + 0.5, grid_search.param_grid['max_depth'][::-1])\n",
    "plt.xlabel('max_features')\n",
    "plt.ylabel('max_depth')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a better sense of the parameter space: a `max_depth > 4` is needed for the decision tree perform reasonably well. Tuning the `max_features` doesn't seem to make as much difference as we achieved the best performance with using `9` out of the 11 features.\n",
    "\n",
    "An alternative way to visualise the performance of your parameter tuning results is to plot [validation curves](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.validation_curve.html).\n",
    "\n",
    "Let's go ahead and use a broader grid search to find the best settings for more parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Instantiate the decision tree classifier\n",
    "decision_tree_classifier = DecisionTreeClassifier(random_state=seed)\n",
    "\n",
    "# Hyperparameter values to explore\n",
    "parameter_grid = {'criterion': ['gini', 'entropy'],\n",
    "                  'splitter': ['best', 'random'],\n",
    "                  'max_depth': [4, 5, 6],\n",
    "                  'max_features': [7, 8, 9, 10, 11],\n",
    "                  'min_samples_split': [2, 4, 6, 8]}\n",
    "\n",
    "# Create a cross validation object using StratifiedKFold\n",
    "cross_validation = StratifiedKFold(y_train, n_folds=10, random_state=seed)\n",
    "\n",
    "# Instantiate the grid search using GridSearchCV\n",
    "grid_search = GridSearchCV(decision_tree_classifier,\n",
    "                           param_grid=parameter_grid,\n",
    "                           cv=cross_validation)\n",
    "\n",
    "# Fit the decision tree using grid search\n",
    "grid_search.fit(X_train, y_train)\n",
    "\n",
    "# Display the best model scores and parameters\n",
    "print('Best score: {}'.format(grid_search.best_score_))\n",
    "print('Best parameters: {}'.format(grid_search.best_params_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the best classifier from the grid search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "decision_tree_classifier = grid_search.best_estimator_\n",
    "# Examine the model\n",
    "decision_tree_classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the learning curve for the best classifier for evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate the learning curve scores\n",
    "\n",
    "# Calculate the result averages and standard deviation\n",
    "# train_mean = \n",
    "# train_std = \n",
    "# val_mean = \n",
    "# val_std = \n",
    "\n",
    "# Plot the training learning curve\n",
    "plt.plot(train_sizes, train_mean,\n",
    "         color='green', marker='o',\n",
    "         markersize=5, label='training accuracy')\n",
    "\n",
    "plt.fill_between(train_sizes,\n",
    "                 train_mean + train_std,\n",
    "                 train_mean - train_std,\n",
    "                 alpha=0.15, color='green')\n",
    "\n",
    "# Plot the validation learning curve\n",
    "plt.plot(train_sizes, val_mean,\n",
    "         color='blue', linestyle='--',\n",
    "         marker='s', markersize=5,\n",
    "         label='validation accuracy')\n",
    "\n",
    "plt.fill_between(train_sizes,\n",
    "                 val_mean + val_std,\n",
    "                 val_mean - val_std,\n",
    "                 alpha=0.15, color='blue')\n",
    "\n",
    "# Plot settings\n",
    "plt.xlabel('Number of training samples')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.legend(loc='lower right')\n",
    "plt.ylim([0.8, 1])\n",
    "plt.tight_layout()\n",
    "# plt.savefig('images/learning_curve.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This learning curve looks better. The gap between the training and validation accuracy has reduced with training score of 95% and a validation score of 93% which suggests the model is no longer overfitting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reporting\n",
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Score the classifier on the test dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Get the best estimator from the grid search\n",
    "decision_tree_classifier = grid_search.best_estimator_\n",
    "\n",
    "# Score the estimator on the test dataset\n",
    "decision_tree_classifier.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our best classifier achieves a test accuracy score of 92.7% which is almost on par with our validation score of 93%. This gives us confidence that our experiment design and trained model is providing realistic results.\n",
    "\n",
    "We can generate a confusion matrix using the `confusion_matrix` function to get more detailed classification results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "# Generate the predictions\n",
    "y_pred = decision_tree_classifier.predict(X_test)\n",
    "\n",
    "# Calculate the confusion matrix\n",
    "cm = confusion_matrix(y_test, y_pred)\n",
    "\n",
    "# Create a DataFrame for plotting\n",
    "labels = ['elliptical', 'spiral']\n",
    "df_cm = pd.DataFrame(cm, columns=labels, index=labels)\n",
    "\n",
    "# Plot the confusion matrix\n",
    "sb.heatmap(df_cm, annot=True, fmt='g', cmap='Blues', annot_kws={\"size\": 15})\n",
    "plt.title('Confusion Matrix', fontsize=20)\n",
    "plt.xlabel('Actual', fontsize=15)\n",
    "plt.xticks(fontsize=13)\n",
    "plt.ylabel('Predicted', fontsize=15)\n",
    "plt.yticks(fontsize=13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The confusion matrix provides more detailed information of where the classifier made mistakes (i.e. where it was confused). Here we can see that the model predicted 12 galaxies as spiral that were actually elliptical and 14 elliptical galaxies that were actually spiral. Overall, it looks like the decision tree classifier has performed quite well on test dataset.\n",
    "\n",
    "We can also visualize the decision tree with [GraphViz](http://www.graphviz.org/) to see how it's making the classifications:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sklearn.tree as tree\n",
    "from sklearn.externals.six import StringIO\n",
    "\n",
    "with open('images/galaxy_tree.dot', 'w') as out_file:\n",
    "    out_file = tree.export_graphviz(decision_tree_classifier, out_file=out_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# The exported decision tree graph displays features by their index. e.g. X[0], ... X[10]\n",
    "# This snippet will replace the index names with the actual feature names\n",
    "replacements = {}\n",
    "for index, feature in enumerate(features):\n",
    "    key = 'X[%d]' % (index)\n",
    "    replacements[key] = feature\n",
    "\n",
    "with open('images/galaxy_tree.dot') as infile, open('images/galaxy_tree_names.dot', 'w') as outfile:\n",
    "    for line in infile:\n",
    "        for src, target in replacements.items():\n",
    "            line = line.replace(src, target)\n",
    "        outfile.write(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have [GraphViz](http://www.graphviz.org/) installed then you can generate a image of the decision tree with the following command:\n",
    "\n",
    "    dot -Tpng images/galaxy_tree_names.dot -o images/galaxy_tree_names.png"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/galaxy_tree_names.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright! We finally have our demo classifier. Let's create some more visuals of its performance so we have something to put in our report."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dt_scores = cross_val_score(decision_tree_classifier, X, y, cv=10)\n",
    "sb.boxplot(dt_scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot is not very useful with just one model so let's train another.\n",
    "\n",
    "A common problem with decision trees is that they're prone to overfitting. They can add complex rules to the point that they classify the training set near-perfectly (as we saw early with our learning curves), but fail to generalise well to data they have not seen before.\n",
    "\n",
    "**Random Forest classifiers** work around this limitation by creating many decision trees (i.e. a forest), each trained on random subsets of training samples and features. These decision trees are then combined to make a more accurate classification. Let's see if a random forest classifier works better for our datatset. \n",
    "\n",
    "The great part about scikit-learn is that the modelling pattern of training, testing, parameter tuning, etc. process is the same for all models. We only need to plug in the new classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Step 1: Import the class\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# Step 2: Instantiate the estimator\n",
    "random_forest_classifier = RandomForestClassifier(random_state=seed)\n",
    "\n",
    "# Hyperparameter values to explore\n",
    "parameter_grid = {'n_estimators': [5, 10, 25, 50],\n",
    "                  'criterion': ['gini', 'entropy'],\n",
    "                  'max_features': [1, 2, 3, 4],\n",
    "                  'warm_start': [True, False]}\n",
    "\n",
    "# Create a cross validation object using StratifiedKFold\n",
    "cross_validation = StratifiedKFold(y_train, n_folds=10, random_state=seed)\n",
    "\n",
    "# Instantiate the grid search using GridSearchCV \n",
    "grid_search = GridSearchCV(random_forest_classifier,\n",
    "                           param_grid=parameter_grid,\n",
    "                           cv=cross_validation)\n",
    "\n",
    "# Step 3 & 4: Fit the estimator on data (i.e. train the model) and generate predictions\n",
    "grid_search.fit(X_train, y_train)\n",
    "\n",
    "# Display the best model scores and parameters\n",
    "print('Best score: {}'.format(grid_search.best_score_))\n",
    "print('Best parameters: {}'.format(grid_search.best_params_))\n",
    "\n",
    "grid_search.best_estimator_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The random forest classifier achieved a 95% cross validation accuracy score. Find out how well it performs on the test dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Get the best estimator from the grid search\n",
    "random_forest_classifier = grid_search.best_estimator_\n",
    "# random_forest_classifier =\n",
    "\n",
    "# Score the estimator on the test dataset\n",
    "random_forest_classifier.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random forests performs better than the decision tree classifier on the test dataset:\n",
    "\n",
    "- decision trees ~ 92%\n",
    "- random forests ~ 96%\n",
    "\n",
    "We can also plot and compare their cross validation performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rf_df = pd.DataFrame({'accuracy': cross_val_score(random_forest_classifier, X_train, y_train, cv=10),\n",
    "                       'classifier': ['Random Forest'] * 10})\n",
    "dt_df = pd.DataFrame({'accuracy': cross_val_score(decision_tree_classifier, X_train, y_train, cv=10),\n",
    "                      'classifier': ['Decision Tree'] * 10})\n",
    "both_df = rf_df.append(dt_df)\n",
    "\n",
    "sb.boxplot(x='classifier', y='accuracy', data=both_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random forests also performs better on average than decisions trees with cross validation. Random forest classifiers can perform particularly well when there's hundreds of possible features to model (we only have 11 features in dataset). In other words, there was not much room for improvement for random forests over decision trees with this dataset.\n",
    "\n",
    "To report more detailed information of our model, let's see what features were considered the most important for both classifiers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Decision tree features\n",
    "sorted(zip(decision_tree_classifier.feature_importances_, features), reverse=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Random forests features\n",
    "sorted(zip(random_forest_classifier.feature_importances_, features), reverse=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both classifiers rank the importance of the features in a similar order with the exception of the magnitudes. The stellar mass, ellipticity, radius and redshift of the galaxies were the most important features. \n",
    "\n",
    "If you want to explore this further you could experiment with removing the `mag_u`, `mag_i`, `mag_g` features as they scored a feature importance of `0` for the decision tree classifier to see how it affects classification performance. Additionally, you can experiment with various [feature selection methods](http://scikit-learn.org/stable/modules/feature_selection.html) for systemtically selecting a reduced feature set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Other algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many [classification algorithms](http://scikit-learn.org/stable/supervised_learning.html#supervised-learning) we can experiment with on our dataset. Try fitting a [Logistic Regression](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html) and [Multi Layer Perceptron (MLP - Neural Network)](http://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPClassifier.html) model on the dataset. It is up to you whether you want to perform parameter tuning on these models. \n",
    "\n",
    "Generate boxplots so we can compare the performance with the decision tree and random forest models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the algorithms and cross validation score function\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "\n",
    "# Initialise the algorithms\n",
    "logreg_classifier = LogisticRegression(random_state=seed)\n",
    "mlp_classifier = MLPClassifier(random_state=seed)\n",
    "\n",
    "# Calculate the cross validation scores\n",
    "logreg_df = pd.DataFrame({'accuracy': cross_val_score(logreg_classifier, X_train, y_train, cv=10),\n",
    "                       'classifier': ['Logistic Regression'] * 10})\n",
    "mlp_df = pd.DataFrame({'accuracy': cross_val_score(mlp_classifier, X_train, y_train, cv=10),\n",
    "                       'classifier': ['Multilayer Perception'] * 10})\n",
    "\n",
    "# Merge all classifier results\n",
    "classifiers_df = rf_df.append([dt_df, logreg_df, mlp_df])\n",
    "\n",
    "# Generate a box plot comparing the different algorithms\n",
    "sb.boxplot(x='classifier', y='accuracy', data=classifiers_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pick your favourite model and use it to predict the morphology of galaxies that have been labelled as  uncertain in the dataset. \n",
    "\n",
    "- How many spiral vs. elliptical galaxies does the model predict?\n",
    "- How do the model predictions compare with the debiased labelled probabilities?\n",
    "- How could you use this model and results in practice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Your colleague has heard about your recent success in applying machine learning for galaxy classification. They mention they have a dataset containing _ugriz_ photomerty of galaxies but no redshift measurements. They would like you to develop a machine learning model to estimate the redshift values.\n",
    "\n",
    "This is a regression problem. Regression is a supervised learning approach but instead of predicting the category of an example, we are predicting a continuous value. i.e. $y_i \\in\\ \\mathbb{R}$. For example, the predicted redshift of galaxy `i` is 0.2312.\n",
    "\n",
    "Developing a regression model follows the same process as building a classification model with the exception of using a different metric for evaluating the model performance. i.e. we can't use classification accuracy for evaluating our regression model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Prepare the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load and prepare the dataset for using the scikit-learn library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Set a random seed number to reproduce our results\n",
    "seed = 11\n",
    "\n",
    "# Clean dataset file\n",
    "clean_dataset = 'data/galaxies-clean.csv'\n",
    "\n",
    "# 1. Load the dataset using pandas read_csv function\n",
    "# df = \n",
    "\n",
    "# 2. Select the columns of interest for modelling\n",
    "# features = \n",
    "\n",
    "# 3. Create the features matrix as X\n",
    "# X = \n",
    "\n",
    "# 4. Create the labels vector as y\n",
    "# y = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Split the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a training and test dataset using the `train_test_split` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the function\n",
    "\n",
    "# Split the dataset into X_train, X_test, y_train, y_test\n",
    "# Use a training dataset size of 80%\n",
    "# X_train, X_test, y_train, y_test = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will fit the regression models without any parameter tuning using the following three models:\n",
    "\n",
    "- [Linear Regression](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html)\n",
    "- [RandomForestRegressor](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html)\n",
    "- [MLPRegressor](http://scikit-learn.org/stable/modules/generated/sklearn.neural_network.MLPRegressor.html)\n",
    "\n",
    "We will measure model performance using the root mean squared error (RMSE) metric. This can be calculated by using the NumPy function `np.sqrt` and the scikit-learn `mean_squared_error` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Step 1: Import the classes\n",
    "\n",
    "# Step 2: Instantiate the estimators\n",
    "# lr =\n",
    "# rf =\n",
    "# mlp = \n",
    "\n",
    "# Step 3: Fit the estimators on data (i.e. train the models)\n",
    "\n",
    "# Step 4: Generate predictions\n",
    "# y_pred_m1 =\n",
    "# y_pred_m2 =\n",
    "# y_pred_m3 =\n",
    "\n",
    "# Calculate the Root Mean Squared Error (RMSE)\n",
    "# np.sqrt(mean_squared_error(...))\n",
    "# m1_score = \n",
    "# m2_score = \n",
    "# m3_score = \n",
    "\n",
    "# Display the model scores\n",
    "print('Linear regression: %.3f' % (m1_score))\n",
    "print('Random forest regressor: %.3f' % (m2_score))\n",
    "print('Multi layer perceptron: %.3f' % (m3_score))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression model is off by an average of 0.022 redshift while the random forest regression has an RMSE of 0.026. The multi layer perceptron (without any parameter tuning) achieves a poorer RMSE of 0.079.\n",
    "\n",
    "Let's examine the correlation between the predicted and actual redshift values for both models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# calculate the correlation between the labels and predictions\n",
    "m1_corr = np.corrcoef(y_test, y_pred_m1)[0][1]\n",
    "m2_corr = np.corrcoef(y_test, y_pred_m2)[0][1]\n",
    "m3_corr = np.corrcoef(y_test, y_pred_m3)[0][1]\n",
    "\n",
    "print('Linear regression: %.3f' % (m1_corr))\n",
    "print('Random forest regressor: %.3f' % (m2_corr))\n",
    "print('Multi layer perceptron: %.3f' % (m3_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's impressive. The linear regression model has a correlation value of 0.889 while the random forests achieves 0.834. However, the multi layer perceptron performs poorly with a correlation of -0.183.\n",
    "\n",
    "Let's look at the generated coefficients and intercept values for the linear regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(lr.intercept_)\n",
    "print(lr.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# pair coefficients with feature names\n",
    "list(zip(features, lr.coef_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that the linear regression model for estimating the redshift of a galaxy is:\n",
    "\n",
    "<p style=\"text-align:center;font-weight:bold\">$redshift = -0.027(mag_u) + 0.255(mag_g) - 0.138 (mag_r) - 0.373(mag_i) + 0.302(mag_z) - 0.294$</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature normalisation\n",
    "\n",
    "An important data preparation step for many machine learning models is to [normalise / standarised the feature values](http://scikit-learn.org/stable/modules/preprocessing.html) in our dataset. e.g. we can scale the magnitude values to a smaller and standardised range of [0, 1]. Note: this is not required for decision trees and random forests as those models are scale invariant.\n",
    "\n",
    "We can normalise our features to a range of [0, 1] using the scikit-learn `MinMaxScaler()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler \n",
    "\n",
    "# Fit the scaler on the training dataset ONLY\n",
    "scaler = MinMaxScaler().fit(X_train)\n",
    "\n",
    "# Transform both the training and test datasets\n",
    "X_train_scaled = scaler.transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)\n",
    "\n",
    "# Print out the new scaled values\n",
    "print('Original scale:')\n",
    "print('  X_train %.3f - %.3f' % (np.min(np.min(X_train)), np.max(np.max(X_train))))\n",
    "print('  X_test %.3f - %.3f' % (np.min(np.min(X_test)), np.max(np.max(X_test))))\n",
    "print('Transformed scale:')\n",
    "print('  X_train %.3f - %.3f' % (np.min(X_train_scaled), np.max(X_train_scaled)))\n",
    "print('  X_test %.3f - %.3f' % (np.min(X_test_scaled), np.max(X_test_scaled)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now fit a linear regression and multi layer perceptron model on the scaled dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Instantiate the estimators\n",
    "# lr_scaled =\n",
    "# mlp_scaled =\n",
    "\n",
    "# Fit the estimators\n",
    "\n",
    "# Generate predictions\n",
    "# y_pred_m1 =\n",
    "# y_pred_m2 =\n",
    "\n",
    "# Calculate the Root Mean Squared Error (RMSE)\n",
    "# np.sqrt(mean_squared_error(...))\n",
    "# m1_score =\n",
    "# m2_score =\n",
    "\n",
    "# calculate the correlation coefficient\n",
    "# m1_corr =\n",
    "# m2_corr =\n",
    "\n",
    "# print the results\n",
    "print('Linear regresssion scaled')\n",
    "print('  RMSE: %.3f' % (m1_score))\n",
    "print('  correlation: %.3f' % (m1_corr))\n",
    "print('Multi layer perception scaled')\n",
    "print('  RMSE: %.3f' % (m2_score))\n",
    "print('  correlation: %.3f' % (m2_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression model achieved the same performance on the scaled and unscaled dataset (0.022 RMSE). However, the multi layer perceptron model has an improved RMSE of 0.040 (previously 0.079) and correlation coefficient of 0.554 (previously -0.183) on the scaled dataset.\n",
    "\n",
    "Normalisation might not always improve performance but it is common practice to do so before training certain types of models. You should also think about your dataset and what features you are normalising. e.g. are the features measured on a linear scale? Should magnitude be normalised of should the raw fluxes be normalised instead?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature engineering\n",
    "\n",
    "Another approach to improve the performance of our models is to try new features based on our knowledge of the domain. This is known as feature engineering.\n",
    "\n",
    "We are currently using the raw magnitude bands as features to estimate a galaxy's redshift. However, it might make more sense to use colour features that measure the ratio of flux in neighbouring filters. This is equivalent to subtracting the magnitudes of the neighbouring filters. The key to photometric red shift is that a red shifted galaxy will have different observed colors to what it would have at red shift zero. i.e. galaxies at higher redshift tend to be redder in colour. \n",
    "\n",
    "With this knowledge, let's build a model with 4 engineered features by subtracting neighbouring bands in the `ugriz` magnitude channels:\n",
    "\n",
    "- $mag_u - mag_g$\n",
    "- $mag_g - mag_r$\n",
    "- $mag_r - mag_i$\n",
    "- $mag_i - mag_z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make a copy of the features matrix\n",
    "X_new = X.copy()\n",
    "\n",
    "# Create new features\n",
    "# subtract neighbouring magnitudes\n",
    "X_new['u-g'] = X_new['mag_u'] - X_new['mag_g']\n",
    "X_new['g-r'] = X_new['mag_g'] - X_new['mag_r']\n",
    "X_new['r-i'] = X_new['mag_r'] - X_new['mag_i']\n",
    "X_new['i-z'] = X_new['mag_i'] - X_new['mag_z']\n",
    "\n",
    "# Remove the old columns\n",
    "X_new.drop(['mag_u', 'mag_g', 'mag_r', 'mag_i', 'mag_z'], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit a linear regression and a multi layer perceptron model with the new features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Split the dataset into X_train, X_test, y_train, y_test\n",
    "# X_train_new, X_test_new, y_train_new, y_test_new = \n",
    "\n",
    "# Instantiate the estimators\n",
    "# lr_new =\n",
    "# mlp_new =\n",
    "\n",
    "# Fit the estimators\n",
    "\n",
    "# Generate predictions\n",
    "# m1_pred =\n",
    "# m2_pred =\n",
    "\n",
    "# Calculate the RMSE\n",
    "# m1_score =\n",
    "# m2_score =\n",
    "\n",
    "# Calculate the correlation coefficient\n",
    "# m1_corr = \n",
    "# m2_corr = \n",
    "\n",
    "# Print the result\n",
    "print('Linear regresssion new features')\n",
    "print('  RMSE: %.3f' % (m1_score))\n",
    "print('  correlation: %.3f' % (m1_corr))\n",
    "print('Multi layer perceptron new features')\n",
    "print('  RMSE: %.3f' % (m2_score))\n",
    "print('  correlation: %.3f' % (m2_corr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear regression model achieved a poorer RMSE of 0.025 (previously 0.022) and correlation coefficient of 0.846 (previously 0.888). However, the multi layer perceptron model achieved improved performance with a RMSE of 0.031 (previously on the scaled dataset 0.040) and a correlation coefficient of 0.753 (previously 0.554) with the new features.\n",
    "\n",
    "It is up to you to decide on how long you want to spend coming up with and experiment with new features. The amount of time you spend on feature engineering is often based on your desired model performance and project requirements (e.g. deadlines)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge\n",
    "\n",
    "Apply the steps used in the [Step 5: Classification](#Step-5:-Classification) section for cross-validation, generating learning curves, parameter tuning and reporting on your regression models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are there any ways you can improve the performance? Can you come up with better features or try other [regression algorithms](http://scikit-learn.org/stable/supervised_learning.html#supervised-learning) for this dataset?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Clustering is an unsupervised learning algorithm that groups similar objects together. How these objects are grouped is based upon what features are provided to the algorithm. Recall that unsupervised learning algorithms model the underlying structure of the data `X`. No corresponding labels `y` are available / provided."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### K-means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "K-means is a simple and popular clustering algorithm. The `K` in its name refers to the number of clusters we want to generate from the dataset. For example, `K = 3` means we want to generate 3 clusters. \n",
    "\n",
    "The K-means algorithm:\n",
    "\n",
    "1. randomly place `K` points into the feature space. These serve as the initial centroids.\n",
    "2. Assign each example to the group with the closest centroid (Euclidian distance).\n",
    "3. When all examples have been assigned, recalculate the positions of the `K` centroids.\n",
    "4. Repeat steps 2 and 3 until the centroids no longer change.\n",
    "\n",
    "An example output of generating 3 clusters for a dataset is shown in figure below. The `x` markers show the location of the final centroids used to group examples.\n",
    "\n",
    "![K-means](images/kmeans.png)\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: K-means clustering example</div>\n",
    "\n",
    "**Note:** K-means does not work for categorical features. If you want to cluster categorical features then consider other methods such as [k-modes](https://pypi.python.org/pypi/kmodes/).\n",
    "\n",
    "Let's try clustering galaxies based on the `u-r` and `g-i` colour features using the K-means algorithm.\n",
    "\n",
    "First, let's load and prepare our dataset again for this section of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Set a random seed number to repoduce our results\n",
    "seed = 19\n",
    "\n",
    "# Clean dataset file\n",
    "clean_dataset = 'data/galaxies-clean.csv'\n",
    "\n",
    "# Load the dataset\n",
    "df = pd.read_csv(clean_dataset)\n",
    "\n",
    "# Calculate the new features\n",
    "# df['u-r'] = \n",
    "# df['g-i'] = \n",
    "\n",
    "# Select columns of interest\n",
    "\n",
    "# Create the features matrix\n",
    "# X = \n",
    "\n",
    "# Create the labels vector\n",
    "# We're not using labels in modelling but rather for plotting\n",
    "# y = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split the dataset. Remember, we need to place the test dataset aside and treat it as it doesn't exist. We will be training a classification model based on insights gained from unsupervised learning later in the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Split the dataset into 80% training, 20% testing\n",
    "# X_train, X_test, y_train, y_test = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the K-means algorithm on the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the algorithm \n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# Specify the number of clusters\n",
    "K = 2\n",
    "\n",
    "# Instantiate the algorithm\n",
    "kmeans = KMeans(n_clusters=K, random_state=seed)\n",
    "\n",
    "# It is important to normalise the data for K-means to have zero mean\n",
    "# To do this we will use the StandardScaler function\n",
    "X_train_scaled = StandardScaler().fit_transform(X_train)\n",
    "\n",
    "# Fit the algorithm\n",
    "kmeans.fit(X_train_scaled) \n",
    "\n",
    "# predict the cluster / group for each example\n",
    "y_pred = kmeans.predict(X_train_scaled)\n",
    "\n",
    "# Generate plots\n",
    "plt.figure(figsize=(13, 5))\n",
    "\n",
    "# Predicted clusters\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.xlim(-3, 3)\n",
    "plt.title('Predicted clusters')\n",
    "plt.scatter(X_train_scaled[:,0], X_train_scaled[:,1], alpha=0.15, c=y_pred, cmap=plt.cm.Paired) \n",
    "# Plot the centroids\n",
    "plt.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1], marker='X', s=500, c='black')\n",
    "\n",
    "# Real labels\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.xlim(-3, 3)\n",
    "plt.title('Labels')\n",
    "plt.scatter(X_train_scaled[:,0], X_train_scaled[:,1], alpha=0.15, c=y_train, cmap=plt.cm.Paired) \n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generated K-means clusters don't match up well with the galaxy labels. That's fine because we aren't explicitly training a model to group galaxies by their morphology. We're asking the model to group galaxies based on the features provided (i.e. colour bands).\n",
    "\n",
    "We can also experiment different values for K and see how the clusters change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Specify a list of cluster numbers\n",
    "K_values = [2, 3, 4]\n",
    "\n",
    "# Initialise plots\n",
    "plt.figure(figsize=(13, 4))\n",
    "\n",
    "for index, K in enumerate(K_values):\n",
    "    kmeans = KMeans(n_clusters=K, random_state=seed)\n",
    "    \n",
    "    # Fit the algorithm\n",
    "    kmeans.fit(X_train_scaled) \n",
    "\n",
    "    # predict the cluster / group for each example\n",
    "    y_pred = kmeans.predict(X_train_scaled)\n",
    "\n",
    "    # plot results\n",
    "    plt.subplot(1, 3, index+1)\n",
    "    plt.xlim(-3, 3)\n",
    "    plt.title('K = %d' % (K))\n",
    "    plt.scatter(X_train_scaled[:,0], X_train_scaled[:,1], alpha=0.15, c=y_pred, cmap=plt.cm.Paired) \n",
    "    # Plot the centroids\n",
    "    plt.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1], marker='X', s=500, c='black')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what `K` value is best suited for our dataset? Since we have 2 classes in our dataset, we naturally want to group the data into 2 clusters. However, if we are open to consider different numbers of clusters, then it might actually make more sense to use `K=3` to generate 3 clusters, where the additional / middle cluster represents the 'green valley' (a transformative phase) where galaxies are not quite elliptical or spiral but a bit of both.\n",
    "\n",
    "If we had no labels, then we don't know how many clusters would best represent our dataset beforehand. Whether these clusters make sense or not will be based on trial and error. We would use our domain knowledge and inspect objects in these clusters to identify similarities.  \n",
    "\n",
    "However, a more systematic way to valiadate the number of clusters is to use the [elbow method](https://bl.ocks.org/rpgove/0060ff3b656618e9136b). This involves calculating the sum of squared errors (distance) between each data point and their assigned cluster centroid. The distance scores are plotted for different values of `K`. If the plot (aka scree plot) looks like an arm then the **elbow on the arm** is the value of `K` that is chosen as shown in the figure below.\n",
    "\n",
    "![K-means](images/elbow_method.png)\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: Elbow method to choose the number of clusters</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have only clustered our dataset using two features: `u-r` and `g-i`. Let's try it with all the original features and with `K = 2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create the features matrix\n",
    "# Be sure to exclude the id and class columns\n",
    "# X =\n",
    "\n",
    "# Split the dataset into 80% training, 20% testing\n",
    "# X_train, X_test, y_train, y_test \n",
    "\n",
    "# Specify the number of clusters\n",
    "# K = \n",
    "\n",
    "# Instantiate the KMeans algorithm\n",
    "# kmeans = \n",
    "\n",
    "# Normalise the dataset using StandardScaler\n",
    "# X_train_scaled = \n",
    "# X_test_scaled = \n",
    "\n",
    "# Fit the algorithm\n",
    "\n",
    "# predict the cluster / group for each example\n",
    "# y_pred = \n",
    "\n",
    "# Print the number of objects grouped in each cluster\n",
    "print('Cluster 0: %d' % (np.sum(y_pred == 0)))\n",
    "print('Cluster 1: %d' % (np.sum(y_pred == 1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like we have more galaxies grouped into Cluster 1. Let's examine a sample of records from each cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# display the first 10 galaxies in cluster 0\n",
    "X.iloc[X_train[y_pred == 0].index].head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# display the first 10 galaxies in cluster 1\n",
    "X.iloc[X_train[y_pred == 1].index].head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have generated 2 clusters using the K-means algorithm using 11 features (dimensions) in our dataset. But how can we visualise a model with more than 3 dimensions? We can explore dimensionality reduction methods to help us with this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit-learn offers many [clustering algorithms](http://scikit-learn.org/stable/modules/clustering.html) that you can easily and quickly experiment with on your datasets as shown below. Experiment with some of these algorithms and compare the generated clusters with the K-means algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Dimensionality reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dimensionality reduction is the process of reducing the number of features (dimensions) while trying to preserve as much structure in the dataset as possible. This is a useful method:\n",
    "\n",
    "- for visualising data with ≥ 3 dimensions\n",
    "- as a pre-processing step for other machine learning algorithms such as regression and clustering\n",
    "\n",
    "There are two approaches to performing dimensionality reduction:\n",
    "\n",
    "1. [Feature selection](http://scikit-learn.org/stable/modules/feature_selection.html): choose a subset of the original feature set\n",
    "2. Feature extraction: construct a new and smaller feature set. e.g. from a combination of the original features.\n",
    "\n",
    "For this notebook, we will experiment with a method called [Principal Components Analysis (PCA)](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html). Scikit-learn provides feature extraction methods in two modules: \n",
    "\n",
    "- [Feature extraction](http://scikit-learn.org/stable/modules/feature_extraction.html) - methods targetted towards extracting features from text and images \n",
    "- [Matrix decomposition](http://scikit-learn.org/stable/modules/classes.html#module-sklearn.decomposition) - dimensionality reduction techniques such as PCA, factor analysis etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Principal Components Analysis (PCA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Principal Components Analysis (PCA) is a popular and commonly used dimensionality reduction method. The objective of PCA is to reduce a dataset to it's basic (principal) components stripping away unnecessary parts. These components represent the underlying structure of the dataset, indicating the direction of most variance (the direction where data is most spread out).\n",
    "\n",
    "Here is a excellent [blog article explaining the intuition and inner workings of PCA](https://georgemdallas.wordpress.com/2013/10/30/principal-component-analysis-4-dummies-eigenvectors-eigenvalues-and-dimension-reduction/). \n",
    "\n",
    "The PCA algorithm:\n",
    "\n",
    "1. Normalise the dataset to have zero mean\n",
    "2. Calculate the co-variance matrix\n",
    "3. Calculate the eigenvalues and eigenvectors\n",
    "4. Choose the top `k` components to form the new feature set\n",
    "5. Transform / project the original dataset based on these `k` eigenvectors\n",
    "\n",
    "An example of performing PCA on a 2-dimensional dataset is shown in the figure below.\n",
    "\n",
    "<img src=\"images/pca.png\" />\n",
    "<div style=\"text-align:center;font-weight:bold\">Figure: PCA example - 1st component (red) and 2nd component (green)</div>\n",
    "\n",
    "A commonly used heurestic is to choose the `top k` components that can represent (explain) 90%, 95% or 99% of the variance of the dataset. In the above example, the 1st principal component (red) explains 99.6% of the variance so we could just use the `top` component and discard the 2nd. We can then use this principal component to transform (project) our original 2D dataset into a 1D dataset.\n",
    "\n",
    "### Visualisation\n",
    "\n",
    "First, let's use PCA to visualise our K-means model by only using the top 2 components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "components = 2\n",
    "K = 2\n",
    "\n",
    "# Initialise the PCA algorithm\n",
    "pca = PCA(n_components=components)\n",
    "\n",
    "# Fit and transform the PCA on the dataset\n",
    "pca.fit(X_train_scaled) \n",
    "X_train_pca = pca.transform(X_train_scaled)\n",
    "\n",
    "# Initialise the K-means algorithm\n",
    "kmeans = KMeans(n_clusters=K, random_state=seed)\n",
    "\n",
    "# Fit the on the transformed dataset\n",
    "kmeans.fit(X_train_pca) \n",
    "\n",
    "# Predict / assign the examples to a cluster\n",
    "y_pred = kmeans.predict(X_train_pca)\n",
    "\n",
    "# Generate plots\n",
    "plt.figure(figsize=(13, 5))\n",
    "plt.suptitle('%d components: explained variance = %.3f' % (components, np.sum(pca.explained_variance_ratio_)))\n",
    "\n",
    "# Predicted clusters\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.title('Predicted clusters')\n",
    "plt.scatter(X_train_pca[:,0], X_train_pca[:,1], alpha=0.8, c=y_pred, cmap=plt.cm.Paired) \n",
    "# Plot the centroids\n",
    "plt.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1], marker='X', s=500, c='black')\n",
    "\n",
    "# Real labels\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.title('Labels')\n",
    "plt.scatter(X_train_pca[:,0], X_train_pca[:,1], alpha=0.8, c=y_train, cmap=plt.cm.Paired) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to before, the generated clusters doesn't really match up with our labelled data. That's fine but let's examine how much variance is explained by the two principal components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(pca.explained_variance_ratio_)\n",
    "print(np.sum(pca.explained_variance_ratio_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The components are able to explain 52% and 27% of the variance in the dataset respectively. The first two components explain a total of 79% variance of the dataset. \n",
    "\n",
    "We can also evaluate the weights of the principal components with our features. Here, we are interested in the absolute values of the weights. Let's calculate percentages so it is easier to see which features contribute the most for each component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Total sum of each component\n",
    "component_totals = np.sum(np.abs(pca.components_), axis=1)\n",
    "\n",
    "# Calculate the percentage for each component and feature \n",
    "component_percentage = np.abs(pca.components_) / component_totals.reshape(2, 1) * 100\n",
    "components_df = pd.DataFrame(component_percentage, columns=features).round(2)\n",
    "\n",
    "components_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the magnitudes contribute the most for the first principal component. Their values range from 12.62% to 13.96% for a total of 66.26%. The second component has the most contribution from the `deVAB_r` (16.73%), `expAB_r` (16.77%) and `stellar_mass` (14.41%) features. \n",
    "\n",
    "This table is a bit difficult to read so let's create a heatmap plot based on its values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Transpose the dataframe to make it easier to read the feature names\n",
    "sb.heatmap(components_df.T, cmap='Blues')\n",
    "plt.ylabel('Feature')\n",
    "plt.xlabel('Principal component')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's much easier to identify the contribution of the features for each of the components. Now, let's find out how many components is required to explain 99% of the variance in the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initialise the algorithm\n",
    "# PCA keeps all components by default if you don't pass the argument n_components \n",
    "pca = PCA(random_state=seed)\n",
    "\n",
    "# Fit the PCA on the dataset\n",
    "pca.fit(X_train_scaled) \n",
    "\n",
    "# Calculate the cumulative sum of the explained variance of each component\n",
    "print(pca.explained_variance_ratio_.cumsum())\n",
    "\n",
    "# Plot the chart showing the explained variance ratio\n",
    "plt.ylim(0, 1.1)\n",
    "plt.xlabel('Principal component', fontsize=15)\n",
    "plt.ylabel('Explained variance', fontsize=15)\n",
    "plt.plot(pca.explained_variance_ratio_.cumsum(), '--o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are able to explain 99% of the variance in the dataset we reach the 6th principal component. \n",
    "\n",
    "Based on this result, we could transform our feature set using these 6 principal components without losing much information. This is particularly useful when we are imposed with engineering requirements where we may to reduce the training / prediction time of our models and the storage of our datasets."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocessing\n",
    "Let's use PCA as a preprocessing step for training a Logistic Regression model. Recall that PCA should only be fitted on the training dataset (`X_train`). The fitted model however will be used to transform both training and testing datasets. We can also experiment with the number of principal components used for transforming the datasets. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "# Initialise the PCA algorithm\n",
    "pca = PCA()\n",
    "\n",
    "# Fit PCA model on the training dataset\n",
    "pca.fit(X_train_scaled) \n",
    "\n",
    "# Use PCA to transform the training and test datasets\n",
    "X_train_pca = pca.transform(X_train_scaled)\n",
    "X_test_pca = pca.transform(X_test_scaled)\n",
    "\n",
    "# Principal components list\n",
    "# 1 to 11\n",
    "top_k_values = range(1, len(features) + 1)\n",
    "\n",
    "train_scores = []\n",
    "test_scores = []\n",
    "\n",
    "for top_k in top_k_values:\n",
    "    # Build the feature sets based on the top_k components\n",
    "    X_train_features = X_train_pca[:, 0:top_k]\n",
    "    X_test_features = X_test_pca[:, 0:top_k]\n",
    "\n",
    "    # Initialise our model\n",
    "    logreg_pca = LogisticRegression(random_state=seed)\n",
    "\n",
    "    # Fit the model\n",
    "    logreg_pca.fit(X_train_features, y_train)\n",
    "    \n",
    "    # Generate predictions\n",
    "    y_train_pred = logreg_pca.predict(X_train_features)\n",
    "    y_test_pred = logreg_pca.predict(X_test_features)\n",
    "    \n",
    "    # Calculate the scores\n",
    "    train_scores.append(accuracy_score(y_train, y_train_pred))\n",
    "    test_scores.append(accuracy_score(y_test, y_test_pred))\n",
    "\n",
    "# Plot the chart showing the explained variance ratio\n",
    "plt.xlabel('Principal components', fontsize=15)\n",
    "plt.ylabel('Accuracy', fontsize=15)\n",
    "plt.ylim(0, 1)\n",
    "plt.plot(train_scores, '--o', label='Training accuracy')\n",
    "plt.plot(test_scores, '--o', label='Test accuracy')\n",
    "plt.legend(fontsize=13)\n",
    "plt.show()\n",
    "\n",
    "# Print results to two decimal places\n",
    "print('Train: %s' % (['%.3f' % x for x in train_scores]))\n",
    "print('Test: %s' % (['%.3f' % x for x in test_scores]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the performance of using the full feature set compare with using the PCA components? What is the number of components achieved the best performance in your experiments?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initialise the model\n",
    "# logreg = \n",
    "\n",
    "# Fit the model using the scaled training dataset\n",
    "\n",
    "# Generate predictions\n",
    "# y_train_pred = \n",
    "# y_test_pred = \n",
    "\n",
    "# Calculate scores\n",
    "print('Training accuracy: %.2f' % (accuracy_score(y_train, y_train_pred)))\n",
    "print('Test accuracy: %.2f' % (accuracy_score(y_test, y_test_pred)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. Reproducibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "Ensuring that our work is reproducible is an important step in any analysis. **As a rule, we shouldn't place much weight on a discovery that can't be reproduced**.\n",
    "\n",
    "Notebooks like this one go a long way toward making our work reproducible. Since we documented every step as we moved along, we have a written record of what we did and why we did it — both in text and code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "We hope you found this example notebook useful for your own work and learned at least one new trick by reading through it.\n",
    "\n",
    "If you've spotted any errors or would like to contribute to this notebook, please don't hestitate to get in touch. We can be reached in the following ways:\n",
    "\n",
    "* E-mail us\n",
    "\n",
    "* [Submit an issue](https://github.com/CurtinIC/https://github.com/CurtinIC/adacs-ml-workshop/issues) on GitHub\n",
    "\n",
    "* Fork the [notebook repository](https://github.com/CurtinIC/https://github.com/CurtinIC/adacs-ml-workshop/), make the fix/addition yourself, then send over a pull request"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Further reading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[[ go back to the top ]](#Table-of-contents)\n",
    "\n",
    "This notebook covers a broad variety of topics but skips over many of the specifics. If you're looking to dive deeper into a particular topic, here's some recommended reading.\n",
    "\n",
    "**Data Science**: William Chen compiled a [list of free books](http://www.wzchen.com/data-science-books/) for newcomers to Data Science, ranging from the basics of R & Python to Machine Learning to interviews and advice from prominent data scientists.\n",
    "\n",
    "**Machine Learning**: /r/MachineLearning has a useful [Wiki page](https://www.reddit.com/r/MachineLearning/wiki/index) containing links to online courses, books, data sets, etc. for Machine Learning. There's also a [curated list](https://github.com/josephmisiti/awesome-machine-learning) of Machine Learning frameworks, libraries, and software sorted by language.\n",
    "\n",
    "**scikit-learn** has a [bunch of tutorials](http://scikit-learn.org/stable/tutorial/index.html) for those looking to learn Machine Learning in Python. Andreas Mueller's [scikit-learn workshop materials](https://github.com/amueller/scipy_2015_sklearn_tutorial) are top-notch and freely available.\n",
    "\n",
    "**pandas** has [several tutorials](http://pandas.pydata.org/pandas-docs/stable/tutorials.html) covering its myriad features.\n",
    "\n",
    "**matplotlib** has many [books, videos, and tutorials](http://matplotlib.org/resources/index.html) to teach plotting in Python.\n",
    "\n",
    "**Seaborn** has a [basic tutorial](http://stanford.edu/~mwaskom/software/seaborn/tutorial.html) covering most of the statistical plotting features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}