tutorials/daniela_machine_learning/SOSTAT2021_ML_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6c62ef0b",
   "metadata": {},
   "source": [
    "# Practical Machine Learning with SDSS Data\n",
    "\n",
    "In this tutorial, we are going to use SDSS data to get some hands-on experience with machine learning. In the first exercise, we're going to train a classifier to distinguish stars from galaxies. This is an important problem, because many scientific questions depends on a clean data sets, and stars can be hard to distinguish from far-away galaxies. We don't want our study of galaxy evolution contaminated by stars in our own galaxy!\n",
    "\n",
    "As a note in advance: none of the results you'll get out of this are science-worthy. This tutorial is meant to give you a first idea for how to set up your own machine learning model. But the first, and most important lesson, is this: **don't blindly trust your ML results.** \n",
    "As with any other science project, reporting or using results from a machine learning classifier or regressor requires careful understanding of the biases and caveats, assumptions and limitations that come with the data and algorithms chosen. Because the data sets you'll be using come straight out of the SDSS catalogue, you can expect there to be funny effects (both subtle and not) that may mess up any classification you'd want to do, and in a real-world setting, this would include understanding the limitations of the instrument and the data processing, before drawing any scientific conclusions from your procedure. \n",
    "\n",
    "With that out of the way, let's have some fun with machine learning! In this tutorial, we will use python and a library called `scikit-learn` to do our machine learning, `pandas` to deal with data structures, and `matplotlib` and `seaborn` to do our plotting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d38a75cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make plots interactive and import plotting functionality\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# pretty plotting\n",
    "import seaborn as sns\n",
    "\n",
    "# my standard styles for plots\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_context(\"talk\")\n",
    "\n",
    "# Always need numpy\n",
    "import numpy as np\n",
    "\n",
    "# data array operations\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39a4dac9",
   "metadata": {},
   "source": [
    "### Load the Data\n",
    "\n",
    "Our first task is loading the data. For this exercise, we'll work with `sdss_dataset1.csv`. Your task is to find the correct file in this folder and load the data into a `pandas.DataFrame` (if you've never worked with pandas, take a look at the `read_csv` function):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b98cf7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "data =    # add your code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbc4b43b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "679ec284",
   "metadata": {},
   "source": [
    "The `head` method on your loaded `DataFrame` gives you a quick overview of what's in your data.\n",
    "\n",
    "**Exercise**: What columns do you recognize? Which ones are new to you? Which columns do you think will be su"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c96ff65",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "15206406",
   "metadata": {},
   "source": [
    "Some quick lingo: In machine learning, the things we are trying to learn are often called **labels**, and the quantities we can use to learn them are **features**. For example, in some of the data sets, you're going to try and separate stars and galaxies by their magnitudes and colours. Here, for each **sample** in your data set, you have a bunch of magnitude and colour measurements, your features, and you're trying to predict whether that sample is a galaxy or a star, its label. For the photometric redshift estimation case, you similarly have magnitudes and colours as features, and you're trying to predict the redshifts (your labels). This is called **supervised learning**. \n",
    "\n",
    "Note that in this case, we always need examples where we *know* the ground truth: we need to know the class really well, or we need to know the redshift beyond a reasonable doubt (in our case here e.g. through precise spectroscopic measurements). This is often not the case in astronomy (or, indeed, science): we often don't know exactly what our labels should be. In these cases, **unsupervised learning** can be really helpful. Some of you have data sets without labels. You'll be playing around with clustering algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcc7478b",
   "metadata": {},
   "source": [
    "## Machine Learning With Messy Data\n",
    "\n",
    "Let's start with something I told you *not* to do in the earlier class: we're just going to build a classifier and see how it does, without knowing too much of what's in the data. \n",
    "\n",
    "Normally, you wouldn't *start* by doing a classification, but for most of your data sets, there are some points we're going to make throughout this tutorial, so having a classification without knowing much about the data serves as a useful baseline. In general, though, running an ML algorithm comes at the end of *many* important steps, which is part of the point of this entire tutorial.\n",
    "\n",
    "### Splitting the Data into a Training and Test Set\n",
    "\n",
    "Our first task is to split the data into a training and a test data set. Pick the first 51% of the data and designate this as our training data set, and the remaining 49% of the DataFrame will be our test data set. This assigns a way higher fraction of data to the test data set than one normally would, but bear with me here. This entire example is a little contrived in the service of being instructional. :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3738bae2",
   "metadata": {},
   "outputs": [],
   "source": [
    "ntraining = # number of training examples\n",
    "ntest = # number of test examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96c46539",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d44057f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8c5fd5df",
   "metadata": {},
   "source": [
    "Note: the thing you're classifying on should *not* be part of the array you use to classify, so remove the `class` column from the training arrays (but make sure you store it in a separate array first! (Hint: the pandas `drop` method comes in handy here)\n",
    "\n",
    "**Advice**: It may be tempting to store the column with the label `class` in a variable `class`. Please don't! `class` is one of python's protected variable names, because it's used to set up a class. Similarly, you should not name a variable `def`: python might let you do it, but then you can never ever make a function again, because you used the keyword Python uses to create a function (i.e. when you type `def myfunction(...): ...` into a variable! There are a few of those to look out for, other examples include `list` and `lambda`. You will not believe how often I've named a list, well `list`, and could then no longer make any lists! Don't make my mistakes! :) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d60dda32",
   "metadata": {},
   "outputs": [],
   "source": [
    "classes = # store class labels here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3f10f40",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "01ffb4d0",
   "metadata": {},
   "source": [
    "**Exercise**: Now let's set up a simple classifier. For this exercise, we're going to use a *Logistic Regression* classifier, as you might have encountered during the earlier lecture.\n",
    "\n",
    "For details, you can take a look at the scikit-learn documentation for [logistic regression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html). For now, use the default parameters, then train the algorithm on the training data and explore how well it does using the test data (take a look at the `score` method of the Logistic Regression classifier). Do you think it produces good results? \n",
    "\n",
    "**Hint**: Basically all algorithms implemented in scikit-learn have the same interface. Basically all of them have a `fit` method that will fit your data, a `predict` method that will predict the classes/values of new samples, and a `score` method that tells you something about how good your algorithm is at making predictions. Some algorithms also have a `transform` and a `fit_transform` method, which allows you to transform your features (e.g. dimensionality reduction algorithms like Principal Component Analysis). Many machine learning libraries outside of scikit-learn have adopted the same structure, which is super helpful when using these libraries for algorithms that are not implemented in scikit-learn. Learning how to do a workflow in scikit-learn is well worth the investment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70b2fe24",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aaeb467",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a8c00a9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebe6c39b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "069e8dcb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4b4315cd",
   "metadata": {},
   "source": [
    "So the accuracy on the training is very low, but it seems like the model does a *perfect* job on the test data set and manages to classify *every single* training example correctly (at least in my version). Pretty cool, right? \n",
    "\n",
    "Well, no. In general, a model *cannot* do better on examples that it has never seen before than on samples it *has* seen before. A test (or validation) score that is perfect (or at least significantly better than the training score) should make us suspicious. That the optimizer complained at us doesn't help, either.\n",
    "\n",
    "**Exercise**: Take five minutes and thing about why this might happen. Take a look at the training and the test data. Are there any notable differences? \n",
    "\n",
    "**Hint**: You might want to look at the class probabilities that the Logistic Regression classifier returns. You can also plot distributions of samples and labels for the two different datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02aa2b9f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9893189b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e3c236c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a037669a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef120997",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dd80185",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0b4b58b0",
   "metadata": {},
   "source": [
    "So it looks like the model only returns probabilities of 0.5 for all of our examples! This is obviously not great; it suggests that the model knows *nothing* about what the samples in our test data are. Whether you get a prediction for a STAR or a GALAXY depends *entirely* on what the programmers of scikit-learn decided what the algorithm should return in case that the probabilities are 0.5. We might just as well have gotten back a score of 0 if they'd made a different choice. \n",
    "\n",
    "Obviously, this is bad: we *thought* that the algorithm had learned something useful, but it hadn't. But why? \n",
    "\n",
    "Well, when I generated the data, I did a silly thing: \n",
    "* I downloaded 50,000 stars and 50,000 galaxies from SDSS\n",
    "* I concatenated both together and stored them in a file\n",
    "* Then we took the first 51,000 examples from that file and called it training data, and called the rest test data.\n",
    "\n",
    "What does that do? Display the number of samples labelled \"STAR\" and those labelled \"GALAXY\" in both the training and test data set to find out. :) \n",
    "Hint: There are many ways you could do this, but casting the training and test labels into a `pandas.Series` object and using the `value_counts()` function gives you a pretty straightforward way of doing this.\n",
    "\n",
    "You can do this by hand, of course, but the scikit-learn function `train_test_split` is very useful here. Let's use that instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "beaa3609",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f61998",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d25f6ac5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f5ad43a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4307a507",
   "metadata": {},
   "source": [
    "**Exercise**: Let's train our logistic regression classifier again and see how it goes this time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f12fc67f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac4e68e5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c502c4a0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c0c1d9e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "76a59020",
   "metadata": {},
   "source": [
    "That looks a lot better, but still not great. This, in part, is because we left out a bunch of crucial steps. Remember how I said **look at the data first**? Yeah, let's do a better job of this and start from the beginning:\n",
    "\n",
    "\n",
    "\n",
    "# Machine Learning: From Start To Finish\n",
    "\n",
    "I set a really bad example above by (1) asking you to run the algorithm on the data before we'd even looked at it, and (2) using the test data set before we were ready. The idea was to give you an idea to explore some of the things that can go wrong with your machine learning if you're not careful. We're now going to do better, and run through a typical machine learning workflow from start to finish.\n",
    "\n",
    "\n",
    "## Figuring out your goal: Asking why!\n",
    "\n",
    "The first step in *any* research project, but certainly in any machine learning project, is to define your goal. What are you going to do with your results? Are you trying to learn something about physics with your data set? Are you just trying to separate out good signals from bad ones? \n",
    "\n",
    "Unfortunately, many physics questions don't necessarily easily translate to a machine learning problem (this is often true for statistics as well). Going from \"I want to know what dark matter is\" to \"run a random forest on SDSS photometric measurements\" is really hard, and requires multiple iterations of reducing your *physics* question down to one that ML or statistics (or a combination of the two) can answer. It's well worth spending significant time at this stage, because this process will give you crucial insights about every step of your analysis procedure, from which columns in your data might be useful, to which dimensionality reduction algorithms might work, to the type of appropriate ML algorithm to use. \n",
    "\n",
    "One important question I've already implicitly answered for you above: Do you have ground-truth labels in your training data? That is, do you have data for which you are reasonably sure that the labels assigned to each sample are correct? Another question I've already implicitly answered: Do you have a **classification** or a **regression** problem. In a classification problem, you have categorical labels (e.g. \"star\", \"galaxy\", ...) that you try to assign to new samples. In a regression context, the variable you try to assign is continuous (e.g. redshift). \n",
    "\n",
    "Here are some additional questions you might want to think about at the very start, but also keep in mind throughout your analysis:\n",
    "* Do you only care about the *predictions*, or do I also care about the *parameters*? That is, are you trying to learn something about the structure of the problem itself (e.g. physics), or do you just care that the predictions are right?\n",
    "* How well does your training data match the samples for which you don't know the labels? Are they from the same survey/instrument? Are there significant differences between the data you can train on, and the data you want to classify? \n",
    "* What biases do you already know if in your training data? Is your survey flux-limited? Did the team making the catalogue only care about a particular funny type of Cataclysmic Variable and leave out another you might be interested in? Is there a part of feature space that's just not covered? \n",
    "* What physical knowledge do you have about your measurements? How can that physics knowledge guide you in selecting or constructing good features? \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbc0cdec",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58ea2121",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d9cf886",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc4c1b33",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c165864",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4e756550",
   "metadata": {},
   "source": [
    "## Feature Selection and Engineering\n",
    "\n",
    "One of the crucial parts of machine learning, and the part that you will likely spend most of the time on, is selecting and engineering features for training on. Features are, essentially, meaningful summaries of your data that are ideally designed such that they make classification and regression easy. In terms of the problems you're considering today, notice that the magnitudes we've extracted from the SDSS catalogue are not actual data. They're measurements derived from the *images* that the telescope took. \n",
    "\n",
    "There are machine learning algorithms that can take raw data, e.g. the pixel values in an image. The most popular type used especially in image recognition in recent years is the Convolutional Neural Network (CNN), which takes raw data and essentially internally learns what good representations of the data are.\n",
    "\n",
    "CNNs have been enormously successful for a whole number of tasks, but it's worth pointing out that they're not always the ideal solution. They're big and very, very expensive to train (some of them can take weeks even on supercomputers!). If you have knowledge about the structure of your data and your problem, then it may be more efficient and reliable to use that knowledge to extract meaningful features. For example, we know that photons hit a CCD following a point spread function (PSF). We might not know that PSF very well, but if we do, there's no point making a neural network *learn* that a PSF exists; we can just extract magnitudes and work with them directly. So for any ML problem, it's worth thinking about what you know, and what you don't know. In cases where you *don't* know your PSF very well (or any number of things that may affect your measurements), it might be worth having your model learn that structure, but just take note that that's not universally true.\n",
    "\n",
    "Let's now turn to the data you extracted. Earlier, you've taken a look at columns and might already have some ideas about which columns are useful for the star/galaxy classification project.\n",
    "\n",
    "In any data set, you might have columns that are more useful than others! Let's find the ones we want to use for classification!\n",
    "\n",
    "**Exercise**: Print the columns in your features table. What columns are there in your table? Are all of them useful for solving your problem? Discuss with your team which ones might be useful, and which ones might not be! Reminder: you can look up what some of those columns mean [here](http://skyserver.sdss.org/dr14/en/help/browser/browser.aspx#&&history=description+SpecPhoto+V). \n",
    "\n",
    "**Exercise**: drop all columns that you don't think are useful.\n",
    "\n",
    "\n",
    "**Note**: For this exercise, I've made you columns that contain the colours as potentially useful classification features, but these might not exist in data sets that you encounter in research contexts. \n",
    "\n",
    "**Hint**: Feature engineering might include *combining* measurements into new features based on your domain knowledge of the problem you're trying to solve.\n",
    "\n",
    "**Note**: for now, let's leave the `class` column in the data, which will help us with visualization later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68daaa24",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82e55c9f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b241e6a0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bb752a3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "966c919f",
   "metadata": {},
   "source": [
    "\n",
    "### Visualizing Feature Spaces\n",
    "\n",
    "One of the most useful things you can do to star with is visualize your data! There are, of course, many different ways to visualize high-dimensional data. For example, you can make a histogram of the distributions for each feature (for classification problems, colour-coded by class works well), or you can make scatter plots or heatmaps of one feature against one another. You can also do both at the same time in what's called a *corner plot* or a *pair plot*. In python, the package [`seaborn`](https://seaborn.pydata.org/index.html) has a lot of nice visualizations pre-built to cut down the time you need to deal with nitty-gritty details of making complicated plots work.\n",
    "\n",
    "**Exercise**: Discuss and try out ways to visualize your features. Take a look at the [`pairplot`](https://seaborn.pydata.org/examples/scatterplot_matrix.html) function in seaborn. In particular, it can be useful to set `hue=\"class\"` to automatically plot separate distributions for `STAR` and `GALAXY` samples. For the classification problems, look for features that separate stars and galaxies well. For the regression, try to gauge by looking at the features how complex a model you'll likely need.\n",
    "\n",
    "**Note**: a pairplot for *all* of the features will be pretty big: pick a subset you think might be informative with respect to classification, and then do a pairplot of those. It's also helpful to start by just plotting the first 1000 samples in your DataFrame if `pairplot` runs very slowly, and to set `diag_kind = \"kde\"`\n",
    "\n",
    "**Hint**: This is a good time to scramble the order of the data, so let's do that first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "a4d879cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.sample(frac=1, replace=False).reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f9b6f32",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66494966",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94965334",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5ee7d61",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "653c0b34",
   "metadata": {},
   "source": [
    "You should see that there are values around -10000 for the colour and some of the magnitudes. This is because SDSS uses `-10000` as their placeholder for bad measurements, equivalent to the `NaN` (\"not a number\") you might see in other datasets. \n",
    "\n",
    "Let's remove all samples that have any -10000 value (there are ways to deal with missing data in machine learning, but we'll not learn about those today):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c52478ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90528ea0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1c6143e1",
   "metadata": {},
   "source": [
    "### Looking for Weird Things\n",
    "\n",
    "One of the key tasks during this stage of your analysis is to *sanity-check* your data. Are there weird things? Instrumental artifacts? Things that don't look right? This is where you explore your data and try to find (and explain) these things, and potentially remove them before training.\n",
    "\n",
    "Our star/galaxy data looks reasonably ok at this point, there are no longer any magnitudes that are wildly different from what we'd expect from either stars or galaxies in an optical survey like SDSS.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ac9c5e8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca56dcab",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f54f607f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "adb85bd4",
   "metadata": {},
   "source": [
    "### Separate Out the Classes Again + Setting up training/testing splits\n",
    "\n",
    "At this point, we should remove the 'class' column from the data again and set up our train/test splits. \n",
    "\n",
    "From now on, we'll *only* explore things on the training data, and ignore the test data until the very end.\n",
    "\n",
    "Another important thing to think about is **stratification**: if your data set is very imbalanced, randomized splitting into training and test sets can lead to all examples of the smaller class ending up in one of the two sets, but not in both. Stratified splitting uses the class labels to ensure that examples of all classes end up in the training and test sets.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "8be88027",
   "metadata": {},
   "outputs": [],
   "source": [
    "classes = data_clean[\"class\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "5b120d74",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, l_train, l_test = train_test_split(data_clean, classes, train_size=0.75, shuffle=True, stratify=None)\n",
    "X_train = X_train.drop([\"class\"], axis=1)\n",
    "X_test = X_test.drop([\"class\"], axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9a8c253",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c719d41b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cefc5f9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50f1f8eb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "91993537",
   "metadata": {},
   "source": [
    "### Dimensionality Reduction\n",
    "\n",
    "Another way to try and make sense of high-dimensional feature spaces is to perform *dimensionality reduction*. There are a lot of different ways to reduce the dimensionality of your features, and some are more useful than others, depending on the structure of your data and your problem. One idea with dimensionality reduction is to find the combination of features that gives you the most useful information in terms of your regression or classification problem. Some other times, all you want is a handy visualization in 2D, since humans in general aren't very good at thinking in higher-dimensional spaces.\n",
    "\n",
    "[Principal Component Analysis (PCA)](https://en.wikipedia.org/wiki/Principal_component_analysis) is one of the most straightforward ways of reducing dimensionality (see also [here](http://scikit-learn.org/stable/modules/classes.html#module-sklearn.decomposition) for alternative methods implemented in scikit-learn). One useful thing to try might be to project your features into a smaller space and see whether you still capture the information relevant in order to perform a good classification or regression.\n",
    "\n",
    "**Exercise**: Use PCA to project your features into 2 dimensions, and compare with the full data space. Do you think most of the information in your features is captured by the 2D representation? \n",
    "\n",
    "**Hint**: The `fit_transform` method will come in handy here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff73d6a4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f79900c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb28b8af",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee7a8d08",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b929e396",
   "metadata": {},
   "source": [
    "Another way to visualize high-dimensional data is called *t-distributed stochastic neighbour embedding*, or, for short, t-SNE. [The paper](https://lvdmaaten.github.io/publications/papers/JMLR_2008.pdf) is surprisingly readable for a computer science paper, and the method is generally pretty good at visualizing complex feature spaces (and you can spend some entertaining minutes letting your brain find fun things in the patterns it produces). \n",
    "\n",
    "One thing you have to be aware of, however, is that t-SNE *does not generalize to new examples*. Wheras methods like PCA can be trained on some data, and then the trained model applied to new samples, this is *not true* for t-SNE. So this is a method that's useful for visualization, but it doesn't necessarily produce features you want to use in your classification. Note that t-SNE is generally quite a slow algorithm, so running this on very large data sets might leave you waiting for a while!\n",
    "\n",
    "**Exercise**: Let's try it out! Visualize your feature space in 2 dimensions using t-SNE (hint: this is also implemented in scikit-learn)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b252958b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7f7816d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36d0b7f4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c1c04fb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16ac8e5d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0857d82f",
   "metadata": {},
   "source": [
    "### Scaling Features\n",
    "\n",
    "Sometimes, your features vary wildly in order of magnitude. You may have photometric magnitudes that all lie between 13 and 20, but maybe you have a feature that's between 1 and 1,000,000, so that your different dimensions have vastly different scales. Some algorithms (e.g. random forests) can deal with that pretty well, but others can't. It's worth knowing enough about the algorithms you use whether they will deal with this kind of issue or not. If the method can't, one common solution is to re-scale the features such that they all have a mean of zero and a variance of one. \n",
    "\n",
    "Scikit-learn has a few ways of scaling features and other potentially useful steps for pre-processing data. Take a look at [this tutorial](http://scikit-learn.org/stable/modules/preprocessing.html). \n",
    "\n",
    "**Exercise**: Do you think scaling will help in your problem? Try scaling your features and re-running the PCA and t-SNE algorithms. Did your results change?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cae8c41",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a33689d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a32d09e5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c95b7b6e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eabe239",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38d3e32d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0ddfd275",
   "metadata": {},
   "source": [
    "## Picking a Machine Learning Algorithm to Start With\n",
    "\n",
    "Which machine learning algorithm you choose for your problem depends, as we've discussed above, strongly on the type of problem you're trying to solve and the structure of the data. But even once you've decided whether you have a regression or classification problem, and whether it's an unsupervised or supervised one, there remain many algorithms to choose from!\n",
    "\n",
    "For classification, bs you learned earlier this week, there are two different types of algorithms to keep in mind: *generative* algorithms and *discriminative* algorithms. Generative algorithms are named this way because they can *generate* data. A popular example is Gaussian Mixture Models (GMMs). Discriminative models draw more or less complex functions in a multi-dimensional space (the number of dimensions corresponds to your number of features). Broadly, in the regression case, these methods try to find the best function to draw *through* the data points in order to model them. In the classification case, these methods try to find a surface that *separates* the different classes from another (this is also called a decision boundary). \n",
    "\n",
    "One big question is how *interpretable* the model ought to be. Simpler algorithms are often easier to understand and the results more straightforward to interpret than from, say, a random forest. A general good suggestion is to start with the simplest model you think you can get away with, and only move to more complex models if your problem demands them.\n",
    "\n",
    "Neural networks have been hugely successful in solving complex machine learning problems, in part because they can *learn* features rather than require the user to hand-craft them. While they do very well in many circumstances, be aware that this isn't *always* the case! These networks work really well (especially for image recognition) when you have no idea what your features ought to be, or there are things in your data that you can't model very easily (or that might be very expensive to model). However, if you understand the underlying structures in your data really well, then extracting physically meaningful features related to your problem can make your algorithm outperform those that need to learn these structures from scratch, or be faster to run.\n",
    "\n",
    "There is another side to this, though: if you *don't* know structures in your data, then neural networks can be very powerful *emulators* of the process you're trying to model. This has been used, for example, for modelling detectors, where running physical simulations of these detectors would be hugely expensive.\n",
    "\n",
    "**Exercise**: You've earlier tried out logistic regression as an algorithm. It's time to revisit this choice. Do you think it was a reasonable one? Given your data explorations earlier, do you think you could get away with fitting a linear model (i.e. a method linear in the parameters)? Or is your function complex enough that you'll need to fit a more complicated function? Note that for many problems, there isn't a single right answer. [This slightly tongue-in-cheek flow chart](http://scikit-learn.org/stable/tutorial/machine_learning_map/index.html) might give you a starting point which methods to explore (Note that in my science projects, I very often end up in the field labelled \"tough luck\"). \n",
    "\n",
    "Some fun methods to try:\n",
    "* Logistic Regression (linear)\n",
    "* Support Vector Machines (linear)\n",
    "* Decision Trees + Random Forests\n",
    "\n",
    "\n",
    "### Hyperparameters\n",
    "\n",
    "Before we go into actually applying a machine learning algorithm to our data, there is one more thing we need to think about: the parameters! Here's where it gets a little confusing, though, so bear with me. All maachine learning algorithms have parameters. For example, a neural network has *weights* for each of the connections between network nodes. However, basically all algorithms also have *hyperparameters*. For the K-Nearest Neighbour algorithm, one hyperparameter is the number of neighbours to use in order to determine the model value. For a random forest, one hyperparameter is the number of trees. \n",
    "\n",
    "**Exercise**: Take a look at the hyperparameters for the algorithm you picked above. Which do you think are particularly useful? Are there any that don't make sense to you at all? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41c70c4d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d01bc00",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b413c10e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "37493dce",
   "metadata": {},
   "source": [
    "### K-Fold Cross Validation\n",
    "\n",
    "Okay, so we've established that basically all machine learning methods have *hyperparameters* we somehow need to deal with.\n",
    "\n",
    "What to do with all those hyperparameters? One way to deal with them is *model selection*. Each different set of hyperparameters defines a different models, which you can compare. To compare, you'll to *score* your model in some way, i.e. determine how good it is. What \"good\" means depends on your problem at hand. Above, you used the `score` method for the K-Nearest Neighbour method, which by default uses *accuracy*, i.e. the fraction of correctly identified samples in your data set. We've also seen above that accuracy isn't always the best option, depending on what you're trying to find out.\n",
    "\n",
    "Another key piece of information is that there's a reasonably high *variance* in whatever score you compute. That is, a random slice into training and test data will produce different scores than a second random slice. One way to deal with this is called cross-validation. There are several kinds, but the most common one we'll discuss here is called **K-Fold Cross Validation**. Under this scheme, you split your training data set (*after* you've already set aside a test set!) into $K$ different slices. In the first instance, you train on nine of these slices, and test on the tenth. In the second instance, you train on slices 1-8 and 10 and test on 9. You continue this process until each of your slices has been used to test the performance. To formally distinguish these slices from your test set, these are usually called **validation** sets, because you use them to *validate* your model, and choose hyperparameters. \n",
    "\n",
    "**Note**: For problems with imbalanced class distributions, there's also a useful scikit-learn class called `StratifiedKFold`. Using this to set up your k-fold cross validation will ensure that all the subsets of data created by the k-fold cross validation always contain a proportional fraction of all classes. Otherwise, especially with classes that are rare in your data set, you might end up with training or validation sets that contain *no* examples of a class you're interested in.\n",
    "\n",
    "**Exercise**: Check the hyperparameters for the algorithm you intend to use. Normally, you should try to learn enough about your algorithm in order to make some educated guesses for what those hyperparameters should be based on your problem, your data and the algorithm itself. In this case, however, we are going to use cross-validation in order to figure out which ones you should use. So pick a parameter you think might be important to your model's performance, and run cross-validation for multiple different values of this. Print or plot the resulting scores.\n",
    "\n",
    "**Hint**: In `scikit-learn`, you can use e.g. [`cross_val_score`](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_score.html#sklearn.model_selection.cross_val_score) in order to evaluate models using cross-validation. You can learn more about evaluating the estimator performance with cross validation [here](https://scikit-learn.org/stable/modules/cross_validation.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b673464",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1ee826b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4dd8ff0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "15e9c17b",
   "metadata": {},
   "source": [
    "So if you run a 5-fold cross-validation, this function returns 4 accuracy scores for the subset of data it trained on. This means that the computer split the data set inot 5 subsets. It then trained on the first four, and calculated a validation score on the last. Then it did the same again, but left out the second-to-last subset to use as a validation set, and so on, until every subset had been used once as a validation set. You can use more or fewer folds within the cross validation, but it's useful to remember: fewer folds of cross-validation means that you will be more uncertain about your scores (it's harder to see what the variance in scores is), but with more folds, there are fewer training examples per fold, which might degrade the performance.\n",
    "\n",
    "\n",
    "**Exercise**: Write a loop that calculates cross-validated scores for a list of different values of `C`:\n",
    "\n",
    "**Hint**: You'll want explore values of `C` over multiple orders of magnitude. That's usually better done on a logarithmic scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8faba2a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e720a8ed",
   "metadata": {},
   "source": [
    "## Exploring the Results\n",
    "\n",
    "Now is the time to explore the results: where did the algorithm do well? Are there parts of parameter space that the algorithm systematically gets wrong?\n",
    "\n",
    "Here, you'll want to split your training data into a single training and single validation set using `test_train_split`, instead of doing cross-validation, so that you can explore the performance on a single validation set in more detail. \n",
    "\n",
    "Things to explore:\n",
    "* plots of different features coloured by whether the algorithm successfully classified samples or not --> are there parts of parameter space where it does badly?\n",
    "* make a [confusion matrix](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.confusion_matrix.html) (this is especially useful on problems with more than two classes :) )\n",
    "* We've not talked a lot about scoring functions in this tutorial, but in reality, you might find that accuracy, AKA \"how many samples were classified correctly\" can be a bit of a limiting performance estimator in real circumstances. Take a look at [this tutorial](https://scikit-learn.org/stable/modules/model_evaluation.html) on scoring functions for more information.\n",
    "\n",
    "Based on what you learn during this exploration, you might want to \n",
    "1. do some more hyper-parameter optimization\n",
    "2. explore a different set of features (maybe a subset, or some scientifically meaningful combination of features)\n",
    "3. try a different algorithm\n",
    "4. Celebrate your success. :)\n",
    "\n",
    "### Share Your Results!\n",
    "\n",
    "You can play around with algorithms, hyperparameters, subsets or combinations of features. When you're done, add your results to the [SOSTAT 2021 ML Exercise Leaderboard](https://board.net/p/sostat2021-ml-leaderboard). Who will get the highest accuracy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "387aec74",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Daniela's PyTorch",
   "language": "python",
   "name": "dhuppenk_pytorch"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}