3 - SVMs, bias-variance, and kernels/Support Vector Machines, Bias-Variance tradeoff and an Introduction to Kernel Theory.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-nc-sa/4.0/80x15.png\" /></a>&nbsp;| [Emmanuel Rachelson](https://personnel.isae-supaero.fr/emmanuel-rachelson?lang=en) | <a href=\"https://supaerodatascience.github.io/machine-learning/\">https://supaerodatascience.github.io/machine-learning/</a>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "209863ed-cab4-4427-9861-05e0e4b27167"
    }
   },
   "source": [
    "<div style=\"font-size:22pt; line-height:25pt; font-weight:bold; text-align:center;\">Support Vector Machines,<br> the Bias-Variance tradeoff<br> and a short introduction to kernel theory</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "209863ed-cab4-4427-9861-05e0e4b27167"
    }
   },
   "source": [
    "SVM is a powerfull method, used mainly for **Classification**. As we will see, we use usually a combination of two steps :\n",
    "\n",
    "- Draw a line that sits as far as possible from the data points (Support Vector Machines)\n",
    "- Send all data points in a higher dimension space where they are linearly separable (Kernel trick)\n",
    "\n",
    "\n",
    "<img id=\"fig1\" src=\"https://imgs.xkcd.com/comics/ai_methodology.png\"> \n",
    "\n",
    "1. [A geometrical approach to Machine Learning: margin maximization](#sec1)\n",
    "    1. [The notion of margin](#sec1-1)\n",
    "    2. [Optimal separating hyperplanes](#sec1-2)\n",
    "    3. [Support vectors and sparsity](#sec1-3)\n",
    "    4. [The non-linearly separable case](#sec1-4)\n",
    "    5. [Making predictions](#sec1-5)\n",
    "2. [Support Vector Machines in scikit-learn](#sec2)\n",
    "3. [When using linear separators makes no more sense](#sec3)\n",
    "4. [A word on the bias-variance tradeoff](#sec4)\n",
    "5. [The kernel trick](#sec5)\n",
    "    1. [The intuition of mapping to higher dimensions](#sec5-1)\n",
    "    2. [The kernel trick](#sec5-2)\n",
    "    3. [Positive definite kernels](#sec5-3)\n",
    "6. [SVMs and kernels](#sec6)\n",
    "7. [What about other uses?](#sec7)\n",
    "8. [Examples](#sec8)\n",
    "    1. [Spam or ham?](#sec8-1)\n",
    "    2. [NIST](#sec8-2)\n",
    "\n",
    "<div class=\"alert alert-success\"><b>In a nutshell:</b>\n",
    "<ul>\n",
    "<li> Support Vector Machines try to separate data by maximizing a geometrical margin\n",
    "<li> They are computed offline\n",
    "<li> They offer a sparse, robust to class imbalance, and easy to evaluate predictor\n",
    "<li> Kernels are a way of enriching (lifting) the data representation so that it becomes linearly separable\n",
    "<li> SVMs + kernels offer a versatile method for classification, regression and density estimation\n",
    "<li> [Documentation in scikit-learn](http://scikit-learn.org/stable/modules/svm.html)\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "34baad98-0948-40a5-a24d-3b0a38e91d67"
    }
   },
   "source": [
    "# <a id=\"sec1\"></a> 1. A geometrical approach to Machine Learning\n",
    "\n",
    "## <a id=\"sec1-1\"></a> 1.1 The notion of margin\n",
    "\n",
    "Suppose we are given the dataset $\\{(x_i,y_i)\\}$ below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "66f910f5-e213-4475-9479-db956975fa33"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "#%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "fig_size=(8, 8)\n",
    "\n",
    "res = np.loadtxt(\"sep_lin.csv\", delimiter=',')\n",
    "X = res[:,0:-1]\n",
    "y = res[:,-1].astype(int)\n",
    "Xblue = X[y==-1]\n",
    "Xred = X[y==1]\n",
    "\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b',s=20)\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r',s=20);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "ecaa380f-3794-460c-a806-eef1424a862f"
    }
   },
   "source": [
    "<div class=\"alert alert-warning\"><b>Question (collective answer)</b><br>\n",
    "<ul>\n",
    "<li> What constitutes a \"good\" separating line between red and blue dots?\n",
    "<li> What about the lines below?\n",
    "</ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "bacb254c-c598-4484-b12f-d9c3b9c3a1d0"
    }
   },
   "outputs": [],
   "source": [
    "fig_size=(10,10)\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b',s=20)\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r',s=20)\n",
    "XX = np.arange(-1.,12.,0.1)\n",
    "YY = 0.7*XX+0.4\n",
    "plt.plot(XX,YY,c='g')\n",
    "YY = 0.7*XX-0.1\n",
    "plt.plot(XX,YY,c='c')\n",
    "YY = 0.9*XX-0.2\n",
    "plt.plot(XX,YY,c='k')\n",
    "YY = 0.2*XX+4\n",
    "plt.plot(XX,YY,c='lawngreen')\n",
    "YY = 0.2*XX+8\n",
    "plt.plot(XX,YY,c='m');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best separating hyperplane is the one that maximizes the distance to the closest points.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Question**  \n",
    "Suppose $w$ is the hyperplane's normal vector, $w_0$ is its intercept, and $\\|w\\|=1$. What is the expression of this distance?\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details class=\"alert alert-danger\">\n",
    "    <summary markdown=\"span\"><b>Ready to see the answer? (click to expand)</b></summary>\n",
    "\n",
    "Let's write $M$ this distance and call it **the margin**.\n",
    "\n",
    "$$M = \\min_i y_i(w^Tx_i + w_0)$$\n",
    "</details>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "9a9f48b8-4a8a-4846-bc1e-fe4da415d5e4"
    },
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-warning\"><b>Question</b><br>\n",
    "What is the margin of the dark green line in the example above?\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "To answer this question, recall that given a hyperplane of equation $w^Tx + w_0=0$, with $\\|w\\|=1$, the (signed) distance between $x_i$ and the hyperplane is precisely $w^Tx_i + w_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "a077da5d-794e-45b6-8349-db14f8f26612"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/code1.py\n",
    "### WRITE YOUR CODE HERE\n",
    "# If you get stuck, uncomment the line above to load a correction in this cell (then you can execute this code).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "2c0d731f-7912-491d-af07-48c2187481b2"
    }
   },
   "source": [
    "## <a id=\"sec1-2\"></a> 1.2 Optimal separating hyperplanes\n",
    "\n",
    "We want to find the hyperplane that maximizes the margin, so:\n",
    "$$\\max_{M,w,w_0} M$$\n",
    "$$y_i(w^Tx_i + w_0) \\geq M, \\forall i \\in [1,N]$$\n",
    "$$\\|w\\| = 1$$\n",
    "Let's use some optimization routines for that."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice**  \n",
    "The `fmin_slsqp` function from `scipy.optimize` solves constrained non-linear optimization problems by sequential least squares programming. It takes as inputs:\n",
    "- a function `func` of variable vector `x`,\n",
    "- an initial guess `x0` for `x`,\n",
    "- a function `eq_const` that returns a vector of values corresponding to the equality constraints `eq_const(x)==0`,\n",
    "- a function `ineq_const` that returns a vector of values corresponding to the inequality constraints `ineq_const(x)>=0`.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "a263ab46-84c6-44ee-b063-a2d86d75f4c5"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/code2.py\n",
    "### WRITE YOUR CODE HERE\n",
    "# If you get stuck, uncomment the line above to load a correction in this cell (then you can execute this code).\n",
    "\n",
    "import scipy.optimize as sopt\n",
    "\n",
    "def func(x):\n",
    "    \"\"\"x = (M,w_1,w_2,b)\"\"\"\n",
    "    return -x[0]\n",
    "\n",
    "def ineq_constr(x):\n",
    "    \"\"\"x = (M,w_1,w_2,b)\"\"\"\n",
    "    return y*(x[1]*X[:,0] + x[2]*X[:,1] + x[3]*np.ones(X.shape[0])) - x[0]*np.ones(X.shape[0])\n",
    "\n",
    "def eq_const(x):\n",
    "    return x[1]*x[1]+x[2]*x[2]-1\n",
    "\n",
    "x0 = np.array([0.,0.7,-1.,0.4])\n",
    "\n",
    "\n",
    "\n",
    "res = sopt.fmin_slsqp(func, x0, f_eqcons=eq_const, f_ieqcons=ineq_constr)\n",
    "\n",
    "\n",
    "\n",
    "print(res)\n",
    "print(\"margin =\", res[0])\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b')\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r')\n",
    "XX = np.arange(-1.,12.,0.1)\n",
    "YY = -(res[1]*XX+res[3])/res[2]\n",
    "plt.plot(XX,YY,c='g')\n",
    "YY = -(res[1]*XX+res[3]+res[0])/res[2]\n",
    "plt.plot(XX,YY,'g--')\n",
    "YY = -(res[1]*XX+res[3]-res[0])/res[2]\n",
    "plt.plot(XX,YY,'g--');\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "8d050713-d7a7-44ae-93b5-d75d810aca9b"
    }
   },
   "source": [
    "We can simplify the optimization problem above.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**  \n",
    "Use the variable change $w' = \\frac{w}{M}$ and $w'_0 = \\frac{w_0}{M}$ to get rid of $M$ in the optimization problem above.\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "8d050713-d7a7-44ae-93b5-d75d810aca9b"
    }
   },
   "source": [
    "<details class=\"alert alert-danger\">\n",
    "    <summary markdown=\"span\"><b>Ready to see the answer? (click to expand)</b></summary>\n",
    "\n",
    "We have $w=M w'$ and $w_0 = M w'_0$, so the problem becomes:\n",
    "$$\\max_{M,w',w'_0} M$$\n",
    "$$y_i(M w'^Tx_i + M w'_0) \\geq M, \\forall i \\in [1,N]$$\n",
    "$$\\|M w'\\| = 1$$\n",
    "And more simply:\n",
    "$$\\max_{M,w',w'_0} M$$\n",
    "$$y_i(w'^Tx_i + w'_0) \\geq 1, \\forall i \\in [1,N]$$\n",
    "$$M = \\frac{1}{\\|w'\\|}$$\n",
    "And finally, by noting that $\\arg\\max \\frac{1}{\\|w'\\|} = \\arg\\min \\|w'\\|$:\n",
    "$$\\min_{w',w'_0} \\frac{1}{2}\\|w'\\|^2$$\n",
    "$$y_i (w'^T x_i + w'_0) \\geq 1, \\forall i \\in [1,N]$$\n",
    "Let's just rename $w'$ into $w$ and $w'_0$ into $w_0$:\n",
    "$$\\min_{w,w_0} \\frac{1}{2}\\|w\\|^2$$\n",
    "$$y_i (w^T x_i + w_0) \\geq 1, \\forall i \\in [1,N]$$\n",
    "</details>\n",
    "\n",
    "Let's check we obtain the same result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "768b9afa-8f8b-4c20-bbf3-09b3dcb501aa"
    }
   },
   "outputs": [],
   "source": [
    "def func(w):\n",
    "    \"\"\"w=(w_0,w_1,w_2)\"\"\"\n",
    "    return .5*(w[1]**2 + w[2]**2)\n",
    "\n",
    "def ineq_constr(w):\n",
    "    \"\"\"w=(w_0,w_1,w_2)\"\"\"\n",
    "    return y*(w[1]*X[:,0] + w[2]*X[:,1] + w[0]*np.ones(X.shape[0])) - np.ones(X.shape[0])\n",
    "\n",
    "x0 = np.array([0.4,0.7,-1.])\n",
    "res = sopt.fmin_slsqp(func, x0, f_ieqcons=ineq_constr)\n",
    "print(res)\n",
    "margin = 1. / np.linalg.norm(np.array([res[1],res[2]]))\n",
    "print(\"margin =\", margin)\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b')\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r')\n",
    "XX = np.arange(-1.,12.,0.1)\n",
    "YY = -(res[1]*XX+res[0])/res[2]\n",
    "plt.plot(XX,YY,c='g')\n",
    "YY = -(res[1]*XX+res[0]+1.)/res[2]\n",
    "plt.plot(XX,YY,'g--')\n",
    "YY = -(res[1]*XX+res[0]-1.)/res[2]\n",
    "plt.plot(XX,YY,'g--');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "37c49d2d-9fc6-4fed-a2b7-82bf5ea660cd"
    }
   },
   "source": [
    "## <a id=\"sec1-3\"></a> 1.3 Support Vectors and sparsity\n",
    "\n",
    "This is a quadratic programming problem. Let's try to solve it analytically.\n",
    "\n",
    "We can introduce the Lagrange multipliers $\\alpha_i$. Then the Lagrangian of the problem is:\n",
    "$$L(w,w_0,\\alpha) = \\frac{1}{2}\\|w\\|^2 - \\sum\\limits_{i=1}^N \\alpha_i\\left(y_i\\left(w^T x_i + w_0 \\right)-1\\right)$$\n",
    "\n",
    "The solution must respect the first order Karush-Kuhn-Tucker conditions:\n",
    "$$\\left\\{ \\begin{array}{l}\n",
    "\\frac{\\partial L}{\\partial w} = 0 \\Rightarrow w = \\sum\\limits_{i=1}^N \\alpha_i y_i x_i\\\\\n",
    "\\frac{\\partial L}{\\partial w_0} = 0 \\Rightarrow 0 = \\sum\\limits_{i=1}^N \\alpha_i y_i\\\\\n",
    "\\forall i=1..N, \\ \\alpha_i\\left(y_i\\left(w^T x_i + w_0 \\right)-1\\right) = 0\\\\\n",
    "\\forall i=1..N, \\alpha_i\\geq 0\n",
    "\\end{array}\\right.$$\n",
    "\n",
    "Let's remark two things:\n",
    "1. $w = \\sum\\limits_{i=1}^N \\alpha_i y_i x_i$, so $w$ is a linear combination of the $x_i$\n",
    "2. $\\alpha_i\\left(y_i\\left(w^T x_i + w_0 \\right)-1\\right) = 0$ so at least one of the elements in this product is zero.\n",
    "\n",
    "So we have two possibilities:\n",
    "- $\\alpha_i>0$, then $y_i\\left(w^T x_i + w_0 \\right)=1$; $x_i$ is on the margin's boundary\n",
    "- $\\alpha_i=0$, then $x_i$ is anywhere on the boundary or further... but does not participate in $w$.\n",
    "\n",
    "Finding the optimal $\\alpha_i$? That requires solving the dual problem (obtained by replacing $w$ in the Lagrangian above):\n",
    "$$\\max\\limits_{\\alpha\\in{\\mathbb{R}^+}^N} L_D(\\alpha) = \\sum\\limits_{i=1}^N \\alpha_i - \\frac{1}{2}\\sum\\limits_{i=1}^N\\sum\\limits_{j=1}^N \\alpha_i \\alpha_j y_i y_j x_i^T x_j\\\\\\text{such that } \\sum\\limits_{i=1}^N \\alpha_i y_i = 0$$\n",
    "\n",
    "This dual problem is usually solved by a dedicated algorithm called SMO (John Platt (1998), Sequential Minimal Optimization: A Fast Algorithm for Training Support Vector Machines)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once $\\alpha$ is found, $w=\\sum\\limits_{i=1}^N \\alpha_i y_i x_i$ follows, and computing $w_0$ is almost trivial \n",
    "\n",
    "We have seen how to compute $w$ but what about $w_0$? Take any support vector $x_i$, its $\\alpha_i$ is non-zero, so $y_i\\left(w^T x_i + w_0 \\right)=1$, so $w_0 = y_i-w^Tx_i$. In practice, to avoid numerical errors, one often averages over all support vectors with:\n",
    "$$w_0 = \\frac{\\sum\\limits_{i \\ st. \\ \\alpha_i\\neq 0} y_i-w^Tx_i}{\\sum\\limits_{i \\ st. \\ \\alpha_i\\neq 0} 1}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "37c49d2d-9fc6-4fed-a2b7-82bf5ea660cd"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**Overall:**  \n",
    "$$w = \\sum\\limits_{i=1}^N \\alpha_i y_i x_i$$\n",
    "The $x_i$ for which $\\alpha_i>0$ are called **Support Vectors**. They lay on the margin's boundary.\n",
    "</div>\n",
    "\n",
    "In the end, the optimal separating hyperplane's normal vector is defined as a linear combination of the support vectors' coordinates. This linear combination is **sparse**, it has few non-zero coefficients.\n",
    "\n",
    "That's it, we've built our first SVM!"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "37c49d2d-9fc6-4fed-a2b7-82bf5ea660cd"
    }
   },
   "source": [
    "## <a id=\"sec1-4\"></a> 1.4 The non-linearly separable case\n",
    "\n",
    "Ok, that's all very nice, but what if it is not possible to separate the data with a line, just as in the example below?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "47f522cf-664e-498b-850c-9410ab580ed6"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "from sklearn.utils import shuffle\n",
    "\n",
    "X1, y1 = datasets.make_gaussian_quantiles(cov=2.0, n_samples=300, n_features=2, n_classes=1)\n",
    "X1[:,0] = 3. + X1[:,0]\n",
    "X1[:,1] = 6. + X1[:,1]/2.5\n",
    "X2, y2 = datasets.make_gaussian_quantiles(cov=1.5, n_samples=300, n_features=2, n_classes=1)\n",
    "X2[:,0] = 8. + X2[:,0]\n",
    "X2[:,1] = 4. + X2[:,1]\n",
    "X3, y3 = datasets.make_gaussian_quantiles(cov=1.5, n_samples=300, n_features=2, n_classes=1)\n",
    "X3[:,0] = 7. + X3[:,0]\n",
    "X3[:,1] = 8. + X3[:,1]\n",
    "X = np.concatenate((X1, X2, X3))\n",
    "y = np.concatenate((y1, - y2 + 1, y3))\n",
    "y = 2*y-1\n",
    "X, y = shuffle(X, y)\n",
    "Xblue = X[y==-1]\n",
    "Xred = X[y==1]\n",
    "\n",
    "# Display\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b',s=20)#, edgecolors='k')\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r',s=20);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "55633a07-f69f-4692-9eb3-18d5787ef78e"
    }
   },
   "source": [
    "We still would like to fit a line through the data, since we can assume the overlapping of the red and blue dots is due to noise. In this case, the optimal separating hyperplane can be found by making a **compromise** between the amount of misclassifications and the value of the margin.\n",
    "\n",
    "For each pair $(x_i,y_i)$, we introduce the slack variable $\\xi_i$ that represents how far $x_i$ is on the wrong side of the margin's boundary, as illustrated by the figure below.\n",
    "\n",
    "<img width=\"300px\" src=\"img/non_lin_sep3.png\">"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "55633a07-f69f-4692-9eb3-18d5787ef78e"
    }
   },
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**<br>\n",
    "Let's write $y_i(w^T x_i + w_0) \\geq M(1-\\xi_i)$ and note that $x_i$ is misclassified if $\\xi_i\\geq1$.\n",
    "Reformulate the SVM optimization problem so that we have a maximum \"budget\" of $K$ misclassifications.\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "55633a07-f69f-4692-9eb3-18d5787ef78e"
    }
   },
   "source": [
    "<details class=\"alert alert-danger\">\n",
    "    <summary markdown=\"span\"><b>Ready to see the answer? (click to expand)</b></summary>\n",
    "\n",
    "The problem simply turns into:\n",
    "$$\\min_{w,w_0} \\frac{1}{2}\\|w\\|^2$$\n",
    "$$y_i (w^T x_i + w_0) \\geq 1-\\xi_i, \\forall i \\in [1,N]$$\n",
    "$$\\xi_i \\geq 0$$\n",
    "$$\\sum\\limits_{i=1}^N \\xi_i \\leq K$$\n",
    "</details>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "55633a07-f69f-4692-9eb3-18d5787ef78e"
    }
   },
   "source": [
    "This last formulation with $K$ is not satisfying: instead of having to impose a constraint on the number of misclassifications, we would rather make a compromise between the number of misclassifications and the value of the margin. So the optimization problem turns into:\n",
    "\n",
    "$$ \\min\\limits_{w,w_0,\\xi} \\displaystyle\\frac{1}{2}\\|w\\|^2 + C\\sum\\limits_{i=1}^N \\xi_i$$\n",
    "\n",
    "$$\\text{ such that } \\forall i=1..N, $$\n",
    "$$y_i\\left(w^T x_i + w_0 \\right)\\geq 1-\\xi_i, $$\n",
    "$$\\xi_i\\geq 0 $$\n",
    "\n",
    "Note that $C$ is the parameter that controls the trade-off between having a large margin and having few misclassified points.\n",
    "\n",
    "This is again a quadratic programming problem. Using the $\\alpha_i$ and $\\mu_i$ Lagrange multipliers, the KKT conditions yield:\n",
    "\n",
    "\n",
    "$$L_P = \\frac{1}{2}\\|w\\|^2 + C\\sum\\limits_{i=1}^N \\xi_i - \\sum\\limits_{i=1}^N \\alpha_i\\left(y_i\\left(w^T x_i + w_0 \\right)- \\left(1-\\xi_i\\right)\\right) - \\sum\\limits_{i=1}^N \\mu_i \\xi_i$$\n",
    "\n",
    "$$\\text{ KKT conditions :} $$\n",
    "$$\\frac{\\partial L_P}{\\partial w} = 0 \\Rightarrow w = \\sum\\limits_{i=1}^N \\alpha_i y_i x_i$$\n",
    "$$\\frac{\\partial L_P}{\\partial w_0} = 0 \\Rightarrow 0 = \\sum\\limits_{i=1}^N \\alpha_i y_i$$\n",
    "$$\\frac{\\partial L_P}{\\partial \\xi} = 0 \\Rightarrow \\alpha_i = C-\\mu_i$$\n",
    "$$\\forall i=1..N, \\ \\alpha_i\\left(y_i\\left(w^T x_i + w_0 \\right)-\\left(1-\\xi_i\\right)\\right) = 0$$\n",
    "$$\\forall i=1..N, \\ \\mu_i\\xi_i = 0$$\n",
    "$$\\forall i=1..N, \\alpha_i\\geq 0, \\mu_i \\geq 0$$\n",
    "\n",
    "The dual problem is almost the same as previously (it eliminates all variables except for the $\\alpha_i$) and can be solved using SMO:\n",
    "\n",
    "$$\\max\\limits_{\\alpha\\in{\\mathbb{R}^+}^N} L_D(\\alpha) = \\sum\\limits_{i=1}^N \\alpha_i - \\frac{1}{2}\\sum\\limits_{i=1}^N\\sum\\limits_{j=1}^N \\alpha_i \\alpha_j y_i y_j x_i^T x_j$$\n",
    "$$\\text{such that } \\sum\\limits_{i=1}^N \\alpha_i y_i = 0$$\n",
    "$$\\text{and } 0 \\leq \\alpha_i \\leq C$$\n",
    "\n",
    "Let's follow the same reasoning as previously on the $\\alpha_i\\left(y_i\\left(w^T x_i + w_0 \\right)-\\left(1-\\xi_i\\right)\\right) = 0$ constraint:\n",
    "- $\\alpha_i>0$, then $y_i\\left(w^T x_i + w_0 \\right)=1-\\xi_i$: $x_i$ is a **support vector**.<br>\n",
    "Among the support vectors:\n",
    "    - $\\xi_i=0$, then $0 < \\alpha_i \\leq C$\n",
    "    - $\\xi_i>0$, then $\\alpha_i=C$ (because $\\mu_i=0$, because $\\mu_i\\xi_i=0$)\n",
    "- $\\alpha_i=0$, then $x_i$ does not participate in $w$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "55633a07-f69f-4692-9eb3-18d5787ef78e"
    }
   },
   "source": [
    "## <a id=\"sec1-5\"></a> 1.5 Making predictions\n",
    "\n",
    "So, suppose we solve this optimization problem and get $(w,w_0)$, how do we make a prediction on a new data point $x$? Simple! Just compute $w^T x + w_0$: if it's positive we are on one side of the classifier, otherwise we are on the other side."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "575227d2-0b4e-43fd-a4f3-c32f6163a9ed"
    }
   },
   "source": [
    "# 2. <a id=\"sec2\"></a>Support Vector Machines in scikit-learn\n",
    "\n",
    "Fortunately for us, this optimization problem is solved is three lines in scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "6dce34d9-ee79-4aaa-86c1-274cce55dec2"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn import svm\n",
    "mySVC = svm.SVC(kernel='linear', C=1)\n",
    "mySVC.fit(X,y)\n",
    "\n",
    "# Compute margin and find support vectors\n",
    "w = mySVC.coef_[0]\n",
    "w0 = mySVC.intercept_\n",
    "M = 1./np.linalg.norm(w)\n",
    "print(\"SV per class:\", mySVC.n_support_)\n",
    "print(\"w_i:\", w)\n",
    "print(\"w_0:\", w0)\n",
    "print(\"Margin:\", M)\n",
    "print(\"w^T x0 + w_0:\", np.dot(w,mySVC.support_vectors_[0,:])+w0)\n",
    "print(\"w^T x1 + w_0:\", np.dot(w,mySVC.support_vectors_[1,:])+w0)\n",
    "print(\"w^T x2 + w_0:\", np.dot(w,mySVC.support_vectors_[2,:])+w0)\n",
    "\n",
    "# Plot the separating plane, the margin and the Support Vectors\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b',s=20)\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r',s=20)\n",
    "XX = np.arange(-1.,12.,0.1)\n",
    "YY = -(w[0]*XX+w0)/w[1]\n",
    "plt.plot(XX,YY,'g')\n",
    "YY = -(w[0]*XX+w0+1)/w[1]\n",
    "plt.plot(XX,YY,'g--')\n",
    "YY = -(w[0]*XX+w0-1)/w[1]\n",
    "plt.plot(XX,YY,'g--')\n",
    "plt.scatter(mySVC.support_vectors_[:,0], mySVC.support_vectors_[:,1], s=80, edgecolors='k', facecolors='none');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "23cb012c-0452-4012-8595-eeb54e756485"
    }
   },
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice**  \n",
    "Try different values for $C$ in the code above to see how the margin's boundaries evolve.\n",
    "</div>\n",
    "\n",
    "Remark: since the optimization problem behind SVMs requires the knowledge of the full dataset, SVMs are necessarily an offline method.\n",
    "\n",
    "Let's plot how the number of support vectors change with $C$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "e5ee8e83-c4af-4aba-a915-8282132cc5ac"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Recall, there are\", Xblue.shape[0], \"blue points and\", Xred.shape[0], \"red points.\")\n",
    "logC = np.arange(-4,7,1)\n",
    "nbSV = np.zeros(len(logC))\n",
    "for i in range(len(logC)):\n",
    "    print(\"Training at C =\", 10.**logC[i])\n",
    "    mySVC = svm.SVC(kernel='linear', C=10.**logC[i])\n",
    "    mySVC.fit(X,y)\n",
    "    nbSV[i] = np.sum(mySVC.n_support_)\n",
    "\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.plot(logC,nbSV)\n",
    "plt.xlabel(\"log C\")\n",
    "plt.ylabel(\"nb support vectors\");"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "7fc6c322-c6f0-461a-a5dd-9cf34eb1ba13"
    }
   },
   "source": [
    "Let's make a few predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "92964068-7386-49e4-8fe9-796899ee0576"
    }
   },
   "outputs": [],
   "source": [
    "mySVC = svm.SVC(kernel='linear', C=1)\n",
    "mySVC.fit(X,y)\n",
    "\n",
    "# Compute margin and find support vectors\n",
    "w = mySVC.coef_[0]\n",
    "w0 = mySVC.intercept_\n",
    "M = 1./np.linalg.norm(w)\n",
    "print(\"SV per class:\", mySVC.n_support_)\n",
    "print(\"w_i:\", w)\n",
    "print(\"w_0:\", w0)\n",
    "print(\"Margin:\", M)\n",
    "print(\"w^T x0 + w_0:\", np.dot(w,mySVC.support_vectors_[0,:])+w0)\n",
    "print(\"w^T x1 + w_0:\", np.dot(w,mySVC.support_vectors_[1,:])+w0)\n",
    "print(\"w^T x2 + w_0:\", np.dot(w,mySVC.support_vectors_[2,:])+w0)\n",
    "\n",
    "# Plot the separating plane, the margin and the Support Vectors\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b',s=20)\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r',s=20)\n",
    "XX = np.arange(-1.,12.,0.1)\n",
    "YY = -(w[0]*XX+w0)/w[1]\n",
    "plt.plot(XX,YY,'g')\n",
    "YY = -(w[0]*XX+w0+1)/w[1]\n",
    "plt.plot(XX,YY,'g--')\n",
    "YY = -(w[0]*XX+w0-1)/w[1]\n",
    "plt.plot(XX,YY,'g--')\n",
    "\n",
    "# Testing data\n",
    "Xtest = np.array([[2,4],[6,2],[10,2]])\n",
    "Ypred = mySVC.predict(Xtest)\n",
    "for i in range(Xtest.shape[0]):\n",
    "    print(\"Prediction in\", Xtest[i,:], \"=\", Ypred[i])\n",
    "\n",
    "plt.scatter(Xtest[:,0],Xtest[:,1], c='m', edgecolors=\"k\", s=80);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "83029f02-43d2-428d-b957-2fe7ceceffab"
    }
   },
   "source": [
    "# <a id=\"sec3\"></a>3. When using linear separators makes no more sense\n",
    "\n",
    "Let's take a look at the following data. Would it make sense to try to separate the red points from the blue ones with a straight line?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "33109de0-9ac6-4138-8574-b26f067ac171"
    }
   },
   "outputs": [],
   "source": [
    "X1, y1 = datasets.make_gaussian_quantiles(cov=2.,\n",
    "                                 n_samples=300, n_features=2,\n",
    "                                 n_classes=2, random_state=1)\n",
    "X2, y2 = datasets.make_gaussian_quantiles(mean=(3, 3), cov=1.5,\n",
    "                                 n_samples=700, n_features=2,\n",
    "                                 n_classes=2, random_state=1)\n",
    "X = np.concatenate((X1, X2))\n",
    "y = np.concatenate((y1, - y2 + 1))\n",
    "y = 2*y-1\n",
    "\n",
    "X, y = shuffle(X, y)\n",
    "\n",
    "Xtest,X = np.split(X,[400])\n",
    "ytest,y = np.split(y,[400])\n",
    "\n",
    "Xblue = X[y==-1]\n",
    "Xred = X[y==1]\n",
    "fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "plt.scatter(Xblue[:,0],Xblue[:,1],c='b')\n",
    "plt.scatter(Xred[:,0],Xred[:,1],c='r');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "973f28e4-dc48-47bf-b401-f1c1d14938c5"
    }
   },
   "source": [
    "It does not look like a linear separation even makes sense. But let's try nonetheless, some optimal separating hyperplane must exist (even if it does not have great generalization properties).\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**  \n",
    "Use scikit-learn to compute a linear Support Vector Classifier for the data above. Compute the classification score and conclude on the limits of this approach.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "d08cd339-291e-43cb-8689-0e83ee74843f"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/code3.py\n",
    "### WRITE YOUR CODE HERE\n",
    "# If you get stuck, uncomment the line above to load a correction in this cell (then you can execute this code).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "d8994a10-ee7b-47a2-892e-8b195de373bb"
    }
   },
   "source": [
    "Quite unconvincing, isn't it?"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A short metaphor before we go into more details**\n",
    "\n",
    "<img src=\"img/crumpled-paper.jpg\" width=\"400px\"></img>\n",
    "\n",
    "Take two blank sheets of paper. On the first one, draw a set of blue points. On the second, draw some red points. Now put one page on top of the other and crumple them together. The red and blue points are separable (whatever number of points there were on each sheet) , it's a certainty, since you put them on purpose on two different pieces of paper (two classes). But if you're given the current coordinates of each point, then building a relevant classifier might not be so easy. Some ideas follow from this intuition:\n",
    "- Good classification might be feasible if we could build an inverse transformation (\"uncrumpling\") of the data that makes it linearly separable.\n",
    "- Good classification functions should be able to approximate a large variety of geometric transformations of the data (this inverse transformation should be able to capture a large variety of of crumpling patterns).\n",
    "- A classification method that can only approximate a few geometrical transformations is prone to miss the structure of the data. We call that *underfitting*, which is due to *high bias* in the family of classification functions.\n",
    "- Having a large set of uncrumpling patterns (being able to approximate many geometrical transformations of the data) comes at the price of taking the risk of picking the wrong uncrumpling pattern given the finite amount of data at our disposal. We call that *overfitting*, which is due to *high variance* in the family of classification functions."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "8eea85d1-7b04-4014-a887-40b78c263bd7"
    }
   },
   "source": [
    "# <a id=\"sec4\"></a>4. A word on the bias-variance tradeoff\n",
    "\n",
    "The previous example is unconvincing mainly because we stubbornly look for a separator that has the shape of a straight line. Maybe if we allowed more flexibility, we could get better separators?\n",
    "\n",
    "This compromise between flexibility and prior knowledge is called the **bias-variance** tradeoff in Machine Learning.\n",
    "- If I search for a function that fits my data within a very restricted set of functions, I introduce a lot of expert knowledge, and thus a lot of **bias**, only allowing for little **variance** in the functions I might find.\n",
    "- On the other hand, if I use a very rich family of functions, I abandon the introduction of prior information, allowing for high **variance** and low **bias**.\n",
    "\n",
    "Although appealing, having high variance is not necessarily a good thing: looking for a function in a set that is just too rich leads to the risk of finding several almost equivalent functions that have great statistical properties but no physical meaning.\n",
    "\n",
    "By maximizing the margin, SVMs actually try to minimize the variance of the family of functions considered: that is a way of trying to find the \"simplest\" function possible by minimizing $\\|w\\|$. This process is called **regularization** in Machine Learning. SVMs perform an $L_2$ regularization, but other forms of regularization exist and are commonly being used (e.g. $L_1$ regularization).\n",
    "\n",
    "The bias-variance tradeoff captures the problem of **overfitting** (due to high variance) versus **underfitting** (due to high bias). Let's try to capture this formally. Assume a supervised learning task where one has data $y = f(x) + \\epsilon$, where $f$ is a deterministic function and $\\epsilon$ is a zero-mean noise with variance $\\sigma^2$. The mean square generalization error of a prediction function $\\hat{f}$ is:\n",
    "$$err =  \\mathbb{E}_x\\left( \\mathbb{E}_{y\\sim f(x)+\\epsilon} \\left( \\left( y - \\hat{f}(x)\\right)^2 \\right) \\right) = \\mathbb{E} _x \\left( err(x) \\right)$$\n",
    "with:\n",
    "$$err(x) = \\mathbb{E}_{y\\sim f(x)+\\epsilon} \\left( \\left( y - \\hat{f}(x)\\right)^2 \\right)$$\n",
    "\n",
    "Let's write $\\mathbb{E}=\\mathbb{E}_{y\\sim f(x)+\\epsilon}$, $f=f(x)$ and $\\hat{f} = \\hat{f}(x)$ for simplicity.\n",
    "\n",
    "We have then:\n",
    "\\begin{align}\n",
    "err(x) & = \\mathbb{E} \\left[ y^2 + \\hat{f}^2 - 2y\\hat{f} \\right]\\\\\n",
    "    & = \\mathbb{E} \\left[y^2\\right] + \\mathbb{E}\\left[ \\hat{f}^2\\right] -2 \\mathbb{E}\\left[y\\hat{f} \\right]\n",
    "\\end{align}\n",
    "But, by definition of the variance:\n",
    "- $\\mathbb{E} \\left[y^2\\right] = Var\\left(y\\right) + \\mathbb{E} \\left[y\\right]^2 = \\sigma^2 + f^2$\n",
    "- $\\mathbb{E}\\left[ \\hat{f}^2\\right] = Var\\left(\\hat{f}\\right) + \\mathbb{E} \\left[\\hat{f}\\right]^2$\n",
    "- and $f$ is deterministic so $\\mathbb{E}\\left[y\\hat{f} \\right] = f\\mathbb{E}\\left[\\hat{f}\\right]$.\n",
    "\n",
    "So:\n",
    "\\begin{align}\n",
    "err(x) & = \\sigma^2 + f^2 + Var(\\hat{f}) + \\mathbb{E}\\left[\\hat{f}\\right]^2 - 2 f\\mathbb{E}\\left[\\hat{f}\\right]\\\\\n",
    "    & = \\sigma^2 + Var(\\hat{f}) + \\left(f^2 -2 f\\mathbb{E}\\left[\\hat{f}\\right] + \\mathbb{E}\\left[\\hat{f}\\right]^2 \\right)\\\\\n",
    "    & = \\sigma^2 + Var(\\hat{f}) + \\left( f - \\mathbb{E}\\left[\\hat{f}\\right] \\right)^2\\\\\n",
    "    & = \\textrm{Irreducible error} + \\textrm{Variance} + \\textrm{Bias}^2\n",
    "\\end{align}\n",
    "\n",
    "Let's explain each of these terms:\n",
    "- The irreducible error $\\sigma^2$ is the noise term in the true relationship that cannot fundamentally be reduced by any  learning model.\n",
    "- The bias $f(x) - \\mathbb{E}\\left[\\hat{f}(x)\\right]$ measures how far from $f$ the best fit of $\\hat{f}$ is.\n",
    "- The variance $\\mathbb{E}\\left[ \\hat{f}^2\\right] - \\mathbb{E} \\left[\\hat{f}\\right]^2$ of $\\hat{f}$ measures how much $\\hat{f}$ is able to move around its mean.\n",
    "\n",
    "So, overall, our goal is to reduce both variance and bias to minimize the generalization error. But in a world with finite data and no information on the true shape of $f$, there is a compromise to make between minimizing the bias or the variance.\n",
    "\n",
    "SVMs make this compromise through the $C$ parameter: small $C$ will give more importance to maximizing the margin, thus minimizing the variance (and conversely, large $C$ values will emphasize the data fitting term and thus will minimize the bias)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "42fadc80-3a4c-4a4c-9425-043e104c62f8"
    }
   },
   "source": [
    "# <a id=\"sec5\"></a>5. The kernel trick\n",
    "\n",
    "## <a id=\"sec5-1\"></a> 5.1 The intuition of mapping to higher dimensions\n",
    "\n",
    "Let's get back to the idea of giving more flexibility to our separator, since hyperplanes did not convince us they could always do the job.\n",
    "\n",
    "To get around this problem, let's first introduce an illustrative example. The data below come from a voltage test in electronics. They indicate if a component fails or not under a certain voltage ($U$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "f25a5d30-b13b-4532-85bd-f61b89f8315e"
    }
   },
   "outputs": [],
   "source": [
    "U = np.array([[0.3, 0.7, 1.1, 1.8, 2.5, 3.0, 3.3, 3.5, 3.7]]).T\n",
    "Good = np.array( [[ -1,  -1,  -1,   1,   1,   1,   1,  -1,  -1]]).T\n",
    "plt.figure()\n",
    "plt.scatter(U, np.zeros((U.shape[0],1)), c=Good, cmap = plt.cm.coolwarm);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "217e5e28-6b8d-4e79-9147-cdbb61f790ab"
    }
   },
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Question:**  \n",
    "Does it look like these data are linearly separable?</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "217e5e28-6b8d-4e79-9147-cdbb61f790ab"
    }
   },
   "source": [
    "Some well-experienced engineer already knows that $U$ is no good criterion to split the data points into two categories with a single threshold (recall: in dimension 1, a hyperplane is a threshold). He also knows that other indicators like $V = (U-2)^2$ are not any better, but that the pair $(U,V)$ actually allows to build a separating line between points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "7547c3d6-ccd0-49cd-97a9-a92ccf7b20b2"
    }
   },
   "outputs": [],
   "source": [
    "V = (2-U)**2\n",
    "plt.figure()\n",
    "plt.scatter(U, V, c=Good, cmap = plt.cm.coolwarm)\n",
    "XX = np.arange(np.min(U), np.max(U), 0.05)\n",
    "YY = -.8+.8*XX\n",
    "plt.plot(XX,YY,'g');"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "7ad8607f-7316-49b0-a25c-38924256edaf"
    }
   },
   "source": [
    "Our nine points are now linearly separable. But we have not introduced any additional information!<br>\n",
    "<br>\n",
    "Actually, the only thing we have done, is to *map* our non-separable points $x$ into **another description space**, of higher dimension, where a linear separator exists.<br>\n",
    "\n",
    "## <a id=\"sec5-2\"></a> 5.2 The kernel trick\n",
    "\n",
    "Suppose we happen to know such a relevant mapping for our data $\\varphi(x)=x'$, where $\\varphi$ maps to a $p$-dimensional Euclidean space ($p\\gg n$, possibly infinite). Then we can compute the optimal linear separator in the corresponding higher dimension description space, find its parameters $w'$ and $w'_0$, and whenever we need to make a new prediction in a point $x$, we first compute its image $x'=\\varphi(x)$ and then calculate $w'^T x' + w'_0$ to know on which side of the hyperplane we stand.\n",
    "\n",
    "But recall that: $$w' = \\sum_{i=1}^N \\alpha_i y_i x'_i$$\n",
    "\n",
    "So :\n",
    "\\begin{align*}\n",
    "w'^T x' + w'_0 & = \\left(\\sum_{i=1}^N \\alpha_i y_i x'_i\\right)^T \\varphi(x) + w_0\\\\\n",
    "& = \\sum_{i=1}^N \\alpha_i y_i \\varphi(x_i)^T \\varphi(x) + w'_0\n",
    "\\end{align*}\n",
    "\n",
    "Suppose that, instead of providing us with a mapping $\\varphi(x)=x'$, somebody gave us a function $K(x_1, x_2)$ that takes two points $x_1$ and $x_2$, computes their respective images $\\varphi(x_1)$ and $\\varphi(x_2)$ and returns the dot product:\n",
    "$$K(x_1,x_2)=\\varphi(x_1)^T \\varphi(x_2)$$\n",
    "\n",
    "Then making predictions would boil down to:\n",
    "$$w'^T x' + w'_0 = \\sum_{i=1}^N \\alpha_i y_i K(x_i,x) + w'_0$$\n",
    "\n",
    "The interesting thing is that we don't need to compute $\\varphi$ anymore. The function $K$ is known as a **kernel function** and that's what we call the **kernel trick**.\n",
    "\n",
    "It is actually possible to compute the $\\alpha_i$ and $w'_0$ just using $K$ and never $\\varphi$, so, as long as somebody insures that **the kernel $K$ is a dot product in some other descriptor space** we can compute the optimal separating hyperplane of our data in this space, **without ever requiring a knowledge of this descriptor space and the mapping $\\varphi$**.\n",
    "\n",
    "## <a id=\"sec5-3\"></a> 5.3 Positive definite kernels\n",
    "\n",
    "The main question is thus: when can we guarantee that $K$ is an acceptable kernel? It is so if it is an inner product on a (separable) Hilbert space. So in more general words, we are interested in positive, definite kernels on a Hilbert space:\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "7ad8607f-7316-49b0-a25c-38924256edaf"
    }
   },
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice (off-class):**  \n",
    "Check for yourself that the dual optimization problem for the SVM in the kernel's Hilbert space can be written:\n",
    "$$\\max\\limits_{\\alpha\\in{\\mathbb{R}^+}^N} L_D(\\alpha) = \\sum\\limits_{i=1}^N \\alpha_i - \\frac{1}{2}\\sum\\limits_{i=1}^N\\sum\\limits_{j=1}^N \\alpha_i \\alpha_j y_i y_j K(x_i,x_j)$$\n",
    "$$\\text{such that } \\sum\\limits_{i=1}^N \\alpha_i y_i = 0$$\n",
    "$$\\text{and } 0 \\leq \\alpha_i \\leq C$$\n",
    "So its resolution only requires the knowledge of $K$.\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "$K(\\cdot,\\cdot)$ is a positive definite kernel on $X$ if\n",
    "$$\\forall n\\in\\mathbb{N}, x\\in X^n \\text{ and } c\\in \\mathbb{R}^n, \\ \\sum\\limits_{i,j=1}^n c_i c_j K(x_i,x_j) \\geq 0$$\n",
    "\n",
    "\n",
    "We will admit that Mercer's condition below is a sufficient condition for $K$ to be a positive definite kernel:\n",
    "\n",
    "Given $K(x,y)$, if:\n",
    "$$\\forall g(x) / \\int g(x)^2dx <\\infty, \\iint K(x,y)g(x)g(y)dxdy \\geq 0$$\n",
    "Then, there exists a mapping $h(\\cdot)$ such that:\n",
    "$$K(x,y) = \\langle h(x), h(y) \\rangle$$\n",
    "\n",
    "There are many kernels that have been developped in the litterature. Combined with what we have seen before, it allows to build non-linear SVMs. The nice thing is that some kernels actually map our data to a descriptor space of infinite dimension, where it is presumably a lot easier to find a separating hyperplane.\n",
    "\n",
    "One such kernel is the so-called \"radial basis kernel\" which is very popular and can be written:\n",
    "$$K(x_1,x_2) = e^{-\\gamma \\|x_1-x_2\\|^2}$$\n",
    "\n",
    "Other common kernels:\n",
    "- polynomial $K(x,y)=\\left(1+\\langle x, y\\rangle\\right)^d$\n",
    "- sigmoid $K(x,y) = \\tanh\\left(\\kappa_1 \\langle x, y\\rangle + \\kappa_2\\right)$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "bf8918ea-6de2-4759-a168-32990eb1de38"
    }
   },
   "source": [
    "# <a id=\"sec6\"></a>6. SVMs and kernels\n",
    "\n",
    "Let's practice on the last dataset with a radial basis kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "fa7dbe63-50b7-4b7b-881b-fc1128e18dcc"
    }
   },
   "outputs": [],
   "source": [
    "# Plot separator, margin and support vectors\n",
    "def plot_SVC(mySVC):\n",
    "    fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "    plt.scatter(Xblue[:,0],Xblue[:,1],c='b')\n",
    "    plt.scatter(Xred[:,0],Xred[:,1],c='r')\n",
    "    XX, YY = np.meshgrid(np.arange(np.min(X[:,0]),np.max(X[:,0]),0.1), np.arange(np.min(X[:,1]),np.max(X[:,1]),0.1))\n",
    "    ZZ = mySVC.decision_function(np.c_[XX.ravel(), YY.ravel()])\n",
    "    ZZ = ZZ.reshape(XX.shape)\n",
    "    plt.contour(XX, YY, ZZ, levels=[0],alpha=0.75)\n",
    "    fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "    cont = plt.contour(XX, YY, ZZ, levels=[-1., 0., 1.], alpha=0.75)\n",
    "    plt.clabel(cont, cont.levels, inline=True, fontsize=10)\n",
    "    fig=plt.figure(figsize=fig_size, dpi= 80, facecolor='w', edgecolor='k')\n",
    "    cont = plt.contourf(XX, YY, ZZ, alpha=0.75, cmap = plt.cm.coolwarm)\n",
    "\n",
    "# We can play with C\n",
    "mySVC = svm.SVC(kernel='rbf', C=1.)\n",
    "mySVC.fit(X,y)\n",
    "\n",
    "print(\"SV per class:\", mySVC.n_support_)\n",
    "plot_SVC(mySVC)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "71301663-2487-40a8-87c6-74a4d1183162"
    }
   },
   "source": [
    "A lot better isn't it? Let's play a little with the value of $C$.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**  \n",
    "Let's see how that SVM with the rbf kernel would do on the data from the beginning of this notebook.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "nbpresent": {
     "id": "e4dc7f1d-2d84-46bb-aa86-6571b3b7b9e9"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/code4.py\n",
    "### WRITE YOUR CODE HERE\n",
    "# If you get stuck, uncomment the line above to load a correction in this cell (then you can execute this code).\n",
    "\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "92ff468b-8cd1-4b92-acbb-b13a3e540fee"
    }
   },
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**  \n",
    "And on the \"linearly separable + noise\" case?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "21e389bc-335d-4b93-aca8-0637b709b2df"
    }
   },
   "outputs": [],
   "source": [
    "X1, y1 = datasets.make_gaussian_quantiles(cov=2.0, n_samples=300, n_features=2, n_classes=1)\n",
    "X1[:,0] = 3. + X1[:,0]\n",
    "X1[:,1] = 6. + X1[:,1]/2.5\n",
    "X2, y2 = datasets.make_gaussian_quantiles(cov=1.5, n_samples=300, n_features=2, n_classes=1)\n",
    "X2[:,0] = 8. + X2[:,0]\n",
    "X2[:,1] = 4. + X2[:,1]\n",
    "X3, y3 = datasets.make_gaussian_quantiles(cov=1.5, n_samples=300, n_features=2, n_classes=1)\n",
    "X3[:,0] = 7. + X3[:,0]\n",
    "X3[:,1] = 8. + X3[:,1]\n",
    "X = np.concatenate((X1, X2, X3))\n",
    "y = np.concatenate((y1, - y2 + 1, y3))\n",
    "y = 2*y-1\n",
    "X, y = shuffle(X, y)\n",
    "Xblue = X[y==-1]\n",
    "Xred = X[y==1]\n",
    "\n",
    "# Display\n",
    "mySVC = svm.SVC(kernel='rbf')\n",
    "mySVC.fit(X,y)\n",
    "print(\"SV per class:\", mySVC.n_support_)\n",
    "plot_SVC(mySVC);"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "cc280e3d-b852-4b4b-89a1-238b30cf1441"
    }
   },
   "source": [
    "These two last experiments illustrate that introducing complex kernels such as the rbf one is sometimes detrimental."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "1a581062-808b-4e72-9de2-f0fec07fabeb"
    }
   },
   "source": [
    "# <a id=\"sec7\"></a>7. What about other uses?\n",
    "\n",
    "SVMs can be adapted to deal with **multi-class** classification by training \"one-versus-one\" or \"one-versus-rest\" binary SVMs. In the first case, for $C$ classes, $C(C-1)/2$ classifiers need to be built and a majority vote is implemented. In the second case, only $C$ classifiers are built, at the risk of class imbalance.\n",
    "\n",
    "SVMs can also be used to solve **regression** problems, although they are not often used for that in practice. The optimization problem differs a little but the main ideas remain. In scikit-learn, they can be accessed by using the `svm.SVR` class.\n",
    "\n",
    "In the regression case, we do not want to separate anymore but to *fit* the hyperplane to the data. For this purpose, we define a *loss function* $\\sum_i V(y_i - f(x_i))$ and the objective function becomes:\n",
    "$$\\min\\limits_{w, w_0} \\frac{1}{2} \\|w\\|^2 + C \\sum\\limits_{i=1}^N V(y_i-w^T x_i + w_0))$$\n",
    "\n",
    "Usual loss functions:\n",
    "- $\\epsilon$-insensitive, $V(z) = \\left\\{\\begin{array}{l} 0\\text{ if }|z|\\leq\\epsilon\\\\ |z|-\\epsilon\\text{ otherwise}\\end{array}\\right.$\n",
    "- Laplacian, $V(z) = |z|$\n",
    "- Gaussian, $V(z) = \\frac{1}{2}z^2$\n",
    "- Huber's robust loss, $V(z) = \\left\\{\\begin{array}{l}\\frac{1}{2\\sigma}z^2\\text{ if }|z|\\leq\\sigma \\\\ |z|-\\frac{\\sigma}{2}\\text{ otherwise}\\end{array}\\right.$\n",
    "\n",
    "\n",
    "The $\\epsilon$-insensitive loss:  \n",
    "<img src=\"img/epsilon_insensitive2.png\" width=\"500px\"></img>\n",
    "\n",
    "The optimization problem and its resolution are quite similar to the classification case (but a little heavier, notation-wise) and left to your curiosity, just as the extension to the automation of the search for the trade-off parameter $C$, through $\\nu$-SVR.\n",
    "\n",
    "Similarly, one can derive a formulation of SVMs for **density estimation** by using *One-Class SVMs*. If you are curious to know more, please refer to the references mentionned a few lines below.\n",
    "\n",
    "SVMs easily support class weights and sample weights and can deal with unbalanced problems.\n",
    "\n",
    "The general description of [SVMs in scikit-learn](http://scikit-learn.org/stable/modules/svm.html) is a good reminder on the versatility and the limitations of SVMs."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "16365b90-7747-494a-bdf9-8f23f5c301cb"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "**In a nutshell:**  \n",
    "- Support Vector Machines for classification try to separate data by maximizing a geometrical margin\n",
    "- They are computed offline\n",
    "- They offer a sparse, robust to class imbalance, and easy to evaluate predictor\n",
    "- Kernels are a way of enriching (lifting) the data representation so that it becomes linearly separable\n",
    "- SVMs + kernels offer a versatile method for classification, regression and density estimation\n",
    "- [Documentation in scikit-learn](http://scikit-learn.org/stable/modules/svm.html)\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "16365b90-7747-494a-bdf9-8f23f5c301cb"
    }
   },
   "source": [
    "**References**\n",
    "\n",
    "On the general theory of SVMs for classification:  \n",
    "**A tutorial on Support Vector Machines for Pattern Recognition.**  \n",
    "C. J. C. Burges, *Data Mining and Knowledge Discovery*, **2**, 131-167, (1998).\n",
    "\n",
    "On support vector regression (and its extension to $\\nu$-SVR):  \n",
    "**A tutorial on Support Vector Regression.**  \n",
    "A. J. Smola and B. Schölkopf, *Journal of Statistics and Computing*, **14**(3), 199-222, (2004).  \n",
    "**New support vector algorithms.**\n",
    "B. Schölkopf, A. J. Smola, R. C. Williamson, and P. L. Bartlett. *Neural computation*, **12**(5), 1207-1245,  (2000).\n",
    "\n",
    "On One-Class SVMs:  \n",
    "**Support vector method for novelty detection.**  \n",
    "B. Schölkopf, R. C. Williamson, A. J. Smola, J. Shawe-Taylor, and John C. Platt. *Neural Information Processing Systems*, **12**, 582-588, (1999).\n",
    "\n",
    "On multi-class SVMs:  \n",
    "**On the algorithmic implementation of multiclass kernel-based vector machines.**  \n",
    "K. Crammer and Y. Singer. *Journal of machine learning research*, **2**, 265-292, (2001)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "109dd2fd-c49e-44a9-8a48-911877fce53e"
    }
   },
   "source": [
    "# <a id=\"sec8\"></a>8. Examples\n",
    "\n",
    "## <a id=\"sec8-1\"></a>8.1 Spam or ham?\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**  \n",
    "Your turn to play: let's build a spam classifier using the ling-spam dataset and a linear SVM. Compute it's score on the validation dataset below.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "67bbc140-344b-4510-8635-72df3e783e70"
    }
   },
   "outputs": [],
   "source": [
    "from sys import path\n",
    "path.append('../2 - Text data preprocessing')\n",
    "import load_spam\n",
    "spam_data = load_spam.spam_data_loader()\n",
    "spam_data.load_data()\n",
    "\n",
    "Xtrain, ytrain, Xtest, ytest = spam_data.split(2000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(Xtrain.shape)\n",
    "print(Xtest.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "9e414752-9438-46d7-b036-e2af3e670979"
    }
   },
   "outputs": [],
   "source": [
    "# %load solutions/code5.py\n",
    "### WRITE YOUR CODE HERE\n",
    "# If you get stuck, uncomment the line above to load a correction in this cell (then you can execute this code).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Exercice:**  \n",
    "We've trained our model in the Tf-Idf data. Let's see how the model behaves on raw word counts.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xtrain, ytrain, Xtest, ytest = spam_data.split(2000, feat='wordcount')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# %load solutions/code6.py\n",
    "### WRITE YOUR CODE HERE\n",
    "# If you get stuck, uncomment the line above to load a correction in this cell (then you can execute this code).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "    \n",
    "**Exercice:** Want to try another kernel than linear?</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's keep the tf-idf with linear SVM classifier (which seems to work better) and use it to identify which are the misclassified emails (and find the confusion matrix by the way)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "873f04f7-2a40-45f9-9a03-970f71cd5928"
    }
   },
   "outputs": [],
   "source": [
    "# Retrain\n",
    "Xtrain, ytrain, Xtest, ytest = spam_data.split(2000)\n",
    "spam_svc = svm.SVC(kernel='linear', C=1.)\n",
    "spam_svc.fit(Xtrain,ytrain);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "bab94b21-0f19-4ab5-81c6-7601ced11400"
    }
   },
   "outputs": [],
   "source": [
    "# Find misclassified examples\n",
    "ypredict = spam_svc.predict(Xtest)\n",
    "misclass = np.not_equal(ypredict, ytest)\n",
    "Xmisclass = Xtest[misclass,:]\n",
    "ymisclass = ytest[misclass]\n",
    "misclass_indices = [i for i, j in enumerate(misclass) if j == True]\n",
    "print(\"Misclassified messages indices:\", misclass_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "print(\"Confusion matrix:\")\n",
    "print(confusion_matrix(ytest, ypredict))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check some misclassified mails\n",
    "index = misclass_indices[0]+2000\n",
    "print(\"Prediction:\", spam_svc.predict(spam_data.tfidf[index,:]))\n",
    "spam_data.print_email(index)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "710e351f-00b6-41bc-aff8-a7a203435248"
    }
   },
   "source": [
    "## <a id=\"sec8-2\"></a>8.2. NIST\n",
    "\n",
    "Let's evaluate SVMs on the optical character recognition task of the NIST data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "b5946f3a-cce1-4138-af17-0fb2cb788f38"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "digits = datasets.load_digits()\n",
    "print(digits.data.shape)\n",
    "print(digits.images.shape)\n",
    "print(digits.target.shape)\n",
    "print(digits.target_names)\n",
    "\n",
    "X = digits.data\n",
    "y = digits.target\n",
    "Xtrain,Xtest = np.split(X,[1000])\n",
    "ytrain,ytest = np.split(y,[1000])\n",
    "#Xtrain = X[:1000,:]\n",
    "#ytrain = y[:1000]\n",
    "#Xtest = X[1000:,:]\n",
    "#ytest = y[1000:]\n",
    "\n",
    "print(digits.DESCR)\n",
    "\n",
    "plt.gray();\n",
    "plt.matshow(digits.images[0]);\n",
    "plt.show();\n",
    "plt.matshow(digits.images[15]);\n",
    "plt.show();\n",
    "plt.matshow(digits.images[42]);\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "428a56b4-f981-4118-8ce6-9147576f53e4"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.utils import shuffle\n",
    "\n",
    "def shuffle_and_split(X,y,n):\n",
    "    X0,y0 = shuffle(X,y)\n",
    "    Xtrain,Xtest = np.split(X0,[n])\n",
    "    ytrain,ytest = np.split(y0,[n])\n",
    "    return Xtrain, ytrain, Xtest, ytest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "7aee1580-4a5f-461d-9bdf-9632ae1da643"
    }
   },
   "outputs": [],
   "source": [
    "print(Xtrain.shape)\n",
    "print(ytrain.shape)\n",
    "digits_svc = svm.SVC(kernel='rbf', gamma=1e-3)\n",
    "digits_svc.fit(Xtrain,ytrain)\n",
    "prediction = digits_svc.predict(Xtest)\n",
    "#print(\"Training error:\", np.sum(np.not_equal(prediction,ytest))/len(ytest))\n",
    "print(\"Generalization error:\", np.sum(np.not_equal(prediction,ytest))/len(ytest) )\n",
    "print(\"Generalization score:\", digits_svc.score(Xtest,ytest))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "55d61d0e-45a0-4442-9a27-6ac040e0c599"
    }
   },
   "outputs": [],
   "source": [
    "# Compute cross-validation score\n",
    "nb_trials = 20\n",
    "score = []\n",
    "for i in range(nb_trials):\n",
    "    Xtrain, ytrain, Xtest, ytest = shuffle_and_split(X,y,1000)\n",
    "    digits_svc = svm.SVC(kernel='rbf', gamma=1e-3)\n",
    "    digits_svc.fit(Xtrain,ytrain)\n",
    "    score += [digits_svc.score(Xtest,ytest)]\n",
    "    print('*',end='')\n",
    "print(\" done!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "b9942737-31dd-44c7-a08f-95f7a75a83f5"
    }
   },
   "outputs": [],
   "source": [
    "print(\"Average generalization score:\", np.mean(score))\n",
    "print(\"Standard deviation:\", np.std(score))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's identify which are the misclassified images (and find the confusion matrix by the way)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "51a0fedf-7c31-4527-a389-32d247fa89be"
    }
   },
   "outputs": [],
   "source": [
    "# Retrain\n",
    "Xtrain = X[:1000,:]\n",
    "ytrain = y[:1000]\n",
    "Xtest = X[1000:,:]\n",
    "ytest = y[1000:]\n",
    "digits_svc = svm.SVC(kernel='rbf', gamma=1e-3)\n",
    "digits_svc.fit(Xtrain,ytrain);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "c4d2ac07-1cf7-465f-9089-1b132459f980"
    }
   },
   "outputs": [],
   "source": [
    "# Examples\n",
    "N = 1053\n",
    "plt.matshow(digits.images[N]) \n",
    "plt.show() \n",
    "x = digits.data[N,:]\n",
    "print(\"prediction on image number\", N, \":\", digits_svc.predict([digits.data[N,:]]))\n",
    "print(\"correct label                :\", digits.target[N])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "c7c0ae00-4061-4d28-a9da-419c4a11f21a"
    }
   },
   "outputs": [],
   "source": [
    "# Find misclassified examples\n",
    "ypredict = digits_svc.predict(Xtest)\n",
    "misclass = np.not_equal(ypredict, ytest)\n",
    "Itest = digits.images[1000:,:]\n",
    "Xmisclass = Xtest[misclass,:]\n",
    "ymisclass = ytest[misclass]\n",
    "Imisclass = Itest[misclass,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Confusion matrix:\")\n",
    "print(confusion_matrix(ytest, ypredict))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "4f03a16d-a638-444e-881f-cbb26eac5685"
    }
   },
   "outputs": [],
   "source": [
    "# Display misclassified examples\n",
    "N = 1\n",
    "plt.matshow(Imisclass[N]) \n",
    "plt.show() \n",
    "print(\"prediction on image number\", N, \":\", digits_svc.predict([Xmisclass[N,:]]))\n",
    "print(\"correct label                :\", ymisclass[N])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"sec8\"></a>8.3 : ionosphere dataset\n",
    "\n",
    "This dataset contains data collected by a radar system in Goose Bay, Labrador. This system consists of a phased array of 16 high-frequency antennas designed to detect free electrons in the ionosphere. In general, there are two types of structures in the ionosphere: “Good” and “Bad”. These radars detected these structures and passed the signal. \n",
    "\n",
    "[Download from UCI](https://archive-beta.ics.uci.edu/ml/datasets/ionosphere)\n",
    "\n",
    "**Questions** :\n",
    "- How many variable are there ? Which are independent and which are dependent ?\n",
    "- Is it a good dataset to try classification on ?\n",
    "- Using SVM, try to obtain at least 64% of classification performance (top is 94%). Note : justify your result with k-fold cross validation !\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.read_csv(\"../data/ionosphere/ionosphere.csv\",header=None)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Aucun(e)",
  "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.11.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}