Chapter_9_sec_9.6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9.6 Lab: Support Vector Machines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn.model_selection import train_test_split, GridSearchCV\n",
    "from sklearn.svm import SVC, LinearSVC\n",
    "from sklearn.metrics import accuracy_score,confusion_matrix,roc_curve, auc, classification_report\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "\n",
    "import json\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# support function to plot the decision boundary of svc and highlight the support vectors\n",
    "def plot_decision_boundary(svc, X, y, h=0.021, pad=0.21):\n",
    "    x_min, x_max = X[:, 0].min() - pad, X[:, 0].max() + pad\n",
    "    y_min, y_max = X[:, 1].min() - pad, X[:, 1].max() + pad\n",
    "    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "    Z = svc.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "    Z = Z.reshape(xx.shape)\n",
    "    plt.contourf(xx, yy, Z, cmap=plt.cm.Paired, alpha=0.2)\n",
    "    plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Paired)\n",
    "\n",
    "    # highlight the support vectors\n",
    "    sv = svc.support_vectors_\n",
    "    plt.scatter(sv[:,0], sv[:,1], c='k', marker='*', s=21, linewidths=1)\n",
    "    plt.xlim(x_min, x_max)\n",
    "    plt.ylim(y_min, y_max)\n",
    "    plt.xlabel('X1')\n",
    "    plt.ylabel('X2')\n",
    "    plt.show()\n",
    "    print('Number of support vectors: ', svc.support_.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.6.1 Support Vector Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we start from generating random dataset: following the bookm we generate a dataset with 20 observations,\n",
    "# 2 features. And we divide these into two classes.\n",
    "# set seed \n",
    "np.random.seed(21)\n",
    "X = np.random.randn(20, 2)\n",
    "y = np.repeat([-1,1], 10)\n",
    "X[y==1] = X[y==1] + 1\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], s=70, c=y, cmap=plt.cm.Paired)\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Support Vector Classifier (i.e. support vector machine with linear kernel)\n",
    "svc1 = SVC(C= 10, kernel='linear')\n",
    "svc1.fit(X, y)\n",
    "\n",
    "plot_decision_boundary(svc1, X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# as mentioned before, we could use dir() to see the methods of the class\n",
    "# I did not find a good way to print out the summary of the SVC model.\n",
    "dir(svc1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we could take a look at the defaul parameters of the SVC model\n",
    "svc1.get_params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use a small cost (c = 0.1). A smaller value of the cost parameter is being used, \n",
    "# we obtain a larger number of support vectors, because the margin is now wider. \n",
    "svc2 = SVC(C=0.1, kernel='linear')\n",
    "svc2.fit(X, y)\n",
    "\n",
    "plot_decision_boundary(svc2, X, y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we could also try to tune the cost parameter (C) of the SVC model using GridSearchCV\n",
    "# in this function, we need to specify cross validation folds and the metric to use for evaluation\n",
    "tuned_parameters = [{'C': [0.001, 0.01, 0.1, 1, 5, 10, 100]}]\n",
    "clf = GridSearchCV(SVC(kernel='linear'), tuned_parameters, cv=10, scoring='accuracy', return_train_score=True)\n",
    "clf.fit(X, y)\n",
    "clf.cv_results_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us see the best parameters. \n",
    "# This is different from the results in the book, it is very likely due to the random generation of the datasetof the data\n",
    "clf.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use the same generation process to generate test data\n",
    "X_test = np.random.randn(20, 2)\n",
    "y_test = np.repeat([-1,1], 10)\n",
    "X_test[y_test==1] = X_test[y_test==1] + 1\n",
    "\n",
    "plt.scatter(X_test[:,0], X_test[:,1], s=70, c=y_test, cmap=plt.cm.Paired)\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train a model with the optimal parameters\n",
    "svc3 = SVC(C=1, kernel='linear')\n",
    "svc3.fit(X, y)\n",
    "\n",
    "y_pred = svc3.predict(X_test)\n",
    "pd.DataFrame(confusion_matrix(y_test, y_pred),index=svc3.classes_, columns=svc3.classes_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we make our data linear separable. In the book, they add another 0.5 to seperate the data. \n",
    "# here we start from the data generation process to aviod confusion.\n",
    "np.random.seed(21)\n",
    "X = np.random.randn(20, 2)\n",
    "y = np.repeat([-1,1], 10)\n",
    "X[y==1] = X[y==1] + 2.5\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], s=70, c=y, cmap=plt.cm.Paired)\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test = np.random.randn(20, 2)\n",
    "y_test = np.repeat([-1,1], 10)\n",
    "X_test[y_test==1] = X_test[y_test==1] + 2.5\n",
    "\n",
    "plt.scatter(X_test[:,0], X_test[:,1], s=70, c=y_test, cmap=plt.cm.Paired)\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here seems the data is linear separable. We could use a bigger cost parameter (C = 100) to train the model.\n",
    "svc4 = SVC(C=100, kernel='linear')\n",
    "svc4.fit(X, y)\n",
    "\n",
    "plot_decision_boundary(svc4, X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = svc4.predict(X_test)\n",
    "pd.DataFrame(confusion_matrix(y_test, y_pred),index=svc4.classes_, columns=svc4.classes_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.6.2 Support Vector Machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generating random dataset\n",
    "np.random.seed(21)\n",
    "X = np.random.randn(200,2)\n",
    "X[:100] = X[:100] + 2\n",
    "X[101:150] = X[101:150] - 2\n",
    "y = np.concatenate([np.repeat(-1, 150), np.repeat(1,50)])\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=2)\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1], s=70, c=y, cmap=plt.cm.Paired)\n",
    "plt.xlabel('X1')\n",
    "plt.ylabel('X2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in python, we can use the same svc model abd kernel to specific the kernel\n",
    "# for rbf kernel, we need to specify the gamma parameter\n",
    "svm = SVC(C=1.0, kernel='rbf', gamma=1)\n",
    "svm.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_decision_boundary(svm, X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = svm.predict(X_test)\n",
    "pd.DataFrame(confusion_matrix(y_test, y_pred),index=svm.classes_, columns=svm.classes_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# increasing C parameter which increases more flexibility\n",
    "svm2 = SVC(C=100, kernel='rbf', gamma=1.0)\n",
    "svm2.fit(X_train, y_train)\n",
    "plot_decision_boundary(svm2, X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "The above decision boundary seems overfitting. We can compute the test accuracy of the model to\n",
    "see whether that is the case. \n",
    "\n",
    "The model (c = 1) yields a test accuracy of 0.85; the model(c = 100) yields a test accuracy of 0.77.\n",
    "\"\"\"\n",
    "y_pred = svm2.predict(X_test)\n",
    "pd.DataFrame(confusion_matrix(y_test, y_pred),index=svm2.classes_, columns=svm2.classes_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the parameters by cross-validation\n",
    "tuned_parameters = [{'C': [0.01, 0.1, 1, 10, 100],\n",
    "                     'gamma': [0.5, 1,2,3,4]}]\n",
    "clf = GridSearchCV(SVC(kernel='rbf'), tuned_parameters, cv=10, scoring='accuracy', return_train_score=True)\n",
    "clf.fit(X_train, y_train)\n",
    "clf.cv_results_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let us see the best parameters.\n",
    "clf.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# confusion matrix for the best model\n",
    "confusion_matrix(y_test, clf.best_estimator_.predict(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the test accuracy\n",
    "clf.best_estimator_.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.6.3 ROC Curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svm3 = SVC(C=1, kernel='rbf', gamma=2)\n",
    "svm3.fit(X_train, y_train)\n",
    "\n",
    "# we train another model flexible model\n",
    "svm4 = SVC(C=1, kernel='rbf', gamma=50)\n",
    "svm4.fit(X_train, y_train)\n",
    "\n",
    "y_train_score3 = svm3.decision_function(X_train)\n",
    "y_train_score4 = svm4.decision_function(X_train)\n",
    "\n",
    "false_pos_rate3, true_pos_rate3, _ = roc_curve(y_train, y_train_score3)\n",
    "roc_auc3 = auc(false_pos_rate3, true_pos_rate3)\n",
    "\n",
    "false_pos_rate4, true_pos_rate4, _ = roc_curve(y_train, y_train_score4)\n",
    "roc_auc4 = auc(false_pos_rate4, true_pos_rate4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1,ax2) = plt.subplots(1, 2, figsize=(14,6))\n",
    "ax1.plot(false_pos_rate3, true_pos_rate3, label='SVM $\\gamma = 1$ ROC curve (area = %0.2f)' % roc_auc3, color='b')\n",
    "ax1.plot(false_pos_rate4, true_pos_rate4, label='SVM $\\gamma = 50$ ROC curve (area = %0.2f)' % roc_auc4, color='r')\n",
    "ax1.set_title('Training Data')\n",
    "\n",
    "y_test_score3 = svm3.decision_function(X_test)\n",
    "y_test_score4 = svm4.decision_function(X_test)\n",
    "\n",
    "false_pos_rate3, true_pos_rate3, _ = roc_curve(y_test, y_test_score3)\n",
    "roc_auc3 = auc(false_pos_rate3, true_pos_rate3)\n",
    "\n",
    "false_pos_rate4, true_pos_rate4, _ = roc_curve(y_test, y_test_score4)\n",
    "roc_auc4 = auc(false_pos_rate4, true_pos_rate4)\n",
    "\n",
    "ax2.plot(false_pos_rate3, true_pos_rate3, label='SVM $\\gamma = 1$ ROC curve (area = %0.2f)' % roc_auc3, color='b')\n",
    "ax2.plot(false_pos_rate4, true_pos_rate4, label='SVM $\\gamma = 50$ ROC curve (area = %0.2f)' % roc_auc4, color='r')\n",
    "ax2.set_title('Test Data')\n",
    "\n",
    "for ax in fig.axes:\n",
    "    ax.plot([0, 1], [0, 1], 'k--')\n",
    "    ax.set_xlim([-0.05, 1.0])\n",
    "    ax.set_ylim([0.0, 1.05])\n",
    "    ax.set_xlabel('False Positive Rate')\n",
    "    ax.set_ylabel('True Positive Rate')\n",
    "    ax.legend(loc=\"lower right\")\n",
    "\n",
    "\"\"\" \n",
    "From the plots below, we can see that the model with gamma = 50 is overfitting the training data \n",
    "(i.e. the training metric is much better than the test metric).\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.6.4 SVM with Multiple Classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate the previously used random dataset\n",
    "np.random.seed(21)\n",
    "X = np.random.randn(200,2)\n",
    "X[:100] = X[:100] + 2\n",
    "X[101:150] = X[101:150] - 2\n",
    "y = np.concatenate([np.repeat(-1, 150), np.repeat(1,50)])\n",
    "\n",
    "# adding another class to the dataset, I used a different offset to separate the classes better\n",
    "XX = np.vstack([X, np.random.randn(50,2)])\n",
    "yy = np.hstack([y, np.repeat(0,50)])\n",
    "XX[yy==0, 1] = XX[yy==0, 1] + 6\n",
    "\n",
    "plt.scatter(XX[:,0], XX[:,1], s=70, c=yy, cmap=plt.cm.prism)\n",
    "plt.xlabel('XX1')\n",
    "plt.ylabel('XX2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the svm model \n",
    "svm5 = SVC(C=10, kernel='rbf', gamma=1)\n",
    "svm5.fit(XX, yy)\n",
    "plot_decision_boundary(svm5, XX, yy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9.6.5 Application to Gene Expression Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I saved the gene expression data as a json file, in python we could load the json file using the json library\n",
    "# after reading in the data, we can use the data is same as a dictionary, we can use the keys to access the data\n",
    "# import json\n",
    "f = open('./data/Khan.json',)\n",
    "Khan = json.load(f)\n",
    "print(Khan.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = np.array(Khan['xtrain'])\n",
    "y_train = np.array(Khan['ytrain'])\n",
    "X_test = np.array(Khan['xtest'])\n",
    "y_test = np.array(Khan['ytest'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a look at the data, we will notice there are 4 classes\n",
    "np.unique(y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svm6 = SVC(C = 10, kernel='linear')\n",
    "svm6.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "We see below that the model is perfect on training data. In fact, this is not surprising, \n",
    "because the large number of variables relative to the number of observations implies that \n",
    "it is easy to find hyperplanes that fully separate the classes. We are most interested not \n",
    "in the support vector classifier’s performance on the training observations, but rather its \n",
    "performance on the test observations.\n",
    "\"\"\"\n",
    "print('train accuracy', svm6.score(X_train, y_train))\n",
    "y_pred = svm6.predict(X_test)\n",
    "print('test accuracy', svm6.score(X_test, y_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 9"
   ]
  }
 ],
 "metadata": {
  "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": 2
}