Chapter_5_sec_3.1_3.4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5.3 Lab: Cross-Validation and the Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3.1 The Validation Set Approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import pandas as pd \n",
    "import math\n",
    "import random\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.graphics.regressionplots import *\n",
    "from sklearn import datasets, linear_model\n",
    "from sklearn.model_selection import KFold, cross_val_score\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.pipeline import Pipeline\n",
    "from collections import OrderedDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Auto = pd.read_csv('data/Auto.csv', header=0, na_values='?')\n",
    "Auto = Auto.dropna().reset_index(drop=True) # drop the observation with NA values and reindex the obs from 0\n",
    "print(Auto.shape)\n",
    "print(Auto.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split the data into training and record the index of train samples\n",
    "np.random.seed(1)\n",
    "train = np.random.choice(Auto.shape[0], 196, replace=False)\n",
    "select = np.in1d(range(Auto.shape[0]), train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# start to build the model\n",
    "lm = smf.ols ('mpg~horsepower', data = Auto[select]).fit()\n",
    "print(lm.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to follow the book, get prediction for all the observations in the dataset\n",
    "# here we use ~ select to exclude the result of the training samples\n",
    "preds = lm.predict(Auto)\n",
    "square_error = (Auto['mpg'] - preds)**2\n",
    "print('-------- Test error for 1st order model --------')\n",
    "print(np.mean(square_error[~select]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build a model with 2nd order of features  \n",
    "lm2 = smf.ols ('mpg~horsepower + I(horsepower ** 2.0)', data = Auto[select]).fit()\n",
    "preds = lm2.predict(Auto)\n",
    "square_error = (Auto['mpg'] - preds)**2\n",
    "print('--------Test error for 2nd order--------')\n",
    "print(square_error[~select].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build a model with 3rd order of features  \n",
    "lm3 = smf.ols ('mpg~horsepower + I(horsepower ** 2.0) + I(horsepower ** 3.0)', data = Auto[select]).fit()\n",
    "preds = lm3.predict(Auto)\n",
    "square_error = (Auto['mpg'] - preds)**2\n",
    "print('--------Test rror for 3rd order--------')\n",
    "print(np.mean(square_error[~select]))\n",
    "\n",
    "\"\"\" \n",
    "These results are consistent with our previous findings: a model that predicts mpg using a quadratic function of \n",
    "horsepower performs better than a model that involves only a linear function of horsepower, \n",
    "and there is little evidence in favor of a model that uses a cubic function of horsepower.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if we look at the summmary for 3rd order regression, \n",
    "# the coefficient of the 3rd order term is not statistically significant. \n",
    "# I will use this as Supporting evidence for the above claim. \n",
    "print(lm3.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3.2 Leave-One-Out Cross-Validation\n",
    "The LOOCV estimates only keep one sample in the validation data and use the rest of the data to train the model. This way the training model has similar dataset comparing to the model trained on entire dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OLS fit \n",
    "ols_fit = smf.ols ('mpg~horsepower', data = Auto).fit()\n",
    "print(ols_fit.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GLM fit. Compare with OLS fit, the coeffs are the same\n",
    "glm_fit = smf.glm('mpg~horsepower', data = Auto).fit()\n",
    "print(glm_fit.params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# trying CV in Python is not as easy as that in R. It will require some manual coding.\n",
    "# to use some of implemented function in Python, we use Sklearn for linear model \n",
    "# from sklearn.model_selection import KFold, cross_val_score\n",
    "# from sklearn.preprocessing import PolynomialFeatures\n",
    "# from sklearn.linear_model import LinearRegression\n",
    "# from sklearn.pipeline import Pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us re-train the model in sklearn\n",
    "x = pd.DataFrame(Auto.horsepower)\n",
    "y = Auto.mpg\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(x, y)\n",
    "print(model.intercept_)\n",
    "print(model.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loo use folds equal to # of observations. We could also choose other number of folds.\n",
    "k_fold = KFold(n_splits=x.shape[0]) \n",
    "test = cross_val_score(model, x, y, cv=k_fold,  scoring = 'neg_mean_squared_error', n_jobs=-1)\n",
    "print(np.mean(-test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for higher order polynomial fit, we use pipline tool. \n",
    "# below shows how to fit an order 1 to 20 polynomial data and show the loo results\n",
    "# this step may take a few mins\n",
    "A = OrderedDict()\n",
    "n_split = x.shape[0]\n",
    "for porder in range(1, 21, 2):\n",
    "    model = Pipeline([('poly', PolynomialFeatures(degree=porder)), ('linear', LinearRegression())])\n",
    "    k_fold = KFold(n_splits=n_split) # loo use folds equal to # of observations\n",
    "    test = cross_val_score(model, x, y, cv=k_fold,  scoring = 'neg_mean_squared_error', n_jobs=-1)\n",
    "    A[str(porder)] = np.mean(-test)\n",
    "    \n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3.3 k-Fold Cross-Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# K-fold validation is exactly same as LOO with different n_splits parameter setup. \n",
    "# the computation time is much shorter than that of LOOCV.\n",
    "np.random.seed(2)\n",
    "A = OrderedDict()\n",
    "n_split = 10\n",
    "for porder in range(1, 21, 2):\n",
    "    model = Pipeline([('poly', PolynomialFeatures(degree=porder)), ('linear', LinearRegression())])\n",
    "    k_fold = KFold(n_splits=n_split) \n",
    "    test = cross_val_score(model, x, y, cv = k_fold,  scoring = 'neg_mean_squared_error', n_jobs = -1)\n",
    "    A[str(porder)] = np.mean(-test)\n",
    "    \n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3.4 The Bootstrap\n",
    "Bootstrap means sampling with replacement. To eliminate the effect of sample size, the norm practice is to sample the same size as original dataset with replacement.\n",
    "\n",
    "Bootstrap can be used in a lot of other places, such as estimating the accuracy of a linear regression model coeffcients / Conduct non-parametric testing (permutation test) / Estimate some complicated probability "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Portfolio = pd.read_csv('data/Portfolio.csv', header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to illustrate the use of the bootstrap on this data, we must first create a function, alpha_fn(), \n",
    "# which takes as input the (X, Y) data as well as a vector indicating which observations should be used to estimate alpha.\n",
    "def alpha_fn(data, index):\n",
    "    X = data.X.iloc[index]\n",
    "    Y = data.Y.iloc[index]\n",
    "    return (np.var(Y) - np.cov(X,Y)[0,1])/(np.var(X) + np.var(Y) - 2 * np.cov(X, Y)[0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_fn(Portfolio, range(0,100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate one set of random index with 100 elements. The array has been sorted to show there are repeat elements.\n",
    "np.sort(np.random.choice(range(0, 100), size=100, replace=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# recall the previous function with a random set of input. \n",
    "alpha_fn(Portfolio, np.random.choice(range(0, 100), size=100, replace=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# since I am not aware of boot like function in python, I just defined an ad-hoc function called boot_python()\n",
    "def boot_python(data, input_fun, iteration):\n",
    "    n = Portfolio.shape[0]\n",
    "    idx = np.random.randint(0, n, (iteration, n))\n",
    "    stat = np.zeros(iteration)\n",
    "    for i in range(len(idx)):\n",
    "        stat[i] = input_fun(data, idx[i])\n",
    "    \n",
    "    return {'Mean': np.mean(stat), 'STD': np.std(stat)}\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boot_python(Portfolio, alpha_fn, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 5"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "interpreter": {
   "hash": "4548a0e672c5b3a287feee7b2962606840aa548749d1830ef724408652b0c250"
  },
  "kernelspec": {
   "display_name": "Python 2.7.16 64-bit ('base': conda)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}