distribute_complexity_and_validation.ipynb

{
  "metadata": {
    "kernelspec": {
      "display_name": "ml1-arm64",
      "language": "python",
      "name": "ml1-arm64"
    },
    "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.9.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4,
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "hide": true
      },
      "source": [
        "# Validation and Regularization\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "hide": true
      },
      "outputs": [],
      "source": [
        "%matplotlib inline\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "import matplotlib as mpl\n",
        "import matplotlib.cm as cm\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "hide": true
      },
      "outputs": [],
      "source": [
        "def make_simple_plot():\n",
        "    fig, axes=plt.subplots(figsize=(12,5), nrows=1, ncols=2);\n",
        "    axes[0].set_ylabel(\"$y$\")\n",
        "    axes[0].set_xlabel(\"$x$\")\n",
        "    axes[1].set_xlabel(\"$x$\")\n",
        "    axes[1].set_yticklabels([])\n",
        "    axes[0].set_ylim([-2,2])\n",
        "    axes[1].set_ylim([-2,2])\n",
        "    plt.tight_layout();\n",
        "    return axes\n",
        "def make_plot():\n",
        "    fig, axes=plt.subplots(figsize=(20,8), nrows=1, ncols=2);\n",
        "    axes[0].set_ylabel(\"$p_R$\")\n",
        "    axes[0].set_xlabel(\"$x$\")\n",
        "    axes[1].set_xlabel(\"$x$\")\n",
        "    axes[1].set_yticklabels([])\n",
        "    axes[0].set_ylim([0,1])\n",
        "    axes[1].set_ylim([0,1])\n",
        "    axes[0].set_xlim([0,1])\n",
        "    axes[1].set_xlim([0,1])\n",
        "    plt.tight_layout();\n",
        "    return axes"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## PART 1: Reading in and sampling from the data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [],
      "source": [
        "df=pd.read_csv(\"data/noisypopulation.csv\")\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [],
      "source": [
        "x=df.f.values\n",
        "f=df.x.values\n",
        "y = df.y.values"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [],
      "source": [
        "df.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From 200 points on this curve, we'll make a random choice of 60 points. We do it by choosing the indexes randomly, and then using these indexes as a way of grtting the appropriate samples"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [],
      "source": [
        "indexes=np.sort(np.random.choice(x.shape[0], size=60, replace=False))\n",
        "indexes"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [],
      "source": [
        "samplex = x[indexes]\n",
        "samplef = f[indexes]\n",
        "sampley = y[indexes]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "figure_type": "m"
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(12,8))\n",
        "plt.plot(x,f, 'k-', alpha=0.6, label=\"f\");\n",
        "plt.plot(x[indexes], y[indexes], 's', alpha=0.3, ms=10, label=\"in-sample y (observed)\");\n",
        "plt.plot(x, y, '.', alpha=0.8, label=\"population y\");\n",
        "plt.xlabel('$x$');\n",
        "plt.ylabel('$y$')\n",
        "plt.legend(loc=4);"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [],
      "source": [
        "sample_df=pd.DataFrame(dict(x=x[indexes],f=f[indexes],y=y[indexes]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Part 2: Fit on training set and predict on test set\n",
        "\n",
        "We will do the split of testing and training for you in order to illustrate how this can be done.\n",
        "\n",
        "### Train-test split"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [],
      "source": [
        "from sklearn.model_selection import train_test_split\n",
        "datasize=sample_df.shape[0]\n",
        "print(datasize)\n",
        "#split dataset using the index, as we have x,f, and y that we want to split.\n",
        "itrain,itest = train_test_split(np.arange(60),train_size=0.8)\n",
        "print(itrain.shape)\n",
        "xtrain= sample_df.x[itrain].values\n",
        "ftrain = sample_df.f[itrain].values\n",
        "ytrain = sample_df.y[itrain].values\n",
        "xtest= sample_df.x[itest].values\n",
        "ftest = sample_df.f[itest].values\n",
        "ytest = sample_df.y[itest].values"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [],
      "source": [
        "sample_df.x[itrain].values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We'll need to create polynomial features, ie add 1, x, x^2 and so on."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### The `scikit-learn` interface\n",
        "\n",
        "Scikit-learn is the main python machine learning library. It consists of many learners which can learn models from data, as well as a lot of utility functions such as `train_test_split`. It can be used in python by the incantation `import sklearn`.\n",
        "\n",
        "The library has a very well defined interface. This makes the library a joy to use, and surely contributes to its popularity. As the [scikit-learn API paper](http://arxiv.org/pdf/1309.0238v1.pdf) [Buitinck, Lars, et al. \"API design for machine learning software: experiences from the scikit-learn project.\" arXiv preprint arXiv:1309.0238 (2013).] says:\n",
        "\n",
        ">All objects within scikit-learn share a uniform common basic API consisting of three complementary interfaces: **an estimator interface for building and \ufb01tting models, a predictor interface for making predictions and a transformer interface for converting data**. The estimator interface is at the core of the library. It de\ufb01nes instantiation mechanisms of objects and exposes a `fit` method for learning a model from training data. All supervised and unsupervised learning algorithms (e.g., for classi\ufb01cation, regression or clustering) are o\ufb00ered as objects implementing this interface. Machine learning tasks like feature extraction, feature selection or dimensionality reduction are also provided as estimators.\n",
        "\n",
        "We'll use the \"estimator\" interface here, specifically the estimator `PolynomialFeatures`. The API paper again:\n",
        "\n",
        ">Since it is common to modify or \ufb01lter data before feeding it to a learning algorithm, some estimators in the library implement a transformer interface which de\ufb01nes a transform method. It takes as input some new data X and yields as output a transformed version of X. Preprocessing, feature selection, feature extraction and dimensionality reduction algorithms are all provided as transformers within the library.\n",
        "\n",
        "To start with we have one **feature** `x` to predict `y`, what we will do is the transformation:\n",
        "\n",
        "$$ x \\rightarrow 1, x, x^2, x^3, ..., x^d $$\n",
        "\n",
        "for some power $d$. Our job then is to **fit** for the coefficients of these features in the polynomial\n",
        "\n",
        "$$ a_0 + a_1 x + a_2 x^2 + ... + a_d x^d. $$\n",
        "\n",
        "In other words, we have transformed a function of one feature, into a (rather simple) **linear** function of many features. To do this we first construct the estimator as `PolynomialFeatures(d)`, and then transform these features into a d-dimensional space using the method `fit_transform`.\n",
        "\n",
        "![fit_transform](images/sklearntrans.jpg)\n",
        "\n",
        "Here is an example. The reason for using `[[1],[2],[3]]` as opposed to `[1,2,3]` is that scikit-learn expects data to be stored in a two-dimensional array or matrix with size `[n_samples, n_features]`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [],
      "source": [
        "np.array([1,2,3]).reshape(-1,1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To transform `[1,2,3]` into [[1],[2],[3]] we need to do a reshape.\n",
        "\n",
        "![reshape](images/reshape.jpg)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [],
      "source": [
        "xtest.reshape(-1,1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {},
      "outputs": [],
      "source": [
        "from sklearn.preprocessing import PolynomialFeatures\n",
        "from sklearn.linear_model import LinearRegression\n",
        "from sklearn.metrics import mean_squared_error"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {},
      "outputs": [],
      "source": [
        "PolynomialFeatures(3).fit_transform(xtest)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {},
      "outputs": [],
      "source": [
        "PolynomialFeatures(3).fit_transform(xtest.reshape(-1,1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Creating Polynomial features\n",
        "\n",
        "We'll write a function to encapsulate what we learnt about creating the polynomial features."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {},
      "outputs": [],
      "source": [
        "def make_features(train_set, test_set, degrees):\n",
        "    train_dict = {}\n",
        "    test_dict = {}\n",
        "    for d in degrees:\n",
        "        traintestdict={}\n",
        "        train_dict[d] = PolynomialFeatures(d).fit_transform(train_set.reshape(-1,1))\n",
        "        test_dict[d] = PolynomialFeatures(d).fit_transform(test_set.reshape(-1,1))\n",
        "    return train_dict, test_dict"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Doing the fit\n",
        "\n",
        "We first create our features, and some arrays to store the errors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {},
      "outputs": [],
      "source": [
        "degrees=range(21)\n",
        "train_dict, test_dict = make_features(xtrain, xtest, degrees)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {},
      "outputs": [],
      "source": [
        "error_train=np.empty(len(degrees))\n",
        "error_test=np.empty(len(degrees))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What is the fitting process? We first loop over all the **hypothesis set**s that we wish to consider: in our case this is a loop over the complexity parameter $d$, the degree of the polynomials we will try and fit. That is we start with ${\\cal H_0}$, the set of all 0th order polynomials, then do ${\\cal H_1}$, then ${\\cal H_2}$, and so on... We use the notation ${\\cal H}$ to indicate a hypothesis set. Then for each degree $d$, we obtain a best fit model. We then \"test\" this model by predicting on the test chunk, obtaining the test set error for the best-fit polynomial coefficients and for degree $d$. We move on to the next degree $d$ and repeat the process, just like before. We compare all the test set errors, and pick the degree $d_*$ and the model in ${\\cal H_{d_*}}$ which minimizes this test set error.\n",
        "\n",
        ">**YOUR TURN HERE**: For each degree d, train on the training set and predict on the test set. Store the training MSE in `error_train` and test MSE in `error_test`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {},
      "outputs": [],
      "source": [
        "#for each degree, we now fit on the training set and predict on the test set\n",
        "#we accumulate the MSE on both sets in error_train and error_test\n",
        "for d in degrees:#for increasing polynomial degrees 0,1,2...\n",
        "    Xtrain = train_dict[d]\n",
        "    Xtest = test_dict[d]\n",
        "    #set up model\n",
        "    est = LinearRegression()\n",
        "    #fit\n",
        "    est.fit(Xtrain, ytrain)\n",
        "    #predict\n",
        "    #your code here\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can find the best degree thus:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {},
      "outputs": [],
      "source": [
        "bestd = np.argmin(error_test)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {},
      "outputs": [],
      "source": [
        "plt.figure(figsize=(12,8))\n",
        "plt.plot(degrees, error_train, marker='o', label='train (in-sample)')\n",
        "plt.plot(degrees, error_test, marker='o', label='test')\n",
        "plt.axvline(bestd, 0,0.5, color='r', label=\"min test error at d=%d\"%bestd, alpha=0.3)\n",
        "plt.ylabel('mean squared error')\n",
        "plt.xlabel('degree')\n",
        "plt.legend(loc='upper left')\n",
        "plt.yscale(\"log\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "![m:caption](images/complexity-error-plot.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Validation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What we have done in picking a given $d$ as the best hypothesis is that we have used the test set as a training set. \n",
        "If we choose the best $d$ based on minimizing the test set error, we have then \"fit for\" hyperparameter $d$ on the test set. \n",
        "\n",
        "In this case, the test-set error will underestimate the true out of sample error. Furthermore, we have **contaminated the test set** by fitting for $d$ on it; it is no longer a true test set.\n",
        "\n",
        "Thus, we introduce a new **validation set** on which the complexity parameter $d$ is fit, and leave out a test set which we can use to estimate the true out-of-sample performance of our learner. The place of this set in the scheme of things is shown below:\n",
        "\n",
        "![m:caption](images/train-validate-test.png)\n",
        "\n",
        "We have split the old training set into a **new smaller training set** and a **validation set**, holding the old test aside for FINAL testing AFTER we have \"fit\" for complexity $d$. Obviously we have decreased the size of the data available for training further, but this is a price we must pay for obtaining a good estimate of the out-of-sample risk $\\cal{E}_{out}$ (also denoted as risk $R_{out}$) through the test risk $\\cal{E}_{test}$ ($R_{test}$).\n",
        "\n",
        "![m:caption](images/train-validate-test-cont.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The validation process is illustrated in these two figures. We first loop over the complexity parameter $d$, the degree of the polynomials we will try and fit. Then for each degree $d$, we obtain a best fit model $g^-_d$ where the \"minus\" superscript indicates that we fit our model on the new training set which is obtained by removing (\"minusing\") a validation chunk (often the same size as the test chunk) from the old training set. We then \"test\" this model on the validation chunk, obtaining the validation error for the best-fit polynomial coefficients and for degree $d$. We move on to the next degree $d$ and repeat the process, just like before. We compare all the validation set errors, just like we did with the test errors earlier, and pick the degree $d_*$ which minimizes this validation set error.\n",
        "\n",
        "![caption](images/train-validate-test3.png)\n",
        "\n",
        "Having picked the hyperparameter $d_*$, we retrain using the hypothesis set $\\cal{H}_{*}$ on the entire old training-set to find the parameters of the polynomial of order $d_*$ and the corresponding best fit hypothesis $g_*$. Note that we left the minus off the $g$ to indicate that it was trained on the entire old traing set. We now compute the test error on the test set as an estimate of the test risk $\\cal{E}_{test}$.\n",
        "\n",
        "Thus the **validation** set if the set on which the hyperparameter is fit. This method of splitting the data $\\cal{D}$ is called the **train-validate-test** split."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Fit on training and predict on validation\n",
        "\n",
        "\n",
        "We carry out this process for one training/validation split below. Note the smaller size of the new training set. We hold the test set at the same size."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "collapsed": true,
        "jupyter": {
          "outputs_hidden": true
        }
      },
      "outputs": [],
      "source": [
        "#we split the training set down further\n",
        "intrain,invalid = train_test_split(itrain,train_size=36, test_size=12)\n",
        "# why not just use xtrain, xtest here? its the indices we need. How could you use xtrain and xtest?\n",
        "xntrain= sample_df.x[intrain].values\n",
        "fntrain = sample_df.f[intrain].values\n",
        "yntrain = sample_df.y[intrain].values\n",
        "xnvalid= sample_df.x[invalid].values\n",
        "fnvalid = sample_df.f[invalid].values\n",
        "ynvalid = sample_df.y[invalid].values\n",
        "\n",
        "degrees=range(21)\n",
        "train_dict, valid_dict = make_features(xntrain, xnvalid, degrees)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        ">YOUR TURN HERE: Train on the smaller training set. Fit for d on the validation set.  Store the respective MSEs in `error_train` and `error_valid`. Then retrain on the entire training set using this d. Label the test set MSE with the variable `err`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "collapsed": true,
        "jupyter": {
          "outputs_hidden": true
        }
      },
      "outputs": [],
      "source": [
        "error_train=np.empty(len(degrees))\n",
        "error_valid=np.empty(len(degrees))\n",
        "#for each degree, we now fit on the smaller training set and predict on the validation set\n",
        "#we accumulate the MSE on both sets in error_train and error_valid\n",
        "#we then find the degree of polynomial that minimizes the MSE on the validation set.\n",
        "#your code here\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {},
      "outputs": [],
      "source": [
        "#calculate the degree at which validation error is minimized\n",
        "mindeg = np.argmin(error_valid)\n",
        "mindeg"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "collapsed": true,
        "jupyter": {
          "outputs_hidden": true
        }
      },
      "outputs": [],
      "source": [
        "#fit on WHOLE training set now. \n",
        "##you will need to remake polynomial features on the whole training set\n",
        "#Put MSE on the test set in variable err.\n",
        "#your code here\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We plot the training error and validation error against the degree of the polynomial, and show the test set error at the $d$ which minimizes the validation set error."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {},
      "outputs": [],
      "source": [
        "plt.plot(degrees, error_train, marker='o', label='train (in-sample)')\n",
        "plt.plot(degrees, error_valid, marker='o', label='validation')\n",
        "plt.plot([mindeg], [err], marker='s', markersize=10, label='test', alpha=0.5, color='r')\n",
        "plt.ylabel('mean squared error')\n",
        "plt.xlabel('degree')\n",
        "plt.legend(loc='upper left')\n",
        "plt.yscale(\"log\")\n",
        "print(mindeg)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "> YOUR TURN HERE: Run the set of cells for the validation process again and again. What do you see? The validation error minimizing polynomial degree might change! What happened?\n"
      ]
    }
  ]
}