ml-sup.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(ml-sup)=\n",
    "# Supervised Machine Learning"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, we're going to look at *supervised* machine learning, which we might just as easily call *prediction*. In maths, we are trying to find an $ \\hat{f}(\\mathbf{x}) $ such that\n",
    "\n",
    "$$\n",
    "y = \\hat{y} + \\varepsilon = \\hat{f}(\\mathbf{x}) + \\varepsilon\n",
    "$$\n",
    "\n",
    "for a set of variables $\\mathbf{x}$ and an outcome (continuous or discrete) $y$. It's also possible for this to be a multi-valued problem, eg $\\mathbf{\\hat{y}} = \\hat{\\mathbf{f}}(\\mathbf{x})$.\n",
    "\n",
    "In the introduction to this section, we saw that there are a number of reasons why we might turn to machine learning. For *supervised* machine learning, the most important application is prediction.\n",
    "\n",
    "Supervised machine learning is usually split into two types: *regression*, which covers prediction on a continuous interval, and *classification*, which is about predicting a class from a finite set of possible discrete classes.\n",
    "\n",
    "In supervised machine learning, we will talk a lot about:\n",
    "\n",
    "- in-sample data, the data used to train a model (the data on which a model learns)\n",
    "- out-of-sample data, data held back to assess the performance of a model (data on which a model has not learned but that can be used for prediction)\n",
    "\n",
    "The other key fact to know about supervised machine learning is that the typical way to assess a model is out-of-sample goodness-of-fit. Usually, this is captured by the *mean squared error*,\n",
    "\n",
    "$$\n",
    "\\text{MSE} = \\frac{1}{N} \\displaystyle\\sum_{i=1}^N (y_i - \\hat{y}_i)^2\n",
    "$$\n",
    "\n",
    "Just applying this metric on in-sample metric would lead to a lot of overfitting as machine learning models are universal function approximators, and they are very good at it. Overfitting is a problem because it gives undue confidence that a model represents reality when, in truth, it would perform poorly on new data. Imagine you had only ever seen blue aquatic species that are fish; you might conclude that *only* fish can be blue and swim. But your model of the world would perform poorly in the real world because there are mammals that are blue and swim (not least Blue Whales!)\n",
    "\n",
    "To ensure that we are not overfitting, there are typically two approaches, which can also be combined. The first is to use penalties for complexity in the model, which pushes the algorithms to find a succinct $\\hat{f}(\\mathbf{x})$. The second is to do tests of model performance on 'held out' or 'out-of-sample' data, ie data that is not used to train the model. These give a better reflection of reality. There are a whole host of ways that we can hoodwink ourselves that our model is performing better than it is, so it's important to be on guard for ways that a machine learning model might pick-up clues in ways you might not expect. There's a famous example in which an image recognition machine learning algorithm was identifying which images were from particular scanners as a proxy for which patients had cancer rather than learning to recognise the cancer itself—you can see how the errors can be serious. On the other hand, when used properly, machine learning is extremely powerful.\n",
    "\n",
    "Although we're trying to keep the theory to a minimum here, it would be remiss not to mention the bias-variance trade-off in machine learning models. Throwing darts at a dartboard helps illustrate the concepts:\n",
    "- bias is a measure of how close to the target your shots are\n",
    "- variance is how concentrated in a small area your shots are, regardless of where their average location is\n",
    "\n",
    "if you think of each prediction of a machine learning model being a 'shot' then you can see how it applies. The following equation breaks down the mean-squared error of a prediction problem on ${\\displaystyle y=f(x)+\\varepsilon }$, where we wish to minimise ${\\displaystyle (y-{\\hat {f}}(x))^{2}}$:\n",
    "\n",
    "$$\n",
    "{\\displaystyle {\\begin{aligned}{\\text{MSE}}&=f^{2}+\\sigma ^{2}-2f\\operatorname {E} [{\\hat {f}}]+\\operatorname {Var} [{\\hat {f}}]+\\operatorname {E} [{\\hat {f}}]^{2}\\\\&=(f-\\operatorname {E} [{\\hat {f}}])^{2}+\\sigma ^{2}+\\operatorname {Var} {\\big [}{\\hat {f}}{\\big ]}\\\\[5pt]&=\\operatorname {Bias} [{\\hat {f}}]^{2}+\\sigma ^{2}+\\operatorname {Var} {\\big [}{\\hat {f}}{\\big ]}\\end{aligned}}}\n",
    "$$\n",
    "\n",
    "There are three parts here; an irreducible error, $\\mathbb{E}[\\varepsilon^2] = \\sigma^2$; a bias, $\\operatorname {Bias} [{\\hat {f}}]^{2} = (f-\\mathbb{E}[\\hat{f}])^2$; and finally a variance, $\\operatorname {Var} {\\big [}{\\hat {f}}{\\big ]} = \\mathbb{E}[X^2]-\\mathbb{E}[X]^2$.\n",
    "\n",
    "Ideally, we want to just have the irreducible error. In practice, minimising the variance increases the bias and vice versa, and this is why people talk about the bias-variance trade-off. Overfitting is a great example: it decreases bias (makes predictions more accurate) but at the cost of having a more complex, varying set of guesses (higher variance). If you're hungry for more, the [wikipedia page](https://en.wikipedia.org/wiki/Bias%E2%80%93variance_tradeoff) on this is great."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Machine learning regression\n",
    "\n",
    "The package we'll be leaning on most heavily in this chapter is [**scikit-learn**](https://scikit-learn.org), which has excellent documentation. You can use **scikit-learn** to do regular regression but it's much less intuitive for economists than the tools we saw in {ref}`regression`, so we're going to skip straight onto other algorithms."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear, penalised regression\n",
    "\n",
    "Penalised regression introduces some kind of penalty to reduce the tendency to overfit, and is the first stop-off for machine learning. In the below, $w$ represents weights and $\\alpha$ a *hyper-parameter* that the user must set (we'll talk about how to set it later). The key algorithms are:\n",
    "\n",
    "- Lasso. Think of this as 'lassoing' a small number of variables by shrinking the coefficients on other variables to zero. Gives sparse models. Solves\n",
    "        $$\\min_{w} { \\frac{1}{2n_{\\text{samples}}} ||X w - y||_2 ^ 2 + \\alpha ||w||_1}$$\n",
    "- Ridge. Think of this as similar to the lasso but instead pushing coefficients more naturally to *small* (but not zero) values. Solves\n",
    "        $$\\min_{w} || X w - y||_2^2 + \\alpha ||w||_2^2$$\n",
    "- ElasticNet. Think of this as combining the strengths of Lasso and Ridge; it's essentially a combination of both. Solves\n",
    "        $$\\min_{w} { \\frac{1}{2n_{\\text{samples}}} ||X w - y||_2 ^ 2 + \\alpha \\rho ||w||_1 + \\frac{\\alpha(1-\\rho)}{2} ||w||_2 ^ 2}$$\n",
    "        where $\\rho$ controls how much like Lasso it is vs Ridge.\n",
    "\n",
    "One of the advantages of penalised regression is that it does still produce coefficients, if you happen to need those. But it's essentially a linear model, and comes with the pros and cons of those too.\n",
    "\n",
    "Let's now look at all of these in code!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The libraries we'll be using are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "random_state = 42  # We'll use this throughout to make this page reproducible\n",
    "prng = np.random.default_rng(random_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need some ground truth data. Let's generate some. We can use the handy `make_regression()` function that comes with **scikit-learn**; note that this library is abbreviated to 'sklearn' in all import statements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_regression\n",
    "\n",
    "X, y, coef = make_regression(\n",
    "    n_samples=500, n_features=6, noise=0.5, random_state=random_state, coef=True\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first object, `X`, is an array with dimensions of (`n_samples` x `n_features`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "while the second is a vector with `n_samples` entries."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, you'll mostly be working with dataframes, so we're going to pop these into a dataframe now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_df = pd.DataFrame(X)\n",
    "reg_df[\"y\"] = y\n",
    "reg_df.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also take a look at the data; let's do a quick chart of the first feature in X against y (remember there is more than one feature/regressor though!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_df.plot.scatter(x=0, y=\"y\", s=30);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, let's get onto some actual machine learning. As noted before, we need to split our sample into a test set and a training set. **scikit-learn** has a handy function for this, called `train_test_split()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "train_df, test_df = train_test_split(reg_df, random_state=random_state, test_size=0.2)\n",
    "train_df.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We used `test_size = 0.2`, which uses 20% of the data as the test set, but we could have also passed an explicit number. 20% is not untypical, depending on the overall size of your data.\n",
    "\n",
    "Let's now pop these data into our algorithms. We're going to go through step-by-step to make clear what's happening first, but later we'll show you a way to efficiently run lots of models at once. We'll start with Lasso."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import linear_model\n",
    "\n",
    "reg_lasso = linear_model.Lasso(alpha=0.1, random_state=random_state)\n",
    "reg_lasso.fit(\n",
    "    train_df.iloc[:, :-1], train_df[\"y\"]\n",
    ")  # note that all but final column is passed as X"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well it doesn't look like much happened! But, *au contraire*, the magic has been happening. `reg_lasso` is an instance of the `Lasso` class that has been fit to the data and now carries with it coefficients that we can look at."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_lasso.coef_"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that, in Python, *everything is an object* and *objects can have state*. An instance of a Lasso object can acquire the 'state' of having been fit to some data, and then it carries its fitted state with it.\n",
    "\n",
    "With these linear models, we'll be able to compare coefficients to the known solutions to our made-up data—but bear in mind that i) the mean square error is the metric of success that we're really interested in; and ii) other machine learning algorithms will not produce coefficients.\n",
    "\n",
    "Let's now run the other two models in a similar fashion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ridge\n",
    "reg_ridge = linear_model.Ridge(alpha=0.1, random_state=random_state)\n",
    "reg_ridge.fit(train_df.iloc[:, :-1], train_df[\"y\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# elastic net w/ 50:50 split between ridge and lasso\n",
    "reg_elast = linear_model.ElasticNet(\n",
    "    alpha=0.1, l1_ratio=0.7, random_state=random_state, fit_intercept=False\n",
    ")\n",
    "reg_elast.fit(train_df.iloc[:, :-1], train_df[\"y\"])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can compute MSEs for all of these. Remember that we need to do this *out-of-sample*, but we'll look at both in-sample and out-of-sample so you can see the difference. We'll pop them into a dataframe too.\n",
    "\n",
    "To get out of sample predictions, we use the `predict()` method on the fitted model, passing new data (but only X) to it. Here's what that looks like on a few rows of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_lasso.predict(test_df.iloc[:5, :-1])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the mean squared error, we use the *scitkit-learn* built-in function and plug in the true y and the predicted y, both using the *test* data so that the MSE is out-of-sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "mean_squared_error(y_true=test_df[\"y\"], y_pred=reg_lasso.predict(test_df.iloc[:, :-1]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's great, but it's only one model! In many machine learning applications, you'll be thinking about multiple models at once. But the great thing about code is that *everything* can be scaled up. So let's populate a dataframe with all of the results for both in-sample and out-of-sample tests."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the names of each of the models we're working with programmatically, we're going to use a really good Python trick: we're going to make use of a *dunder method*, which is all methods that look like this: `__methodname__`. Dunder is short for double underscore. Methods are the operations you can do on objects of classes in Python. For example, `reg_lasso` is an object of class `Lasso`. Dunder methods are methods that many classes have by default. A common one is `.__name__` which returns the name of an object. In this case, we don't have `__name__`, but we do have `__class__`, so we can use that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_lasso.__class__"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get just the end bit (the rest is too long), we'll convert this to a string and strip out the bits we don't want, like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "str(reg_lasso.__class__).strip(\">'\").split(\".\")[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} Exercise\n",
    "Why does the above need `.strip(\">'\")`?\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll put this into a function so we can re-use it easily:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_nice_name_algo(object_in):\n",
    "    return str(object_in.__class__).strip(\">'\").split(\".\")[-1]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay now we have a convenient way to name each column of MSE results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "models = [reg_lasso, reg_ridge, reg_elast]\n",
    "model_names = [get_nice_name_algo(model) for model in models]\n",
    "\n",
    "mse_df = pd.DataFrame(index=[\"out-of-sample\", \"in-sample\"], columns=model_names)\n",
    "\n",
    "for model, name in zip(models, model_names):\n",
    "    mse_df.loc[\"out-of-sample\", name] = mean_squared_error(\n",
    "        y_true=test_df[\"y\"], y_pred=model.predict(test_df.iloc[:, :-1])\n",
    "    )\n",
    "\n",
    "mse_df.head()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add in the in-sample data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for model, name in zip(models, model_names):\n",
    "    mse_df.loc[\"in-sample\", name] = mean_squared_error(\n",
    "        y_true=train_df[\"y\"], y_pred=model.predict(train_df.iloc[:, :-1])\n",
    "    )\n",
    "\n",
    "mse_df"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Typically, we expect the in-sample MSE to be lower than the MSE for the out-of-sample data, though clearly something very strange is going on with the Elastic Net—we'll return to this later but for now note that we just used fairly random settings of *hyper-parameters* (`alpha` and `l1_ratio`) but we need to be a bit more clever with setting them!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a treat, let's look at the predicted vs actual values for the three models on the test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_df = pd.DataFrame()\n",
    "\n",
    "for model, name in zip(models, model_names):\n",
    "    results_df[name] = model.predict(test_df.iloc[:, :-1])\n",
    "\n",
    "results_df[\"truth\"] = test_df[\"y\"].values\n",
    "# NB: .values important in line above as test_df has a different index\n",
    "# that will try and align with index on results_df. .values returns only the\n",
    "# numbers with no indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "for column in results_df.columns[:-1]:\n",
    "    ax.scatter(\n",
    "        results_df[\"truth\"],\n",
    "        results_df[\"truth\"] - results_df[column],\n",
    "        alpha=0.9,\n",
    "        label=column,\n",
    "    )\n",
    "\n",
    "ax.legend()\n",
    "ax.set_ylim(-4, 4)\n",
    "ax.set_xlabel(\"y value\")\n",
    "ax.set_title(r\"Error in prediction, $y - \\hat{y}$\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-linear regression\n",
    "\n",
    "There are lots of methods we could look at in this section, but we're going to stick to the absolute classics."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Support Vector Machines\n",
    "\n",
    "Support Vector Machine algorithms do a trick where they try to convert the space of data that is assumed to be non-linear into a space that is linearly separable. This image, from Wikipedia, gives the sense of this:\n",
    "\n",
    "![Illustration of SVMs](https://upload.wikimedia.org/wikipedia/commons/thumb/f/fe/Kernel_Machine.svg/2880px-Kernel_Machine.svg.png)\n",
    "\n",
    "The way that they convert the underlying space depends on what *basis kernel* you use, and there are several available including options like linear and radial (rbf), which is the default. They can be used for regression and classification. You should be aware, however, that they are not scale invariant, so we highly recommended that you scale your data (something that we'll talk more about in due course).\n",
    "\n",
    "Let's see how an SVM does on the data we already generated—rather than the default radial basis function, we'll use a linear one as we created linear data in the first place. The `C=` keyword argument controls the regularisation (inverse proportional)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVR\n",
    "\n",
    "reg_svr = SVR(kernel=\"linear\", C=10)\n",
    "reg_svr.fit(train_df.iloc[:, :-1], train_df[\"y\"])\n",
    "\n",
    "mean_squared_error(y_true=test_df[\"y\"], y_pred=reg_svr.predict(test_df.iloc[:, :-1]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Without much work, we've already gotten close to the best estimate from our previous set of linear models. SVR doesn't scale well to very large datasets, but is incredibly powerful."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gaussian Process Regression\n",
    "\n",
    "Gaussian processes are a big, complex topic that we don't have time to go into here! Their main strength as a universal function approximator is that they are able to give an estimate of their own uncertainty of prediction (unsurprisingly, they use Bayes' rule under the hood). Another way of putting this is that 'Gaussian processes know what they don't know', and that can be useful (particularly in policy). The downside is that Gaussian processes are, like Bayesian approaches generally, computationally expensive. Let's see how one does on our (mostly linear) problem.\n",
    "\n",
    "We'll use an off the shelf kernel from an example in the *scikit-learn* [documentation](https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_prior_posterior.html#dot-product-kernel) on using Gaussian Process Regression, but you should always think about what kernel might suit your problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "from sklearn.gaussian_process.kernels import ConstantKernel, DotProduct\n",
    "\n",
    "kernel = ConstantKernel(0.1, (0.01, 10.0)) * (\n",
    "    DotProduct(sigma_0=1.0, sigma_0_bounds=(0.1, 10.0)) ** 2\n",
    ")\n",
    "\n",
    "reg_gauss = GaussianProcessRegressor(random_state=random_state, kernel=kernel)\n",
    "reg_gauss.fit(train_df.iloc[:, :-1], train_df[\"y\"])\n",
    "\n",
    "mean_predictions_gpr, std_predictions_gpr = reg_gauss.predict(\n",
    "    test_df.iloc[:, :-1],\n",
    "    return_std=True,\n",
    ")\n",
    "\n",
    "mean_squared_error(y_true=test_df[\"y\"], y_pred=mean_predictions_gpr)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As there is a posterior after doing GPR, it would be remiss not to look at it. For the high-dimensional problem we've chosen, it cannot be easily visualised, so we'll create a simple 1D problem, fit it, and look at that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.gaussian_process.kernels import RBF\n",
    "\n",
    "X_gpr_1d = prng.uniform(0, 5, 7).reshape(-1, 1)\n",
    "y_gpr_1d = np.sin((X_gpr_1d[:, 0] - 2.5) ** 2)\n",
    "\n",
    "kernel_gpr_1d = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-1, 10.0))\n",
    "\n",
    "gpr_1d = GaussianProcessRegressor(random_state=random_state, kernel=kernel_gpr_1d)\n",
    "gpr_1d.fit(X_gpr_1d, y_gpr_1d)\n",
    "\n",
    "x_gpr_grid = np.linspace(0, 5, 20)\n",
    "X_gpr_grid = x_gpr_grid.reshape(-1, 1)\n",
    "\n",
    "\n",
    "y_mean_gpr, y_std_gpr = gpr_1d.predict(X_gpr_grid, return_std=True)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(X_gpr_1d, y_gpr_1d)\n",
    "ax.fill_between(\n",
    "    x_gpr_grid,\n",
    "    y_mean_gpr - y_std_gpr,\n",
    "    y_mean_gpr + y_std_gpr,\n",
    "    color=\"tab:green\",\n",
    "    alpha=0.2,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Other Regression Algorithms\n",
    "\n",
    "There are a couple of other non-linear algorithms you might think of using for high-dimensional regression. They're really much better suited to classification, but we mention them here for completeness. They are:\n",
    "\n",
    "- the K nearest neighbours regression, imported via `from sklearn.neighbors import KNeighborsRegressor`. This algorithm is what happens when good neighbours become good friends, at least to econometricians. In the case of regression, this is basically local interpolation.\n",
    "- decision tree regression, imported via `from sklearn import tree.DecisionTreeRegressor`. We'll see more of decision trees in the classification section.\n",
    "- neural nets. These are a whole class of powerful agorithms in themselves! They are available in **scikit-learn** as `sklearn.neural_network.MLPRegressor` but, for any advanced work using neural nets, you should consider using the **pytorch** package instead.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification\n",
    "\n",
    "Classification is the art of predicting *categories* rather than continuous numbers. Many of the algorithms we saw in the regression section above can also be adapted for classification problems—for example there's a support vector classifier in **scikit-learn**. You can find a full list of algorithms and whether they support classification or not [here](https://scikit-learn.org/stable/supervised_learning.html#supervised-learning).\n",
    "\n",
    "The main difference with classification problems is how the scoring is undertaken: rather than a continuous score, it relies on counting the number of correctly predicted values as a ratio. You can read more about metrics for classification [here](https://en.wikipedia.org/wiki/Evaluation_of_binary_classifiers).\n",
    "\n",
    "It's probably easiest to show a simple example. We'll use the Iris dataset and we'll be trying to predict what species a flower is based on some other features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets, svm\n",
    "from sklearn.metrics import ConfusionMatrixDisplay, classification_report\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "iris = datasets.load_iris()\n",
    "X_iris = iris.data\n",
    "y_iris = iris.target\n",
    "class_names = iris.target_names\n",
    "\n",
    "# Split the data into a training set and a test set\n",
    "X_train_i, X_test_i, y_train_i, y_test_i = train_test_split(\n",
    "    X_iris, y_iris, random_state=random_state\n",
    ")\n",
    "\n",
    "# Run classifier, using a model that is too regularized (C too low) to see\n",
    "# the impact on the results\n",
    "classifier = svm.SVC(kernel=\"linear\", C=0.0135).fit(X_train_i, y_train_i)\n",
    "disp = ConfusionMatrixDisplay.from_estimator(\n",
    "    classifier,\n",
    "    X_test_i,\n",
    "    y_test_i,\n",
    "    display_labels=class_names,\n",
    "    cmap=plt.cm.Blues,\n",
    "    normalize=\"true\",\n",
    ")\n",
    "disp.ax_.set_title(\"Confusion matrix for classification\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from this *confusion matrix* that the diagonal has the ratios of labels that were correctly predicted (true positives) and off-diagonals are labels that were mispredicted.\n",
    "\n",
    "There are a large range of potential metrics for classification problems and which will be most appropriate in your situation will depend on the context, and in particular whether you have very imbalanced classes or you care more about different types of error. You can get some of the most popular metrics from the `classification_report()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    classification_report(\n",
    "        y_test_i, classifier.predict(X_test_i), target_names=class_names\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just because it is so ubiquitous, it seems worth showing another example using a popular classification algorithm, the *decision tree classifier*. Decision trees are well-liked partly because, in advanced forms, they are very powerful but partly because, in simple forms, they're very interpretable—as we are about to see:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import tree\n",
    "\n",
    "clf = tree.DecisionTreeClassifier()\n",
    "clf = clf.fit(X_train_i, y_train_i)\n",
    "tree.plot_tree(\n",
    "    clf, class_names=list(class_names), filled=True, feature_names=iris.feature_names\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has been a *very* quick introduction to classification, but hopefully you will take away that it's not so dissimilar from regression, you must be careful with the metrics of success you use depending on your problem, and that **scikit-learn** has many powerful options for you."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "codeforecon",
   "language": "python",
   "name": "python3"
  },
  "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.10.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}