Lab 5 - LDA and QDA in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lab on Logistic Regression is a Python adaptation of p. 161-163 of \"Introduction to Statistical Learning with Applications in R\" by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani. Adapted by R. Jordan Crouser at Smith College for SDS293: Machine Learning (Spring 2016)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.lda import LDA\n",
    "from sklearn.qda import QDA\n",
    "from sklearn.metrics import confusion_matrix, classification_report, precision_score\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4.6.3 Linear Discriminant Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's return to the ${\\tt Smarket}$ data from ${\\tt ISLR}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('Smarket.csv', usecols=range(1,10), index_col=0, parse_dates=True)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will perform LDA on the ${\\tt Smarket}$ data from the ${\\tt ISLR}$ package. In ${\\tt Python}$, we can fit a LDA model using the ${\\tt LDA()}$ function, which is part of the ${\\tt lda}$ module of the ${\\tt sklearn}$ library. As we did with logistic regression and KNN, we'll fit the model using only the observations before 2005, and then test the model on the data from 2005."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_train = df[:'2004'][['Lag1','Lag2']]\n",
    "y_train = df[:'2004']['Direction']\n",
    "\n",
    "X_test = df['2005':][['Lag1','Lag2']]\n",
    "y_test = df['2005':]['Direction']\n",
    "\n",
    "lda = LDA()\n",
    "model = lda.fit(X_train, y_train)\n",
    "\n",
    "print(model.priors_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LDA output indicates prior probabilities of ${\\hat{\\pi}}_1 = 0.492$ and ${\\hat{\\pi}}_2 = 0.508$; in other words,\n",
    "49.2% of the training observations correspond to days during which the\n",
    "market went down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(model.means_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above provides the group means; these are the average\n",
    "of each predictor within each class, and are used by LDA as estimates\n",
    "of $\\mu_k$. These suggest that there is a tendency for the previous 2 days’\n",
    "returns to be negative on days when the market increases, and a tendency\n",
    "for the previous days’ returns to be positive on days when the market\n",
    "declines. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(model.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients of linear discriminants output provides the linear\n",
    "combination of ${\\tt Lag1}$ and ${\\tt Lag2}$ that are used to form the LDA decision rule.\n",
    "\n",
    "If $−0.0554\\times{\\tt Lag1}−0.0443\\times{\\tt Lag2}$ is large, then the LDA classifier will\n",
    "predict a market increase, and if it is small, then the LDA classifier will\n",
    "predict a market decline. **Note**: these coefficients differ from those produced by ${\\tt R}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ${\\tt predict()}$ function returns a list of LDA’s predictions about the movement of the market on the test data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred=model.predict(X_test)\n",
    "print(np.unique(pred, return_counts=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model assigned 70 observations to the \"Down\" class, and 182 observations to the \"Up\" class. Let's check out the confusion matrix to see how this model is doing. We'll want to compare the **predicted class** (which we can find in ${\\tt pred}$) to the **true class** (found in ${\\tt y\\_test})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(confusion_matrix(pred, y_test))\n",
    "print(classification_report(y_test, pred, digits=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get the predicted _probabilities_ using the ${\\tt predict\\_proba()}$ function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred_p = model.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying a 50% threshold to the posterior probabilities allows us to recreate\n",
    "the predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.unique(pred_p[:,1]>0.5, return_counts=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the posterior probability output by the model corresponds to\n",
    "the probability that the market will **increase**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print np.stack((pred_p[10:20,1], pred[10:20])).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we wanted to use a posterior probability threshold other than 50% in\n",
    "order to make predictions, then we could easily do so. For instance, suppose\n",
    "that we wish to predict a market decrease only if we are very certain that the\n",
    "market will indeed decrease on that day—say, if the posterior probability\n",
    "is at least 90%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.unique(pred_p[:,1]>0.9, return_counts=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No days in 2005 meet that threshold! In fact, the greatest posterior probability\n",
    "of decrease in all of 2005 was 54.2%:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "max(pred_p[:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4.6.4 Quadratic Discriminant Analysis\n",
    "We will now fit a QDA model to the ${\\tt Smarket}$ data. QDA is implemented\n",
    "in ${\\tt sklearn}$ using the ${\\tt QDA()}$ function, which is part of the ${\\tt qda}$ module. The\n",
    "syntax is identical to that of ${\\tt LDA()}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "qda = QDA()\n",
    "model2 = qda.fit(X_train, y_train)\n",
    "print model2.priors_\n",
    "print model2.means_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output contains the group means. But it does not contain the coefficients\n",
    "of the linear discriminants, because the QDA classifier involves a\n",
    "_quadratic_, rather than a linear, function of the predictors. The ${\\tt predict()}$\n",
    "function works in exactly the same fashion as for LDA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pred2=model2.predict(X_test)\n",
    "print(np.unique(pred2, return_counts=True))\n",
    "print(confusion_matrix(pred2, y_test))\n",
    "print(classification_report(y_test, pred2, digits=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, the QDA predictions are accurate almost 60% of the time,\n",
    "even though the 2005 data was not used to fit the model. This level of accuracy\n",
    "is quite impressive for stock market data, which is known to be quite\n",
    "hard to model accurately. \n",
    "\n",
    "This suggests that the quadratic form assumed\n",
    "by QDA may capture the true relationship more accurately than the linear\n",
    "forms assumed by LDA and logistic regression. However, we recommend\n",
    "evaluating this method’s performance on a larger test set before betting\n",
    "that this approach will consistently beat the market!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An Application to Carseats Data\n",
    "Let's see how the ${\\tt LDA/QDA}$ approach performs on the ${\\tt Carseats}$ data set, which is\n",
    "included with ${\\tt ISLR}$. \n",
    "\n",
    "Recall: this is a simulated data set containing sales of child car seats at 400 different stores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df2 = pd.read_csv('Carseats.csv')\n",
    "df2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See if you can build a model that predicts ${\\tt ShelveLoc}$, the shelf location (Bad, Good, or Medium) of the product at each store. Don't forget to hold out some of the data for testing!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get credit for this lab, please post your answers to the following questions:\n",
    "\n",
    "- What was your approach to building the model?\n",
    "- How did your model perform?\n",
    "- Was anything easier or more challenging than you anticipated?\n",
    "\n",
    "to Piazza: https://piazza.com/class/igwiv4w3ctb6rg?cid=23"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}