03_outliers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression Mini Lab\n",
    "\n",
    "In this short lab we will highlight the importance of dealing with outliers prior to fitting data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate the \"data\"\n",
    "\n",
    "The values of our input are randomly sampled from an exponentiated weibull distribution, which is often used for failure analysis, but also describes long-tailed populations.\n",
    "\n",
    "The values we are trying to predict are generated by transforming the inputs using a 2nd order polynomial, and injecting some normally distributed noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gen_func(x, coefs=[400, 10, 5], noise=0.):\n",
    "    y = np.zeros(len(x))\n",
    "    for idx, c in enumerate(coefs):\n",
    "        y += c*np.power(x, idx)\n",
    "    y += scipy.stats.norm.rvs(size=len(x), random_state=2)*noise\n",
    "    return y\n",
    "\n",
    "\n",
    "n = 1000\n",
    "f_train = 0.9\n",
    "a = 1.9\n",
    "c = 0.9\n",
    "x = scipy.stats.exponweib.rvs(a, c, size=n, random_state=3)\n",
    "y = gen_func(x, noise=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate the \"outliers\"\n",
    "\n",
    "To simulate outliers in our dataset, we will inject a small number of zero-valued entries. For the purposes of this exercise, we will assume that the outliers affect only the training data, which can happen when historical logs get corrupted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_outlier = 20\n",
    "x_outlier = 4 + np.random.random(size=n_outlier)*2\n",
    "y_outlier = np.zeros(len(x_outlier))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit a regression\n",
    "\n",
    "Now, let's split the data into a training and testing dataset using a 90/10 split. We will use a linear regression with polynomial terms to predict the values and the mean-squared error to evaluate the quality of the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "# Split the data into a training and testing set\n",
    "n_train = int(f_train*n)\n",
    "x_train = x[:n_train]\n",
    "y_train = y[:n_train]\n",
    "x_test = x[n_train:]\n",
    "y_test = y[n_train:]\n",
    "\n",
    "# create a smooth set of inputs so we can plot the prediction curve\n",
    "x_grid = np.linspace(scipy.stats.exponweib.ppf(0.01, a, c),\n",
    "                 scipy.stats.exponweib.ppf(0.999, a, c), 100)\n",
    "\n",
    "fig, axs = plt.subplots(2, 2, figsize=(14,10))\n",
    "axs = np.ravel(axs)\n",
    "\n",
    "# Set the order of the polynomial fit\n",
    "n_poly=4\n",
    "\n",
    "for ax in axs:\n",
    "    ax.plot(x_train, y_train, 'bo', alpha=0.2, label='data')\n",
    "\n",
    "# Fit the data with increasing numbers of ouliters\n",
    "for idx, n_outliers in enumerate((0, n_outlier/2, n_outlier)):\n",
    "    if n_outliers > 0:\n",
    "        ax = axs[idx+1]\n",
    "        ax.plot(x_outlier[:n_outliers], y_outlier[:n_outliers], 'ro', alpha=0.5, label='outliers')\n",
    "\n",
    "    # Inject the outliers into the training data\n",
    "    x_tmp = np.concatenate([x_train, x_outlier[:n_outliers]])\n",
    "    y_tmp = np.concatenate([y_train, y_outlier[:n_outliers]])\n",
    "\n",
    "    # Transform the inputs into a n_poly order polynomial\n",
    "    poly_train = np.array([np.power(x_tmp, i) for i in range(1, n_poly)]).T\n",
    "    poly_test = np.array([np.power(x_test, i) for i in range(1, n_poly)]).T\n",
    "    poly_grid = np.array([np.power(x_grid, i) for i in range(1, n_poly)]).T\n",
    "\n",
    "    # Fit and evaluate the regression\n",
    "    clf = LinearRegression()\n",
    "    clf.fit(poly_train, y_tmp)\n",
    "    fit_r2 = clf.score(poly_train, y_tmp)\n",
    "    y_pred = clf.predict(poly_test)\n",
    "    mse_pred = mean_squared_error(y_test, y_pred)\n",
    "    \n",
    "    print '%2d outliers: Fit R^2 = %.2f, prediction MSE = %4d' % (n_outliers, fit_r2, mse_pred)\n",
    "\n",
    "    for ax in axs[idx+1:]:\n",
    "        y_pred = clf.predict(poly_grid)\n",
    "        ax.plot(x_grid, y_pred, '-', label='Fit: %d outliers' % n_outliers)\n",
    "\n",
    "for ax in axs:\n",
    "    ax.set_ylim(0, 1200)\n",
    "    ax.legend(loc='upper left', ncol=2, framealpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Open Questions\n",
    "\n",
    "* How does the presence of outliers in the training data affect the quality of the fit?\n",
    "* How could we have detected the degredation of the fit quality during the training?\n",
    "* How is the fit quality affected if we substitute invalid values with the mean, median, or remove them entirely?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}