Lab 7 - Cross-Validation in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lab on Cross-Validation is a python adaptation of p. 190-194 of \"Introduction to Statistical Learning\n",
    "with Applications in R\" by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani. Written\n",
    "by R. Jordan Crouser at Smith College for SDS293: Machine Learning (Spring 2016).\n",
    "\n",
    "# 5.3.1 The Validation Set Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.api as sm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we'll explore the use of the validation set approach in order to estimate the\n",
    "test error rates that result from fitting various linear models on the ${\\tt Auto}$ data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'pandas.core.frame.DataFrame'>\n",
      "Int64Index: 392 entries, 0 to 396\n",
      "Data columns (total 9 columns):\n",
      "mpg             392 non-null float64\n",
      "cylinders       392 non-null int64\n",
      "displacement    392 non-null float64\n",
      "horsepower      392 non-null float64\n",
      "weight          392 non-null int64\n",
      "acceleration    392 non-null float64\n",
      "year            392 non-null int64\n",
      "origin          392 non-null int64\n",
      "name            392 non-null object\n",
      "dtypes: float64(4), int64(4), object(1)\n",
      "memory usage: 30.6+ KB\n"
     ]
    }
   ],
   "source": [
    "df1 = pd.read_csv('Auto.csv', na_values='?').dropna()\n",
    "df1.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by using the ${\\tt sample()}$ function to split the set of observations\n",
    "into two halves, by selecting a random subset of 196 observations out of\n",
    "the original 392 observations. We refer to these observations as the training\n",
    "set.\n",
    "\n",
    "We'll use the ${\\tt random\\_state}$ parameter in order to set a seed for\n",
    "${\\tt python}$’s random number generator, so that you'll obtain precisely the same results as those shown below. It is generally a good idea to set a random seed when performing an analysis such as cross-validation\n",
    "that contains an element of randomness, so that the results obtained can be reproduced precisely at a later time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train = df1.sample(196, random_state = 1)\n",
    "test = df1[~df1.isin(train)].dropna(how = 'all')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then use the ${\\tt sm.OLS.from\\_formula()}$ to fit a linear regression to predict ${\\tt mpg}$ from ${\\tt horsepower}$ using only\n",
    "the observations corresponding to the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lm = sm.OLS.from_formula('mpg~horsepower', train)\n",
    "result = lm.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use the ${\\tt predict()}$ function to estimate the response for the test\n",
    "observations, and we use some ${\\tt numpy}$ functions to caclulate the MSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23.361902892587235\n"
     ]
    }
   ],
   "source": [
    "pred = result.predict(test)\n",
    "\n",
    "MSE = np.mean(np.square(np.subtract(test[\"mpg\"], pred)))\n",
    "    \n",
    "print(MSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the estimated test MSE for the linear regression fit is 23.36. We\n",
    "can use the ${\\tt np.power()}$ function to estimate the test error for the polynomial\n",
    "and cubic regressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.252690858350192\n",
      "20.325609366115582\n"
     ]
    }
   ],
   "source": [
    "lm2 = sm.OLS.from_formula('mpg~' + '+'.join(['np.power(horsepower,' + str(i) + ')' for i in [1,2]]), train)\n",
    "print(np.mean(np.square(np.subtract(test[\"mpg\"], lm2.fit().predict(test)))))\n",
    "\n",
    "lm3 = sm.OLS.from_formula('mpg~' + '+'.join(['np.power(horsepower,' + str(i) + ')' for i in [1,2,3]]), train)\n",
    "print(np.mean(np.square(np.subtract(test[\"mpg\"], lm3.fit().predict(test)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These error rates are 20.25 and 20.33, respectively. If we choose a different\n",
    "training set instead, then we will obtain somewhat different errors on the\n",
    "validation set. We can test this out by setting a different random seed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23.214272449679587\n",
      "19.525710963117103\n",
      "19.667628097077426\n"
     ]
    }
   ],
   "source": [
    "train = df1.sample(196, random_state = 2)\n",
    "\n",
    "lm = sm.OLS.from_formula('mpg~horsepower', train)\n",
    "print(np.mean(np.square(np.subtract(test[\"mpg\"], lm.fit().predict(test)))))\n",
    "\n",
    "lm2 = sm.OLS.from_formula('mpg~' + '+'.join(['np.power(horsepower,' + str(i) + ')' for i in [1,2]]), train)\n",
    "print(np.mean(np.square(np.subtract(test[\"mpg\"], lm2.fit().predict(test)))))\n",
    "\n",
    "lm3 = sm.OLS.from_formula('mpg~' + '+'.join(['np.power(horsepower,' + str(i) + ')' for i in [1,2,3]]), train)\n",
    "print(np.mean(np.square(np.subtract(test[\"mpg\"], lm3.fit().predict(test)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this split of the observations into a training set and a validation\n",
    "set, we find that the validation set error rates for the models with linear,\n",
    "quadratic, and cubic terms are 23.21, 19.53, and 19.67, respectively.\n",
    "\n",
    "These results are consistent with our previous findings: a model that\n",
    "predicts ${\\tt mpg}$ using a quadratic function of ${\\tt horsepower}$ performs better than\n",
    "a model that involves only a linear function of ${\\tt horsepower}$, and there is\n",
    "little evidence in favor of a model that uses a cubic function of ${\\tt horsepower}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "indices are out-of-bounds",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-108-09c439889dca>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mloo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mLeaveOneOut\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mtrain_index\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtest_index\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mloo\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m     \u001b[0mdf1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtrain_index\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/jcrouser/anaconda/envs/py3k/lib/python3.5/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   1961\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mSeries\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mIndex\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1962\u001b[0m             \u001b[0;31m# either boolean or fancy integer index\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1963\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1964\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDataFrame\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1965\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_frame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/jcrouser/anaconda/envs/py3k/lib/python3.5/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m_getitem_array\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2006\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2007\u001b[0m             \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mix\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_convert_to_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2008\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtake\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconvert\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2009\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2010\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_getitem_multilevel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/jcrouser/anaconda/envs/py3k/lib/python3.5/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36mtake\u001b[0;34m(self, indices, axis, convert, is_copy)\u001b[0m\n\u001b[1;32m   1369\u001b[0m         new_data = self._data.take(indices,\n\u001b[1;32m   1370\u001b[0m                                    \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_block_manager_axis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1371\u001b[0;31m                                    convert=True, verify=True)\n\u001b[0m\u001b[1;32m   1372\u001b[0m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_constructor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__finalize__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1373\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/jcrouser/anaconda/envs/py3k/lib/python3.5/site-packages/pandas/core/internals.py\u001b[0m in \u001b[0;36mtake\u001b[0;34m(self, indexer, axis, verify, convert)\u001b[0m\n\u001b[1;32m   3617\u001b[0m         \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3618\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mconvert\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3619\u001b[0;31m             \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmaybe_convert_indices\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3620\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3621\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mverify\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/jcrouser/anaconda/envs/py3k/lib/python3.5/site-packages/pandas/core/indexing.py\u001b[0m in \u001b[0;36mmaybe_convert_indices\u001b[0;34m(indices, n)\u001b[0m\n\u001b[1;32m   1748\u001b[0m     \u001b[0mmask\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mindices\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m|\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mindices\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1749\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mmask\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0many\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1750\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mIndexError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"indices are out-of-bounds\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1751\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mindices\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1752\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: indices are out-of-bounds"
     ]
    }
   ],
   "source": [
    "from sklearn.cross_validation import LeaveOneOut\n",
    "\n",
    "loo = LeaveOneOut(10)\n",
    "for train_index, test_index in loo:\n",
    "    df1[train_index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}