class09_Learning-Curves.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Data Mining:<br>Statistical Modeling and Learning from Data\n",
    "\n",
    "## Dr. Ciro Cattuto<br>Dr. Laetitia Gauvin<br>Dr. André Panisson\n",
    "\n",
    "### Exercises - Learning Curves, Bias-Variance\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Illustration of the Bias-Variance Tradeoff with Scikit-Learn\n",
    "\n",
    "For this section, we'll work with a simple 1D regression problem.  This will help us to\n",
    "easily visualize the data and the model, and the results generalize easily to  higher-dimensional\n",
    "datasets.  We'll explore a simple **linear regression** problem.\n",
    "This can be accomplished within scikit-learn with the `sklearn.linear_model` module.\n",
    "\n",
    "We'll create a simple nonlinear function that we'd like to fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def test_func(x, err):\n",
    "    y = np.sin(x*np.pi*2)\n",
    "    if err > 0:\n",
    "        y = np.random.normal(y, err)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create a realization of this dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_data(n=40, error=0.1, random_seed=1):\n",
    "    # randomly sample the data\n",
    "    np.random.seed(random_seed)\n",
    "    X = np.random.random(n)[:, np.newaxis]\n",
    "    y = test_func(X.ravel(), error)\n",
    "    \n",
    "    return X, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = make_data(40)\n",
    "plt.scatter(X.ravel(), y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now say we want to perform a regression on this data.  Let's use the built-in linear regression function to compute a fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn import metrics\n",
    "\n",
    "X_test = np.linspace(-0.1, 1.1, 500)[:, None]\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X, y)\n",
    "y_test = model.predict(X_test)\n",
    "\n",
    "plt.scatter(X.ravel(), y)\n",
    "plt.plot(X_test.ravel(), y_test)\n",
    "print (\"mean squared error:\", metrics.mean_squared_error(model.predict(X), y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have fit a straight line to the data, but clearly this model is not a good choice.  We say that this model is **biased**, or that it **under-fits** the data.\n",
    "\n",
    "Let's try to improve this by creating a more complicated model.  We can do this by adding degrees of freedom, and computing a polynomial regression over the inputs.  Let's make this easier by creating a quick PolynomialRegression estimator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class PolynomialRegression(LinearRegression):\n",
    "    \"\"\"Simple Polynomial Regression to 1D data\"\"\"\n",
    "    def __init__(self, degree=1, **kwargs):\n",
    "        self.degree = degree\n",
    "        LinearRegression.__init__(self, **kwargs)\n",
    "        \n",
    "    def fit(self, X, y):\n",
    "        if X.shape[1] != 1:\n",
    "            raise ValueError(\"Only 1D data valid here\")\n",
    "        Xp = X ** (1 + np.arange(self.degree))\n",
    "        return LinearRegression.fit(self, Xp, y)\n",
    "        \n",
    "    def predict(self, X):\n",
    "        Xp = X ** (1 + np.arange(self.degree))\n",
    "        return LinearRegression.predict(self, Xp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use this to fit a quadratic curve to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PolynomialRegression(degree=2)\n",
    "model.fit(X, y)\n",
    "y_test = model.predict(X_test)\n",
    "\n",
    "plt.scatter(X.ravel(), y)\n",
    "plt.plot(X_test.ravel(), y_test)\n",
    "print (\"mean squared error:\", metrics.mean_squared_error(model.predict(X), y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This reduces the mean squared error, and makes a much better fit.  What happens if we use an even higher-degree polynomial?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = PolynomialRegression(degree=60)\n",
    "model.fit(X, y)\n",
    "y_test = model.predict(X_test)\n",
    "\n",
    "plt.scatter(X.ravel(), y)\n",
    "plt.plot(X_test.ravel(), y_test)\n",
    "plt.ylim(-4, 14)\n",
    "print (\"mean squared error:\", metrics.mean_squared_error(model.predict(X), y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we increase the degree to this extent, it's clear that the resulting fit is no longer reflecting the true underlying distribution, but is more sensitive to the noise in the training data. For this reason, we call it a **high-variance model**, and we say that it **over-fits** the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detecting over-fitting\n",
    "Clearly, computing the error on the training data is not enough (we saw this previously).  But computing the *validation error* can help us determine what's going on: in particular, comparing the training error and the validation error can give you an indication of how well your data is being fit.\n",
    "\n",
    "Let's do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import cross_validation\n",
    "\n",
    "degrees = np.arange(1, 60)\n",
    "\n",
    "X, y = make_data(100)\n",
    "X_train, X_test, y_train, y_test = cross_validation.train_test_split(X, y, test_size=0.3)\n",
    "\n",
    "training_error = []\n",
    "test_error = []\n",
    "mse = metrics.mean_squared_error\n",
    "\n",
    "for d in degrees:\n",
    "    model = PolynomialRegression(d).fit(X_train, y_train)\n",
    "    training_error.append(mse(model.predict(X_train), y_train))\n",
    "    test_error.append(mse(model.predict(X_test), y_test))\n",
    "    \n",
    "# note that the test error can also be computed via cross-validation\n",
    "plt.plot(degrees, training_error, label='training')\n",
    "plt.plot(degrees, test_error, label='test')\n",
    "plt.legend()\n",
    "plt.xlabel('degree')\n",
    "plt.ylabel('MSE')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = degrees[np.argmin(test_error)]\n",
    "\n",
    "model = PolynomialRegression(d).fit(X_train, y_train)\n",
    "y_pred = model.predict(X)\n",
    "\n",
    "plt.scatter(X.ravel(), y)\n",
    "plt.scatter(X.ravel(), y_pred)\n",
    "plt.ylim(-1.5, 1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Illustration of dependency of $E_{\\mathbf{in}}$ and $E_{\\mathbf{out}}$ in relation to $N$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this section, we will see how a classification model behaves when we increase the number of samples. We will work on a 2D classification problem, and the results generalize to higher dimensional datasets.\n",
    "\n",
    "We will generate 1000 samples from a bivariate Gaussian distribution $\\mathcal{ N } ((2, 0)^T , I)$ and will label this class **RED**. Similarly, 1000 more are drawn from $ \\mathcal{N} ((0, 2)^T , I)$ and are labeled **GREEN**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(2)\n",
    "nr_samples = 2000\n",
    "\n",
    "# Generate samples from two bivariate Gaussian distributions\n",
    "samples_red   = np.random.multivariate_normal(mean=(0,2), cov=np.identity(2), size=nr_samples//2)\n",
    "samples_green = np.random.multivariate_normal(mean=(2,0), cov=np.identity(2), size=nr_samples//2)\n",
    "\n",
    "# Join the red and green datasets as X and the class definitions as y\n",
    "X = np.concatenate([samples_red, samples_green])\n",
    "y = np.zeros(nr_samples, dtype=int)\n",
    "y[nr_samples//2:] = 1\n",
    "\n",
    "# plot the red and green class points\n",
    "figure(num=None, figsize=(5, 5))\n",
    "plot(samples_red[:,0], samples_red[:,1], 'o', mec='r', mfc='none')\n",
    "plot(samples_green[:,0], samples_green[:,1], 'o', mec='g', mfc='none')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A 2D classification can be accomplished within Scikit-Learn with the `sklearn.linear_model` module.\n",
    "We will divide the dataset in two parts: a train dataset and a test dataset. This can be accomplished within Scikit-Learn with the function `train_test_split` from the package `cross_validation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn import cross_validation\n",
    "\n",
    "X_train, X_test, y_train, y_test = cross_validation.train_test_split(X, y, test_size=0.5, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's build a Logistic Regression model and evaluate both $E_{\\mathbf{in}}$ and $E_{\\mathbf{out}}$ for the instance trained with the train set, where the error is defined as the mean absolute error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import linear_model, metrics\n",
    "\n",
    "def evaluate(X_train, y_train, X_test, y_test):\n",
    "    model = linear_model.LogisticRegression()\n",
    "    model.fit(X_train, y_train)\n",
    "\n",
    "    ein = metrics.mean_absolute_error(model.predict(X_train), y_train)\n",
    "    eout = metrics.mean_absolute_error(model.predict(X_test), y_test)\n",
    "    \n",
    "    return ein, eout\n",
    "\n",
    "print (evaluate(X_train, y_train, X_test, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then let's create a function that runs the evaluation for a number of realizations, but in every realization selects randomly from the training set a fixed number of samples (`n_samples`) and evaluates the Logistic Regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_many(X_train, y_train, X_test, y_test, n_samples=100, n_realizations=100):\n",
    "    \n",
    "    n_max = X_train.shape[0]\n",
    "    ein_n = []\n",
    "    eout_n = []\n",
    "    \n",
    "    for i in range(n_realizations):\n",
    "    \n",
    "        f = np.random.choice(n_max, n_samples, replace=False)\n",
    "        X_train_filtered = X_train[f]\n",
    "        y_train_filtered = y_train[f]\n",
    "        \n",
    "        ein, eout = evaluate(X_train_filtered, y_train_filtered, X_test, y_test)\n",
    "        \n",
    "        ein_n.append(ein)\n",
    "        eout_n.append(eout)\n",
    "    \n",
    "    return np.mean(ein_n), np.mean(eout_n)\n",
    "\n",
    "print (evaluate_many(X_train, y_train, X_test, y_test, n_samples=50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will evaluate the Logistic Regression model for different datasets, fixing the number of samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random.seed(2)\n",
    "\n",
    "n_list = arange(20,1000,10)\n",
    "\n",
    "n_max = X_train.shape[0]\n",
    "\n",
    "ein_list = []\n",
    "eout_list = []\n",
    "\n",
    "for n in n_list:\n",
    "    \n",
    "    ein_n, eout_n = evaluate_many(X_train, y_train,\n",
    "                                  X_test, y_test, n_samples=n)\n",
    "    \n",
    "    ein_list.append(ein_n)\n",
    "    eout_list.append(eout_n)\n",
    "    \n",
    "plot(n_list, ein_list, label=\"$E_{\\mathbf{in}}$\")\n",
    "plot(n_list, eout_list, label=\"$E_{\\mathbf{out}}$\")\n",
    "xlabel('N')\n",
    "ylabel('Mean Absolute Error')\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise: Bias and Variance on a sine function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this exercise, we will produce data from a target function $f(x) = \\sin(\\pi x)$.\n",
    "\n",
    "For $N$ times, we will choose two points from the interval $[-1,1]$ and evaluate $f(x)$ for these points.\n",
    "Then, we will create two classes of functions: $h_0$ and $h_1$. The functions in $h_0$ are linear regressions of the type $y = b$, while $h_1$ are linear regressions of the type $y = ax + b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = lambda x: sin(pi*x)\n",
    "N = 100\n",
    "h0 = []\n",
    "h1 = []\n",
    "\n",
    "for i in  range(N):\n",
    "    d = random.uniform(-1, 1, 2) \n",
    "    y = f(d)\n",
    "    \n",
    "    g0 = poly1d(polyfit(d, y, 0))\n",
    "    h0.append(g0)\n",
    "    \n",
    "    g1 = poly1d(polyfit(d, y, 1))\n",
    "    h1.append(g1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(5,8))\n",
    "\n",
    "ax1 = fig.add_subplot(2, 1, 1)\n",
    "ax2 = fig.add_subplot(2, 1, 2)\n",
    "\n",
    "x = linspace(-1, 1, 20)\n",
    "\n",
    "ax1.plot(x, f(x), 'b')\n",
    "ylim(-2,2)\n",
    "ax2.plot(x, f(x), 'b')\n",
    "ylim(-2,2)\n",
    "\n",
    "for i in  range(N):\n",
    "    ax1.plot(x, h0[i](x), 'g', alpha=0.1)\n",
    "    ax2.plot(x, h1[i](x), 'g', alpha=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next, you will calculate the **bias** and **variance** for the two classes of functions: $h_0$ and $h_1$.\n",
    "\n",
    "For the solution, proceed with the following steps:\n",
    "\n",
    "1. Create a list `x` of 100 values between -1 and 1. These values will be used to evaluate the two classes of functions at different points.\n",
    "2. Create a list of values $g_0$ that represents the average of the functions $h_0$ in each point of `x`:\n",
    "    $${g_0}_k = \\frac{1}{N}\\sum_i^N{{h_0}_i(x_k)}$$\n",
    "3. Create a list of values $g_1$ that represents the average of the functions $h_1$ in each point of `x`:\n",
    "    $${g_1}_k = \\frac{1}{N}\\sum_i^N{{h_1}_i(x_k)}$$\n",
    "4. Assign to $\\mathbf{bias}_0$ the average of the squares of the differences between $g_0$ and $f$:\n",
    "    $$\\mathbf{bias}_0 = \\frac{1}{N}\\sum_i^N{\\left({g_0}(x_i) - f(x_i)\\right)^2}$$\n",
    "5. Assign to $\\mathbf{bias}_1$ the average of the squares of the differences between $g_1$ and $f$:\n",
    "    $$\\mathbf{bias}_1 = \\frac{1}{N}\\sum_i^N{\\left({g_1}(x_i) - f(x_i)\\right)^2}$$\n",
    "6. Assign to a matrix $\\mathbf{var}_0$ the the differences between the functions $h_0$ and $g_0$ in each point of `x`:\n",
    "    $${\\mathbf{var}_0}_{ij} = {\\left({h_0}_i(x_j) - g_0(x_j)\\right)^2}$$\n",
    "7. Assign to a matrix $\\mathbf{var}_1$ the differences between the functions $h_1$ and $g_1$ in each point of `x`:\n",
    "    $${\\mathbf{var}_1}_{ij} = {\\left({h_1}_i(x_j) - g_1(x_j)\\right)^2}$$"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}