3 Kernel, Bayes and Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Assignment 3 - basic classifiers\n",
    "\n",
    "Math practice and coding application for main classifiers introduced in Chapter 3 of the Python machine learning book. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Weighting\n",
    "\n",
    "Note that this assignment is more difficult than the previous ones, and thus has a higher weighting 3 and longer duration (3 weeks). Each one of the previous two assignments has a weighting 1.\n",
    "\n",
    "Specifically, the first 3 assignments contribute to your continuous assessment as follows:\n",
    "\n",
    "Assignment weights: $w_1 = 1, w_2 = 1, w_3 = 3$\n",
    "\n",
    "Assignment grades: $g_1, g_2, g_3$\n",
    "\n",
    "Weighted average: $\\frac{1}{\\sum_i w_i} \\times \\sum_i \\left(w_i \\times g_i \\right)$\n",
    "\n",
    "Future assignments will be added analogously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# RBF kernel (20 points)\n",
    "\n",
    "Show that a Gaussian RBF kernel can be expressed as a dot product:\n",
    "$$\n",
    "K(\\mathbf{x}, \\mathbf{y}) \n",
    "= e^\\frac{-|\\mathbf{x} - \\mathbf{y}|^2}{2} \n",
    "= \\phi(\\mathbf{x})^T \\phi(\\mathbf{y})\n",
    "$$\n",
    "by spelling out the mapping function $\\phi$.\n",
    "\n",
    "For simplicity\n",
    "* you can assume both $\\mathbf{x}$ and $\\mathbf{y}$ are 2D vectors\n",
    "$\n",
    "x =\n",
    "\\begin{pmatrix}\n",
    "x_1 \\\\\n",
    "x_2\n",
    "\\end{pmatrix}\n",
    ", \\;\n",
    "y =\n",
    "\\begin{pmatrix}\n",
    "y_1 \\\\\n",
    "y_2\n",
    "\\end{pmatrix}\n",
    "$\n",
    "* we use a scalar unit variance here\n",
    "\n",
    "even though the proof can be extended for vectors $\\mathbf{x}$ $\\mathbf{y}$ and general covariance matrices.\n",
    "\n",
    "Hint: use Taylor series expansion of the exponential function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answer\n",
    "We denote $e^x$ as exp($x$). Since $\n",
    "\\mathbf x =\n",
    "\\begin{pmatrix}\n",
    "x_1 \\\\\n",
    "x_2\n",
    "\\end{pmatrix}\n",
    ", \\;\n",
    "\\mathbf y =\n",
    "\\begin{pmatrix}\n",
    "y_1 \\\\\n",
    "y_2\n",
    "\\end{pmatrix}\n",
    "$, we have\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "K(\\mathbf{x}, \\mathbf{y}) = \\text{exp}(\\frac{-||\\mathbf{x} - \\mathbf{y}||^2}{2}) & =  \\text{exp}(\\frac{-||(x_1-y_1, x_2-y_2)||^2}{2} )\\\\\n",
    "& = \\text{exp}(\\frac{-(x_1-y_1)^2-(x_2-y_2)^2}{2})\\\\\n",
    "& = \\text{exp}(\\frac{-{x_1}^2-{x_2}^2-{y_1}^2-{y_2}^2+2 x_1 y_1+2 x_2 y_2}{2}) \\\\\n",
    "& = \\text{exp}(\\frac{-||\\mathbf{x}||^2}{2}) \\text{ exp}(\\frac{-||\\mathbf{y}||^2}{2}) \\text{ exp}(x_1 y_1 + x_2 y_2)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "<br>By Taylor series of $f(x)$ on $a$, $e^x$ at $a=0$ can be expressed as $\\sum_{n=0}^{\\infty} \\frac{x^n}{n!} = 1 + x + \\frac{x^2}{2!} +$ ... for $x \\in \\mathbb{R}^1$. Therefore,<br><br>\n",
    "\n",
    "$$K(\\mathbf{x}, \\mathbf{y}) = \\text{exp}(\\frac{-||\\mathbf{x}||^2}{2}) \\text{ exp}(\\frac{-||\\mathbf{y}||^2}{2}) \\sum_{n=0}^{\\infty} \\frac{(x_1 y_1 + x_2 y_2)^n}{n!}$$\n",
    "\n",
    "<br>By binomial expansion, we have\n",
    "\n",
    "$$ (x_1 y_1 + x_2 y_2)^n = \\sum_{i = 0}^{n} \\binom{n}{i} (x_1 y_1)^{n-i} (x_2 y_2)^i = \\sum_{i = 0}^{n} \\sqrt{\\binom{n}{i}} (x_1^{n-i} x_2^i) \\sqrt{\\binom{n}{i}} (y_1^{n-i} y_2^i)$$\n",
    "\n",
    "<br>We let $\\xi_{n}(\\mathbf{x}) = \\xi_{n}(x_1, x_2) = \\left[\\sqrt{\\binom{n}{i}} (x_1^{n-i} x_2^i) \\right] = \\left[\\sqrt{\\binom{n}{0}} (x_1^{n} x_2^0), \\sqrt{\\binom{n}{1}} (x_1^{n-1} x_2^1), ..., \\sqrt{\\binom{n}{n}} (x_1^{0} x_2^n)\\right] \\in \\mathbb{R}^{\\text{n}}$. <br>\n",
    "\n",
    "<br>Hence we have <br>\n",
    "\n",
    "$$ \\left\\{\n",
    "\\begin{aligned}\n",
    "\\phi(\\mathbf{x}) & = \\text{exp}(\\frac{-||\\mathbf{x}||^2}{2}) \\left[1, \\frac{\\xi_{1}(\\mathbf{x})}{\\sqrt{1!}}, \\frac{\\xi_{2}(\\mathbf{x})}{\\sqrt{2!}} , \\frac{\\xi_{3}(\\mathbf{x})}{\\sqrt{3!}} , ... \\right]^T \\\\\n",
    "\\phi(\\mathbf{y}) & = \\text{exp}(\\frac{-||\\mathbf{y}||^2}{2}) \\left[1, \\frac{\\xi_{1}(\\mathbf{y})}{\\sqrt{1!}}, \\frac{\\xi_{2}(\\mathbf{y})}{\\sqrt{2!}} , \\frac{\\xi_{3}(\\mathbf{y})}{\\sqrt{3!}} , ... \\right]^T\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "<br>The mapping function is therefore $\\phi(\\mathbf{x}) = \\text{exp}(\\frac{-||\\mathbf{x}||^2}{2}) \\left[1, \\frac{\\xi_{1}(\\mathbf{x})}{\\sqrt{1!}}, \\frac{\\xi_{2}(\\mathbf{x})}{\\sqrt{2!}} , \\frac{\\xi_{3}(\\mathbf{x})}{\\sqrt{3!}} , ... \\right]^T$, where $\\xi_{n}(\\mathbf{x}) = \\left[\\sqrt{\\binom{n}{0}} (x_1^{n} x_2^0), \\sqrt{\\binom{n}{1}} (x_1^{n-1} x_2^1), ..., \\sqrt{\\binom{n}{n}} (x_1^{0} x_2^n)\\right]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Kernel SVM complexity (10 points)\n",
    "\n",
    "How would the complexity (in terms of number of parameters) of a trained kernel SVM change with the amount of training data, and why?\n",
    "Note that the answer may depend on the specific kernel used as well as the amount of training data.\n",
    "Consider specifically the following types of kernels $K(\\mathbf{x}, \\mathbf{y})$.\n",
    "* linear:\n",
    "$$\n",
    "K\\left(\\mathbf{x}, \\mathbf{y}\\right) = \\mathbf{x}^T \\mathbf{y}\n",
    "$$\n",
    "* polynomial with degree $q$:\n",
    "$$\n",
    "K\\left(\\mathbf{x}, \\mathbf{y}\\right) =\n",
    "(\\mathbf{x}^T\\mathbf{y} + 1)^q\n",
    "$$\n",
    "* RBF with distance function $D$:\n",
    "$$\n",
    "K\\left(\\mathbf{x}, \\mathbf{y} \\right) = e^{-\\frac{D\\left(\\mathbf{x}, \\mathbf{y} \\right)}{2s^2}}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answer\n",
    "\n",
    "For all examples, we assume $\\mathbf{x}, \\mathbf{y} \\in \\mathbb{R}^\\text{d}$.\n",
    "\n",
    "### Linear:\n",
    "For linear kernal, the mapping function $\\phi(\\mathbf{x}) = \\mathbf{x}$, which mapps $\\mathbb{R}^\\text{d}$ to $\\mathbb{R}^\\text{d}$, therefore the size of data is unchanged.<br>\n",
    "There are not explicit parameters, therefore the time cost increase linearly with the dimension of data, or the amount of data increase $n$ times, the time cost simply increase $O(n)$ time. Both changes in dimension or data amount will not afftect any parameters.<br>\n",
    "\n",
    "### Polynomial with degree $q$:\n",
    "For simplicity we write $1 = x_{d+1} y_{d+1}$. Then\n",
    "\n",
    "$$K\\left(\\mathbf{x}, \\mathbf{y}\\right) =(\\mathbf{x}^T\\mathbf{y} + 1)^q = (\\sum_{i=1}^{d+1} x_i y_i)^q = \\sum_{k_1 + k_2 + ... + k_{d+1} = q} \\binom{q}{k_1, k_2, ..., k_{d+1}} \\prod_{t=1}^{d+1} (x_t y_t)^{k_t} = \\sum_{\\sum_{i=1}^{d+1} k_i = q} \\frac{q!}{\\prod_{i=1}^{d+1} k_i!} \\prod_{t=1}^{d+1} (x_t y_t)^{k_t}$$\n",
    "\n",
    "by Multinomial theorem. Therefore the mapping function is\n",
    "\n",
    "$$\\phi(\\mathbf{x}) = \\left[\\sqrt{\\frac{q!}{\\prod_{i=1}^{d+1} k_i!}} \\prod_{t=1}^{d+1} (x_t)^{k_t}\\right]_{\\sum_{i=1}^{d+1} k_i = q}^T,$$\n",
    "\n",
    "which maps $\\mathbb{R}^\\text{d}$ to $\\mathbb{R}^\\binom{p+(d+1)-1}{(d+1)-1} = \\mathbb{R}^\\binom{p+d}{d} = \\mathbb{R}^\\frac{(p+d)!}{p! d!}$, computed using the stars and bars method. <br>\n",
    "\n",
    "* If $p=1$, only one useless dimension is added, where $x_{d+1} = 1$. In this case the actual dimension remains. <br>\n",
    "* If $p>2$, then the dimension increases from $d$ to $\\binom{p+d}{d}$, where actural dimension is $\\binom{p+d}{d} - 1$ since we always have a $x_{d+1}^q = 1$ term.<br><br>\n",
    "\n",
    "Now we consider the parameters.\n",
    "* For each entry in $K\\left(\\mathbf{x}, \\mathbf{y}\\right)$, we have a parameter $\\frac{q!}{\\prod_{i=1}^{d+1} k_i!} = \\binom{q}{k_1, k_2, ..., k_{d+1}}$, which takes $O(q \\prod_{t=1}^{d+1} k_t)$ to compute in brute force. Considering the dimension analysis we discuss above, the greater the dimension or the greater $q$ is, the more parameters and greater time complexity we will have in the kernal function.<br>\n",
    "* However, since $q$ and $k_i$ are identical for any set of input data, increasing amount of data will not change number of parameter to be calculated (because they only need to be calculated once), although multiplying them to each term of $x$ and $y$ takes constant time.<br>\n",
    "* If we do $\\mathbf{x}^T \\mathbf{y} + 1$ first and then do the power function, then the parameter analysis is the same as the linear function, except that we need an extra power $q$ after the $\\mathbf{x}^T \\mathbf{y} + 1$.<br>\n",
    "\n",
    "### RBF with distance function $D$:\n",
    "\n",
    "Assume $D(\\mathbf{x}, \\mathbf{y}) = \\omega(\\mathbf{x}) \\omega(\\mathbf{y})$. For $K(\\mathbf{x}, \\mathbf{y} ) = e^{-\\frac{D\\left(\\mathbf{x}, \\mathbf{y} \\right)}{2s^2}} = e^{-\\frac{\\omega(\\mathbf{x}) \\omega(\\mathbf{y})}{2s^2}} = e^{-\\frac{1}{2s^2} \\omega(\\mathbf{x}) \\omega(\\mathbf{y})}$, we have the mapping function as\n",
    "\n",
    "$$\n",
    "\\phi(\\mathbf{x}) = e^{-\\frac{1}{4s^2}} \\left[1, \\frac{\\omega(\\mathbf{x})}{\\sqrt {1!}}, \\frac{\\omega(\\mathbf{x})^2}{\\sqrt {2!}}, \\frac{\\omega(\\mathbf{x})^3}{\\sqrt {3!}}, ... \\right]^T,\n",
    "$$\n",
    "\n",
    "which maps $\\mathbb{R}^\\text{d}$ to $\\mathbb{R}^\\infty$. That is, RBF essentially projects the original vector to an infinite dimensional space.<br><br>\n",
    "\n",
    "Now we consider the parameters.\n",
    "* We first clarify that although RBF maps $\\mathbb{R}^\\text{d}$ to $\\mathbb{R}^\\infty$, the dimension actually used is determined by the explicit function $K(\\mathbf{x}, \\mathbf{y})$, because we don't have to separate the mapping function. Instead, we can just compute the kernal $K$ directly.<br>\n",
    "* As calculating exp($\\mathbf{x}$) is simply mapping to an exponential function, the main cost in terms of dimension of data is at the distance function $D(\\mathbf{x}, \\mathbf{y})$, which varies with different distance functions we choose. For example, if we choose simple metrics such as Taxicab distance and Euclidean distance, the cost is relatively small. However, if we choose complex metrics for some reason, then the time cost could be huge.<br>\n",
    "* Also, if all input data share the same set of parameters, then the parameters only need to be computed once, and applyed to each set of input data with constant time. However, if parameters change with different sets of input data in a specific kernal, then the number of parameters as well as time complexity also increase linearly with the increase of amount of data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Gaussian density Bayes (30 points)\n",
    "\n",
    "$$\n",
    "p\\left(\\Theta | \\mathbf{X}\\right)\n",
    "= \n",
    "\\frac{p\\left(\\mathbf{X} | \\Theta\\right) p\\left(\\Theta\\right)}{p\\left(\\mathbf{X}\\right)}\n",
    "$$\n",
    "\n",
    "Assume both the likelihood and prior have Gaussian distributions:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(\\mathbf{X} | \\Theta)\n",
    "&=\n",
    "\\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)\n",
    "\\\\\n",
    "p(\\Theta)\n",
    "&=\n",
    "\\frac{1}{\\sqrt{2\\pi}\\sigma_0} \\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Derive $\\Theta$ from the dataset $\\mathbf{X}$ via the following methods:\n",
    "\n",
    "### ML (maximum likelihood) estimation \n",
    "$$\n",
    "\\Theta_{ML} = argmax_{\\Theta} p(\\mathbf{X} | \\Theta)\n",
    "$$\n",
    "\n",
    "### MAP estimation\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Theta_{MAP} \n",
    "&= \n",
    "argmax_{\\Theta} p(\\Theta | \\mathbf{X})\n",
    "\\\\\n",
    "&=\n",
    "argmax_{\\Theta} p(\\mathbf{X} | \\Theta) p(\\Theta)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "### Bayes estimation\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Theta_{Bayes} \n",
    "&= \n",
    "E(\\Theta | \\mathbf{X})\n",
    "\\\\\n",
    "&= \n",
    "\\int \\Theta p(\\Theta | \\mathbf{X}) d\\Theta\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Answer\n",
    "\n",
    "### 1. ML (maximum likelihood) estimation \n",
    "\n",
    "To maximize $p(\\mathbf{X} | \\Theta)$, we set $\\nabla_\\Theta p(\\mathbf{X} | \\Theta) = \\nabla_\\Theta \\left(\\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)\\right) = 0$. <br>\n",
    "By Chain rule we get<br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\nabla_\\Theta p(\\mathbf{X} | \\Theta) & = \\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\frac{\\partial \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)}{\\partial \\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)} \\frac{\\partial \\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)}{\\partial \\Theta}\\\\\n",
    "0 & = \\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right) \\left( - \\frac{\\sum_{t=1}^N -2(\\mathbf{x}^{(t)} - \\Theta)}{2\\sigma^2} \\right) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "Note that $p(\\mathbf{X} | \\Theta) = \\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)$ is non-zero, because $e^y$ is always positive for $y \\in \\mathbb{R}$, and the constant $\\frac{1}{(2\\pi)^{N/2}\\sigma^N}$ is positive. <br>\n",
    "Then we have<br>\n",
    "$$\n",
    "\\begin{align}\n",
    "0 & = \\frac{\\sum_{t=1}^N 2(\\mathbf{x}^{(t)} - \\Theta)}{2\\sigma^2} = \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)} - \\Theta}{\\sigma^2} \\\\\n",
    "0 & = \\frac{(\\sum_{t=1}^N \\mathbf{x}^{(t)}) - N\\Theta}{\\sigma^2} \\\\\n",
    "N\\Theta & = \\sum_{t=1}^N \\mathbf{x}^{(t)} \\\\\n",
    "\\Theta & = \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}\n",
    "\\end{align}\n",
    "$$\n",
    "Hence $$\\Theta_{ML} = argmax_{\\Theta} p(\\mathbf{X} | \\Theta) = \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}.$$\n",
    "\n",
    "### 2. MAP estimation\n",
    "To maximize $p(\\mathbf{X} | \\Theta) p(\\Theta)$, we set $\\nabla_\\Theta (p(\\mathbf{X} | \\Theta) p(\\Theta)) = p(\\Theta)\\nabla_\\Theta p(\\mathbf{X} | \\Theta) + p(\\mathbf{X} | \\Theta)\\nabla_\\Theta p(\\Theta) = 0$. <br><br>\n",
    "We get $p(\\Theta)\\nabla_\\Theta p(\\mathbf{X} | \\Theta) = - p(\\mathbf{X} | \\Theta)\\nabla_\\Theta p(\\Theta) $ and therefore <br><br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\frac{1}{\\sqrt{2\\pi}\\sigma_0} \\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right) \\nabla_\\Theta \\left(\\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)\\right) & \\\\\n",
    "= - \\frac{1}{(2\\pi)^{N/2}\\sigma^N} \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right) \\nabla_\\Theta & \\left(\\frac{1}{\\sqrt{2\\pi}\\sigma_0} \\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right)\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "<br><br>By Removing the constants we get<br><br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right) \\nabla_\\Theta \\left(\\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)\\right) & = - \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right) \\nabla_\\Theta \\left(\\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right)\\right) \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "<br>By Chain rule we get<br><br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right) \\frac{\\partial \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)}{\\partial \\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)} \\frac{\\partial \\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right)}{\\partial \\Theta} & = - \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right) \\frac{\\partial \\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right)}{\\partial \\left(-\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2}\\right)} \\frac{\\partial \\left(-\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2}\\right)}{\\partial \\Theta}\\\\\n",
    "\\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right) \\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right) \\left( - \\frac{\\sum_{t=1}^N -2(\\mathbf{x}^{(t)} - \\Theta)}{2\\sigma^2} \\right) & = -\\exp\\left(-\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)^2}{2\\sigma^2}\\right) \\exp\\left( -\\frac{(\\Theta - \\mu_0)^2}{2\\sigma_0^2} \\right) \\left(-\\frac{2(\\Theta - \\mu_0)}{2\\sigma_0^2}\\right)\n",
    "\\end{align}\n",
    "$$\n",
    "<br>Since $e^y$ is always positive for $y \\in \\mathbb{R}$, the first two terms of both sides are non-zero. Dividing them on both sides we have<br>\n",
    "$$\n",
    "\\begin{align}\n",
    "- \\frac{\\sum_{t=1}^N -2(\\mathbf{x}^{(t)} - \\Theta)}{2\\sigma^2} & = \\frac{2(\\Theta - \\mu_0)}{2\\sigma_0^2}\\\\\n",
    "\\frac{\\sum_{t=1}^N (\\mathbf{x}^{(t)} - \\Theta)}{\\sigma^2} & = \\frac{\\Theta -\\mu_0}{\\sigma_0^2}\\\\\n",
    "\\frac{(\\sum_{t=1}^N \\mathbf{x}^{(t)}) - N\\Theta}{\\sigma^2} & = \\frac{\\Theta -\\mu_0}{\\sigma_0^2}\\\\\n",
    "(\\sum_{t=1}^N \\mathbf{x}^{(t)})\\sigma_0^2 - N\\Theta\\sigma_0^2 & = \\sigma^2 \\Theta - \\sigma^2 \\mu_0 \\\\\n",
    "(\\sum_{t=1}^N \\mathbf{x}^{(t)})\\sigma_0^2 + \\sigma^2 \\mu_0 & = (N\\sigma_0^2 + \\sigma^2)\\Theta \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "<br>Hence <br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Theta_{MAP} & = argmax_{\\Theta} p(\\Theta | \\mathbf{X}) \\\\\n",
    "& = argmax_{\\Theta} p(\\mathbf{X} | \\Theta) p(\\Theta) \\\\\n",
    "& = \\frac{(\\sum_{t=1}^N \\mathbf{x}^{(t)})\\sigma_0^2 + \\sigma^2 \\mu_0}{N\\sigma_0^2 + \\sigma^2} \\\\\n",
    "\\end{align}\n",
    "$$<br>\n",
    "Furthurmore, since $\\Theta_{ML} = \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}$, we have<br>\n",
    "$$\\Theta_{MAP} = \\frac{(\\sum_{t=1}^N \\mathbf{x}^{(t)})\\sigma_0^2 + \\sigma^2 \\mu_0}{N\\sigma_0^2 + \\sigma^2}\n",
    "= \\frac{N\\Theta_{ML}\\sigma_0^2 + \\sigma^2 \\mu_0}{N\\sigma_0^2 + \\sigma^2} = \\frac{N/\\sigma^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\Theta_{ML} +\\frac{1/\\sigma_0^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\mu_0$$\n",
    "\n",
    "### 3. Bayes estimation\n",
    "\n",
    "For $\\Theta_{Bayes} = E(\\Theta | \\mathbf{X}) = \\int \\Theta p(\\Theta | \\mathbf{X}) d\\Theta \\\\$, since $p(\\Theta | \\mathbf{X}) = \\frac{p(\\mathbf{X}| \\Theta) p(\\Theta)}{p(\\mathbf{X})}$ and $p(\\mathbf{X})$ is a constant for given $\\mathbf{X}$, our interest is in $p(\\mathbf{X}| \\Theta) p(\\Theta)$. We denote $\\phi(x, \\mu, \\sigma^2)$ as the normal distribution with input $x$, mean $\\mu$ and standard deviation $\\sigma$. Then  $$p(\\mathbf{X}| \\Theta) p(\\Theta) = \\phi(\\Theta, \\mu_0, \\sigma_0^2) \\prod_{i=1}^N \\phi(\\Theta, \\mathbf{x}^{(i)}, \\sigma^2).$$\n",
    "\n",
    "Notice that $$\\phi(x, \\mu_1, \\sigma_1^2) \\phi(x, \\mu_2, \\sigma_2^2) = \\phi(\\mu_1, \\mu_2, \\sigma_1^2 + \\sigma_2^2) \\phi(x, \\mu_i, \\sigma_i^2)$$ where $$\\mu_i = \\frac{1 / \\sigma_1^2}{1/ \\sigma_1^2 + 1/ \\sigma_2^2}\\mu_1 + \\frac{1 / \\sigma_2^2}{1/ \\sigma_1^2 + 1/ \\sigma_2^2}\\mu_2 \\text{ and } \\sigma_i^2 = \\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2 + \\sigma_2^2}.$$<br>\n",
    "\n",
    "We will prove the formula later on. Using this formuma, we have $$\\phi(\\Theta, \\mathbf{x}^{(a)}, \\sigma^2) \\phi(\\Theta, \\mathbf{x}^{(b)}, \\sigma^2) = \\phi(\\mathbf{x}^{(a)}, \\mathbf{x}^{(b)}, 2\\sigma^2) \\phi(\\Theta, \\frac{\\mathbf{x}^{(a)}+\\mathbf{x}^{(b)}}{2}, \\frac{\\sigma^2}{2}) = C_0 \\phi(\\Theta, \\frac{\\mathbf{x}^{(a)}+\\mathbf{x}^{(b)}}{2}, \\frac{\\sigma^2}{2}),$$ where $C_0$ is some constant since all variables of $\\phi(\\mathbf{x}^{(a)}, \\mathbf{x}^{(b)}, 2\\sigma^2)$ are set. Following similar steps we get\n",
    "\n",
    "$$\\prod_{i=1}^N \\phi(\\Theta, \\mathbf{x}^{(i)}, \\sigma^2) = C_1 \\phi(\\Theta,  \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}, \\frac{\\sigma^2}{N}),$$ where $C_1$ is some constant.<br>\n",
    "\n",
    "Hence, $$p(\\mathbf{X}| \\Theta) p(\\Theta) = C_1 \\phi(\\Theta, \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}, \\frac{\\sigma^2}{N}) \\phi(\\Theta, \\mu_0, \\sigma_0^2) =  C_2 \\phi(\\Theta, \\mu_\\text{new}, \\sigma_\\text{new}^2),$$ where $C_2$ is some constant and by the formula, $\\mu_\\text{new} = \\frac{N/\\sigma^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\Theta_{ML} +\\frac{1/\\sigma_0^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\mu_0 $, where $\\Theta_{ML} = \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}.$ <br>\n",
    "\n",
    "Notice that for a given normal distribution, multiplying the probability density function by a constant will not change its mean value. Therefore the expectation of $p(\\Theta | \\mathbf{X})$ is exactly the expectation of the non-constant normal distribution part. Hence,\n",
    "$$\n",
    "E(\\Theta | \\mathbf{X}) = \\mu_\\text{new} = \\frac{N/\\sigma^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\Theta_{ML} +\\frac{1/\\sigma_0^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\mu_0 = \\Theta_{MAP},\n",
    "$$\n",
    "where $\\Theta_{ML} = \\frac{\\sum_{t=1}^N \\mathbf{x}^{(t)}}{N}.$ <br>\n",
    "\n",
    "Finally we want to prove the formula $\\phi(x, \\mu_1, \\sigma_1^2) \\phi(x, \\mu_2, \\sigma_2^2) = \\phi(\\mu_1, \\mu_2, \\sigma_1^2 + \\sigma_2^2) \\phi(x, \\mu_i, \\sigma_i^2)$:<br><br>\n",
    "$$\n",
    "\\begin{align}\n",
    "\\phi(x, \\mu_1, \\sigma_1^2) \\phi(x, \\mu_2, \\sigma_2^2) & = \\frac{1}{\\sqrt{2\\pi}\\sigma_1} \\exp\\left( -\\frac{(x - \\mu_1)^2}{2\\sigma_1^2} \\right) \\frac{1}{\\sqrt{2\\pi}\\sigma_2} \\exp\\left( -\\frac{(x - \\mu_2)^2}{2\\sigma_2^2} \\right) \\\\\n",
    "& =  \\frac{1}{2\\pi \\sigma_1 \\sigma_2} \\exp\\left( -\\frac{(x - \\mu_1)^2}{2\\sigma_1^2} - \\frac{(x - \\mu_2)^2}{2\\sigma_2^2}\\right) \\\\\n",
    "& = \\frac{1}{2\\pi \\sigma_1 \\sigma_2} \\exp\\left( -\\frac{(\\sigma_1^2+\\sigma_2^2) x^2 -2(\\mu_1 \\sigma_2^2 +\\mu_2 \\sigma_1^2)x +(\\mu_1^2 \\sigma_2^2+\\mu_2^2 \\sigma_1^2)}{2\\sigma_1^2 \\sigma_2^2}    \\right)\\\\\n",
    "& = \\frac{1}{2\\pi \\sigma_1 \\sigma_2} \\exp\\left( -\\frac{x^2 -2\\frac{\\mu_1 \\sigma_2^2 +\\mu_2 \\sigma_1^2}{\\sigma_1^2+\\sigma_2^2}x + \\frac{\\mu_1^2 \\sigma_2^2+\\mu_2^2 \\sigma_1^2}{\\sigma_1^2+\\sigma_2^2}}{2\\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}}    \\right) \\\\\n",
    "& = \\frac{1}{2\\pi \\sigma_1 \\sigma_2} \\exp\\left( -\\frac{x^2 -2\\frac{\\mu_1 \\sigma_2^2 +\\mu_2 \\sigma_1^2}{\\sigma_1^2+\\sigma_2^2}x + \\frac{\\mu_1^2 \\sigma_2^4+\\mu_2^2 \\sigma_1^4 + 2\\mu_1 \\sigma_2^2 \\mu_2 \\sigma_1^2}{(\\sigma_1^2+\\sigma_2^2)^2}}{2\\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}} \\right) \\\\\n",
    "& \\times \\exp\\left(\\frac{ - (\\sigma_1^2+\\sigma_2^2)(\\mu_1^2 \\sigma_2^2+\\mu_2^2 \\sigma_1^2) +  (\\mu_1^2 \\sigma_2^4+\\mu_2^2 \\sigma_1^4 + 2\\mu_1 \\sigma_2^2 \\mu_2 \\sigma_1^2) }{2 \\sigma_1^2 \\sigma_2^2 (\\sigma_1^2+\\sigma_2^2)}\\right) \\\\\n",
    "& = \\frac{1}{2\\pi \\sqrt{\\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}} \\sqrt{\\sigma_1^2+\\sigma_2^2}} \\exp\\left( -\\frac{(x -\\frac{\\mu_1 \\sigma_2^2 +\\mu_2 \\sigma_1^2}{\\sigma_1^2+\\sigma_2^2})^2}{2\\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}} \\right) \\exp\\left( -\\frac{(\\mu_1 - \\mu_2)^2}{2 \\sigma_1^2 + \\sigma_2^2} \\right)\\\\\n",
    "& = \\frac{1}{\\sqrt{2\\pi} \\sqrt{\\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}}} \\exp\\left( -\\frac{(x -\\frac{\\mu_1 \\sigma_2^2 +\\mu_2 \\sigma_1^2}{\\sigma_1^2+\\sigma_2^2})^2}{2\\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}} \\right) \\frac{1}{\\sqrt{2\\pi} \\sqrt{\\sigma_1^2+\\sigma_2^2}} \\exp\\left( -\\frac{(\\mu_1 - \\mu_2)^2}{2 \\sigma_1^2 + \\sigma_2^2} \\right)\\\\\n",
    "& = \\phi(x, \\mu_i, \\sigma_i^2) \\phi(\\mu_1, \\mu_2, \\sigma_1^2 + \\sigma_2^2)\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "where $$\\mu_i = \\frac{1 / \\sigma_1^2}{1/ \\sigma_1^2 + 1/ \\sigma_2^2}\\mu_1 + \\frac{1 / \\sigma_2^2}{1/ \\sigma_1^2 + 1/ \\sigma_2^2}\\mu_2 \\text{ and } \\sigma_i^2 = \\frac{\\sigma_1^2 \\sigma_2^2}{\\sigma_1^2 + \\sigma_2^2}.$$<br>\n",
    "Hence we complete the proof and we validate that\n",
    "$$\\Theta_{Bayes} = \\mu_\\text{new} = \\frac{N/\\sigma^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\Theta_{ML} +\\frac{1/\\sigma_0^2}{N/\\sigma^2 + 1/\\sigma_0^2} \\mu_0.$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Hand-written digit classification (40 points)\n",
    "\n",
    "In the textbook sample code we applied different scikit-learn classifers for the Iris data set.\n",
    "\n",
    "In this exercise, we will apply the same set of classifiers over a different data set: hand-written digits.\n",
    "Please write down the code for different classifiers, choose their hyper-parameters, and compare their performance via the accuracy score as in the Iris dataset.\n",
    "Which classifier(s) perform(s) the best and worst, and why?\n",
    "\n",
    "The classifiers include:\n",
    "* perceptron\n",
    "* logistic regression\n",
    "* SVM\n",
    "* decision tree\n",
    "* random forest\n",
    "* KNN\n",
    "* naive Bayes\n",
    "\n",
    "The dataset is available as part of scikit learn, as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "last updated: 2016-10-17 \n",
      "\n",
      "CPython 3.5.2\n",
      "IPython 4.2.0\n",
      "\n",
      "numpy 1.11.1\n",
      "pandas 0.18.1\n",
      "matplotlib 1.5.1\n",
      "scipy 0.17.1\n",
      "sklearn 0.18\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -a '' -u -d -v -p numpy,pandas,matplotlib,scipy,sklearn\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Added version check for recent scikit-learn 0.18 checks\n",
    "from distutils.version import LooseVersion as Version\n",
    "from sklearn import __version__ as sklearn_version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1797, 64)\n",
      "(1797,)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.datasets import load_digits\n",
    "digits = load_digits()\n",
    "\n",
    "X = digits.data # training data\n",
    "y = digits.target # training label\n",
    "\n",
    "print(X.shape)\n",
    "print(y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Visualize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEZCAYAAADCJLEQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuQFfWVB/DvV3CFqMjD0kWMVMASjRI1PFbRBcux4lqa\nGpDgKxVxWAtL1FSC2SRbimsUrLgllGEFSpMFcRd5LSJl4TOUqwmPZdQki1liBAwqGFEYBFESlbN/\n3IuZ0L/fTN87ffv36zvfTxUVOXOm769Peu6hb5/ppplBREQkNoeFXoCIiIiLGpSIiERJDUpERKKk\nBiUiIlFSgxIRkSipQYmISJSCNSiS80jeVf7v80luTPl9qXM7E9UzW6pn9lTTbHWGekZxBmVmvzSz\n06rJJfkGyQvb+h6SDSQ3kvyQ5CqSJ3V0zTGrZT1JHk5yaTnvAMmRWaw5ZjWu59+RfJbkTpLvklxM\n8m+zWHfMalzT00g2k9xVruuzJFO9VlHV+j20Ve4d5Z/7VPkdFUWDqiWSfQAsA3AbgN4AXgawOOii\niu8XAL4J4J3QC6kDvQA8CKB/+c+HAOYFXVHxbQNwhZn1BnAsgCcALAq7pOIjOQDANwBsz+s1c2tQ\nJM8m+TLJD0guAtCt1ddGkXyr1d+/SvKVcu4Skotancp+nkvyEQAnAXiC5B6S33O89OUAXjWzx8zs\nzwDuBHAmyVNqt7e1F6qeZvaJmc00szUADtR6P/MSsJ5Pm9kyM/vQzPYDeADAiBrvbi4C1nSPmb1R\n/msXlI7TgbXb03wEfA89aBaA7wP4pBb755JLgyJ5OIDlAOajdBazFMDYQ9KsVe5jAOaWcxcCGOPK\nNbNrAbwJ4DIz62Fm9zle/nQAv/n8G80+ArCpHC+kwPWsO5HVcxSA31a3J/GIoaYkWwB8BOAnAKZ1\ncJeCCl1PkuMA7DezpzPZoZS65vQ65wDoamYzy39fRrLZk3sugC5m9kD578tJrm9n+2zja0cB2HFI\nbA+Ao9vZZsxC1rMeRVFPkl8BMAXA19PkRy54Tc2sF8nuAMaj9CZcZMHqSfIolBp8QyULzkJeDeoE\nlD4Xbm2rJ7evI/ctV2JKHwLocUjsGAB7O7DN0ELWsx4FryfJkwE8CeCW8senRRe8pgBgZh+TfBDA\neyRPNbP3s9huACHreSeAR8ws9/eNvK5BvQOg3yEx3ySdK/eLbWy7vdux/xbAWQf/QvJIlD6PLvLH\nKCHrWY+C1pNkfwDPAfiRmT3aXn5BxHSMdgHwBcdrFEnIejYA+DbJd0i+U97WEpL/1M73dVheDWot\ngE9J3kKyK8nLAQxvI/czkjeR7EKysY1cAPgjgAFtfH05gNNJjiF5BIB/AfBrM/t9FfsRi5D1BMm/\nIXnwAu0R5boWWbB6kuwHYBWAfzOzn1a5/hiFrOlFJM8ieRjJHgBmANgFoBC/++MR8mf+QgBnADiz\n/Gc7gIkoDU3UVC4Nysw+QWmargnATgDjUBr9biv3egAtAK5BaUz0T57N/xjAFJZ+52GyY3vvo3Qx\n8R6UDtKhAK7qyP6EFrKeZa8B2IfSxw5PA/iIBf7dssD1/EcAXwJwZ3mKai/JPR3ZnxgErmlPlAYD\ndgN4HaX6/kN5ireQAr+HtpjZjoN/AHwKYHd54KymWIQHFpJcB2COmc0PvZZ6oHpmS/XMnmqaraLW\nM8pf1CU5kuTx5dPT8QAGo/QvdamC6pkt1TN7qmm26qWeeU3xVWoQgCUoXdjcAmCsmb0bdkmFpnpm\nS/XMnmqarbqoZyE+4hMRkc6nzTMokupeKZlZu784qHqml6aegGqaluqZLdUze66atvsRXx5nWGvX\nrnXGJ0yYkIhdfvnlztwpU6Y44926dXPGs0Smv/FCyDPWMWMOvdsJsGPHoTfZKLn//vud8WHDhmW6\nJpdK6gmEq+lrr73mjJ977rnO+KhRoxKx5cuXZ7omlxjrOX9+8lr9dddd58w99dRTnfFf/epXiVhs\nP+9A2J/5/fv3J2KTJk1y5s6dO7fWy/Hy1TTKIQkRERE1KBERiZIalIiIRCmKMXPXtSYA+N3vfpeI\n7dq1y5nbvXt3Z3zNGvd9N33XCepZr169ErHHH3/cmfvMM88443lcg4rRtm2H3nvTf23EVWcA2LBh\nQ6ZrKoLp06c74z/72c8SsZUrVzpzL730Umd8y5YtidiXv/zlClZX/1asWJGIDR06NMBKqqMzKBER\niZIalIiIREkNSkREoqQGJSIiUVKDEhGRKOU6xffWW+4nBrum9QD3xJ5vQso33dcZp/hcE2eAf2LP\npZ7rUw3XNNSIESOcud/85jed8ZtuuinTNRWBb0LXVYuzzz7bmeubltTE3l+47hgBADNnzkzE7rrr\nLmfu7t27K3rNnj17VpRfDZ1BiYhIlNSgREQkSmpQIiISJTUoERGJUq5DEnv37nXGL7jgAmfcNxDh\nMnz48GqWVGiLFy92xm+88UZnvKWlJfW2hwwZUtWa6pXrYv+gQYOcuePGjXPGm5qaMl1TEfh+hl3H\nom9Y6oorrnDGXYMBeTxuI0auIR4A2LhxYyLW0NDgzJ06daoz3rt3b2fc99iOLOkMSkREoqQGJSIi\nUVKDEhGRKKlBiYhIlNSgREQkSrlO8X3wwQfO+GWXXdbhbftudeSbQKkHV155pTPe2NjojPse6uiy\nb98+ZzyP25uE5LtlzNy5cxOxBQsWVLTt2bNnV7WmeuSa7vv444+duZdccknq+FNPPeXMrZfpvubm\nZmf8qquucsYnT56cettTpkxxxn/+85+n3kbWdAYlIiJRUoMSEZEoqUGJiEiU1KBERCRKalAiIhKl\nXKf4jjnmGGd8/fr1qbfhm7LyPZjwuuuuS71t+QvffdH69euX80rydd999znjvgknF9/xXC+TZLXi\nq49vMu+73/1uIjZr1ixn7q233lr9wiLSo0cPZ9x3z8MZM2YkYuvWravoNc8777yK8rOkMygREYmS\nGpSIiERJDUpERKKkBiUiIlFSgxIRkSjlOsXXt29fZ3zVqlXO+Nq1axOxRx55pKLXHD9+fEX50rn5\nnnrrmiTzTY76nu7s2rbv6cfDhg3zLbEuTJ8+PRHz3XPPdw/PpUuXJmI33HBDxxYWOd9TnH33It22\nbVsiNnjwYGeu7759IadPdQYlIiJRUoMSEZEoqUGJiEiU1KBERCRKuQ5J+G7H4Rt8mDBhQiJ2wQUX\nOHOff/75qtdVb3wXNV0X6efNm+fMffLJJ53xhoaG6hdWAL5bOa1evToRc12ABvy3RXLVesCAAc7c\neh+SOPbYYxOxsWPHVrQN10DEtGnTql5TPTryyCMTsZaWFmfuxIkTa72ciukMSkREoqQGJSIiUVKD\nEhGRKKlBiYhIlNSgREQkSjQz/xdJ/xflr5gZ28tRPdNLU09ANU1L9cyW6pk9V03bbFAiIiKh6CM+\nERGJkhqUiIhESQ1KRESipAYlIiJRUoMSEZEoqUGJiEiU1KBERCRKalAiIhIlNSgREYlSsAZFch7J\nu8r/fT7JjSm/L3VuZ6J6Zkv1zJ5qmq3OUM8ozqDM7Jdmdlo1uSTfIHmhL59kf5IHSO4hubf8v7dl\nse5Y1bKe5ZzuJGeTfI9kC8n/7uCSo1bj4/OaVsflHpL7ysfr2VmsPVY5HKNXkPw/kh+QfJVkY0fX\nHLMc6nk9ydfLx+iTJPt2dM1pRNGgcmAAjjGzo82sh5npudAd81MAPQEMAtAbwHfDLqe4zOzRVsdl\nDwCTAGw2s1+FXltRkTwBwH8A+I6ZHQPg+wAeJZl8zry0i+QFAKYB+DpKP+9/ALAwj9fOrUGRPJvk\ny+V/0SwC0K3V10aRfKvV379K8pVy7hKSi1qdyn6eS/IRACcBeKLc2b/ne3nUWTMOVU+SgwBcBmCi\nme2yksK/mQY+PlsbD+CRTHcukIA1PRFAi5k9CwBm9iSAfQAG1mxncxCwnpcCWGpmvzOzTwHcDWAk\nyS/VcHcB5PSmTfJwAMsBzEepAy8FMPaQNGuV+xiAueXchQDGuHLN7FoAbwK4rPwv0Ps8SzAAfyD5\nJsm5JPt0fK/CCVzP4QC2Arir/BHfb0hensmOBRLB8XlwHf0B/D3qoEEFrulLADaSvIzkYSRHA9gP\n4H+z2LcQYjlGyw72jTMq35PK5HVWcQ6ArmY208w+M7NlAJo9uecC6GJmD5RzlwNY387223o2y/sA\nhgHoD2AIgKMBLKhs+dEJWc8TAQwG0AKgL4BbAMwvn1kVVch6tnYtgF+Y2daU+TELVlMzO4DSR3wL\nAfwJwH8CuMHMPq54L+IR8hh9GsA4kmeQ7A7gDgAHAHyhwn2oWF4N6gQA2w6J+X4I+zpy33IlpmFm\n+8zsFTM7YGbvAbgZwNdIHlntNiMQrJ4APgbwZwBTzexTM3sRwPMAvtaBbYYWsp6tfQvAwxltK7Rg\nNSV5EYB/BTDSzA4HcAGAfyf5lWq3GYGQ76GrANyJ0lnZlvKfvQDernabaeXVoN4B0O+Q2EkV5H6x\njW1X88RFQ7GvSYWs58GPSVr/i6voT70MfnySPA+lN5ZlafILIGRNzwTwwsFro2b2EoD/AXBRO98X\ns6DHqJnNMbNTzKwvSo2qK4BX2/u+jsrrTXotgE9J3kKya/maxfA2cj8jeRPJLiyNh/pyAeCPAAb4\nvkhyOMlTWNIHwE8APG9me6vclxgEqyeAF1H6zPqfy9s7D6V/oT5T8V7EI2Q9DxoPYJmZ7ato5fEK\nWdNmAOeTPBMoDRcAOB8FvgaFsO+hR5A8vfzfJwF4CMD9ZvZBVXtSgVwalJl9AuByAE0AdgIYB8+/\nFFvlXo/SdY5rADyB0mfJLj8GMIXkLpKTHV8fgNJnqHtQOkD3l7dZWCHrWZ7iaURpsmc3gAcBfMvM\nft+RfQop8PEJkkcA+Abq5+O90MfoiwB+BOC/SH6A0kDBNDP7eYd2KqDAx2g3lMb09wJYB2A1Steh\nao5m8X86Q3IdgDlmNj/0WuqB6pkt1TN7qmm2ilrPKK/DkBxJ8vjy6el4lKbGng69rqJSPbOlemZP\nNc1WvdSza+gFeAwCsASlMcYtAMaa2bthl1Roqme2VM/sqabZqot6FuIjPhER6XzaPIMiqe6Vkpm1\n+8uYqmd6aeoJqKZpqZ7ZUj2z56ppux/x5XGGNWbMoXfhKBkwIDn5OH369Fovp2Jk2hsF5FNPH1ed\nd+zY4cxdvXp1rZfjVUk9gXxqunjx4kRs586dztwFC9w3KlmzZk0i1qtXL2fu9u3bE7GuXbuia9fK\nP5WPsZ5Tp05NxB5++GFn7uTJzuFHTJgwIRHr1q2bIzNbMdbTVQsAaGlpScSWL19e6+VUzFfTKIck\nRERE1KBERCRKalAiIhKlNqf4SFoen5+efPLJzvjmzZtTb2PgQPejXjZt2lTVmipBMvWQRB71bG52\n3+R4+PDk3U5mzZrlzJ00aVKma6pE2nqWc3OpqesalM9ZZ53ljN97772JmOsaAZDtdYIY6+m6Hrph\nw4aKtjF48OBELI/rKyHruXv3bmfcdy2zEiNGjHDG87ge7aupzqBERCRKalAiIhIlNSgREYmSGpSI\niEQpinvxHX/88c64a0jCdzGwsbHRGd+/f78znscv9IXyne98J3Wur27y16688srUubNnz3bGX3vt\ntURs1apVVa+pyIYMGZKIuX4xH/D/cn7v3r0TMVeNAWDQoEEVrC5e+/ZV9riw0aNHJ2K+Oq9YsaKq\nNdWSzqBERCRKalAiIhIlNSgREYmSGpSIiERJDUpERKIUxRSfb8LG9XgC361hXLfxAep7Ws/n3Xfd\nD8503cqkX79+tV5OofimwCqZtrv99ttT5/puI9PQ0JB6G0XU1NSUiJ144onO3C1btjjjrik+30Rw\nvejTp09F+QsXLkzErr76amfurl27qlpTLekMSkREoqQGJSIiUVKDEhGRKKlBiYhIlNSgREQkSlFM\n8c2dO9cZ/8EPfpCI/frXv3bmXnXVVRW9ZiX3Visa3zSO6wFvvgfxXXzxxc54z549q19YAfimwF56\n6aVE7PHHH69o22vXrk3E6uUecZX68MMPU+f66uya6K3349M3lex72GD37t0TsbvvvtuZ+8ILLzjj\nvock5lFrnUGJiEiU1KBERCRKalAiIhIlNSgREYlSFEMSPllcQH799dczWEmxnHbaac6462Lzjh07\nnLm+oZO3337bGa+XWyb5Lvy6BnnmzZvnzF2/fr0z3hkHIrZt2+aMn3rqqYnYrFmznLmuB5cCwKWX\nXpqIrVy50plb78MTvltmuepf6c/q5MmTnXHfcFuWdAYlIiJRUoMSEZEoqUGJiEiU1KBERCRKalAi\nIhKlKKb4mpubnfEePXokYj/84Q8r2va4ceOqWlORffvb33bGXQ+A9E2Wbdy40RlfsWKFMz5p0qSU\nqyumqVOnJmK9evVy5rpuKdVZ+R6w56rdhAkTnLk7d+50xl0POHz00UedufV+fPq4JvZcxzIAzJgx\nwxl33aIrLzqDEhGRKKlBiYhIlNSgREQkSmpQIiISpSiGJJ555hlnfMqUKam34bsdR2e8vUxjY6Mz\n7noOjO/C6OjRoyvadr176qmnEjHfcet7Zk9n5KuF6/hyPbsI8A+jNDU1JWK+QYt65xt8ePnllxMx\n3+3NNmzY4IyHvI2ZzqBERCRKalAiIhIlNSgREYmSGpSIiERJDUpERKJEM/N/kfR/Uf6KmbG9HNUz\nvTT1BFTTtFTPbKme2XPVtM0GJSIiEoo+4hMRkSipQYmISJTUoEREJEpqUCIiEiU1KBERiZIalIiI\nREkNSkREoqQGJSIiUVKDEhGRKAVrUCTnkbyr/N/nk9yY8vtS53Ymqme2VM/sqabZ6gz1jOIMysx+\naWanVZNL8g2SF7b1PSQbSG4k+SHJVSRP6uiaY1bLepI8nOTSct4BkiOzWHPMalzPvyP5LMmdJN8l\nuZjk32ax7pjVuKankWwmuatc12dJpnqtoqr1e2ir3DvKP/ep8jsqigZVSyT7AFgG4DYAvQG8DGBx\n0EUV3y8AfBPAO6EXUgd6AXgQQP/ynw8BzAu6ouLbBuAKM+sN4FgATwBYFHZJxUdyAIBvANie12vm\n1qBInk3yZZIfkFwEoFurr40i+Varv3+V5Cvl3CUkF7U6lf08l+QjAE4C8ATJPSS/53jpywG8amaP\nmdmfAdwJ4EySp9Rub2svVD3N7BMzm2lmawAcqPV+5iVgPZ82s2Vm9qGZ7QfwAIARNd7dXASs6R4z\ne6P81y4oHacDa7en+Qj4HnrQLADfB/BJLfbPJZcGRfJwAMsBzEfpLGYpgLGHpFmr3McAzC3nLgQw\nxpVrZtcCeBPAZWbWw8zuc7z86QB+8/k3mn0EYFM5XkiB61l3IqvnKAC/rW5P4hFDTUm2APgIwE8A\nTOvgLgUVup4kxwHYb2ZPZ7JDKXXN6XXOAdDVzGaW/76MZLMn91wAXczsgfLfl5Nc387223o2y1EA\ndhwS2wPg6Ha2GbOQ9axHUdST5FcATAHw9TT5kQteUzPrRbI7gPEovQkXWbB6kjwKpQbfUMmCs5BX\ngzoBpc+FW9vqye3ryH3LlZjShwB6HBI7BsDeDmwztJD1rEfB60nyZABPAril/PFp0QWvKQCY2cck\nHwTwHslTzez9LLYbQMh63gngETPL/X0jr2tQ7wDod0jMN0nnyv1iG9tu74mLvwVw1sG/kDwSpc+j\ni/wxSsh61qOg9STZH8BzAH5kZo+2l18QMR2jXQB8wfEaRRKyng0Avk3yHZLvlLe1hOQ/tfN9HZZX\ng1oL4FOSt5DsSvJyAMPbyP2M5E0ku5BsbCMXAP4IYEAbX18O4HSSY0geAeBfAPzazH5fxX7EImQ9\nQfJvSB68QHtEua5FFqyeJPsBWAXg38zsp1WuP0Yha3oRybNIHkayB4AZAHYBKMTv/niE/Jm/EMAZ\nAM4s/9kOYCJKQxM1lUuDMrNPUJqmawKwE8A4lEa/28q9HkALgGtQGhP9k2fzPwYwhaXfeZjs2N77\nKF1MvAelg3QogKs6sj+hhaxn2WsA9qH0scPTAD5igX+3LHA9/xHAlwDcWZ6i2ktyT0f2JwaBa9oT\npcGA3QBeR6m+/1Ce4i2kwO+hLWa24+AfAJ8C2F0eOKspmsX/iQ7JdQDmmNn80GupB6pntlTP7Kmm\n2SpqPaP8RV2SI0keXz49HQ9gMEr/UpcqqJ7ZUj2zp5pmq17qmdcUX6UGAViC0oXNLQDGmtm7YZdU\naKpntlTP7Kmm2aqLehbiIz4REel82jyDIqnulZKZtfuLg6pnemnqCaimaame2VI9s+eqabsf8eVx\nhrV//35n/L77knfdmDFjhjN39OjRzvjcuXOrX1hKZPobL8R2xnryySc748cff7wzvmrVKme8W7du\nzng1KqknkE9Nm5uTv7R/zz33OHMXLlzojGdZo0qErOfu3bud8QceeCAR8/1s9+7d2xm/7rrrErGm\npiZnbr9+2f0KVIzHp8/s2bMTsdtvv92Zu327+x6weRy3vppGOSQhIiKiBiUiIlFSgxIRkShFMWY+\nadIkZ3zevORz22bNct9dw/f5te+aSUND7jfmDc51HWXz5s3OXF/cd70w1PWVvFx88cWJmO/ayIoV\nK5zxK6+8MtM1FcG777onm5966qlEbOrUqc7cXbt2OeNTpkxJxHz/n/jeY+qF7+fS9b542mmVPVw4\n5M+8zqBERCRKalAiIhIlNSgREYmSGpSIiEQp1yEJ3y/tuYYhAGDy5OSd9H0XO30XUteuXeuMd8Yh\niauvvjp1ru8Xn3v27JnVcgrFdWHZN4Djq3NnHJIYNGiQM7569epEzFfPG264wRnv1atXItbY2FjB\n6urHbbfd5oy73hdfeOEFZ+4JJ5zgjIe8CYLOoEREJEpqUCIiEiU1KBERiZIalIiIREkNSkREopTr\nFF+lt8aYOHFi6lzfLU7qme8WJL6JHt/ti+QvfJOm55xzTiLmO543bNiQ6Zo6iwULFlSUv2XLlkSs\n3qdMFy9e7Iz7bvW2aNGiRKxPnz7O3JaWFmd86NChKVeXPZ1BiYhIlNSgREQkSmpQIiISJTUoERGJ\nkhqUiIhEKdcpvq1bt+b5cnVv586dzrhrugkABg4cmIj5JvuGDBlS/cIKzDcF5no4nk8lD3us9wc9\nVsI3iTZgwABn3HWvzjzuDxfS66+/XlH+zJkzEzHflK/PsGHDKsrPks6gREQkSmpQIiISJTUoERGJ\nkhqUiIhEiWbm/yJpbX29Ur5b83Tv3t0ZX79+fSI2ePBgZ67vQYZ33323M96vXz9nvBokYWZMkZdp\nPSvV3NyciA0fPtyZ63oYHOB/MGSW0taznBuspr4H7I0bN84Zz6N2LkWpp4/v9lOu4QnfA0p9D06s\nRsh6Vnp7M9fDYH23NHINUQHApk2bUq6uer6a6gxKRESipAYlIiJRUoMSEZEoqUGJiEiU1KBERCRK\nUTywcPTo0c74Pffck4j5bnvimzrLclqv6Hr06JE6tzM+ALItU6dOTcR8tz/yHYuubfjqfM011yRi\n3bt3xxFHHNHWMqPjmzpzPdRxz549ztw77rjDGXdNo7399tvO3Cyn+ELyvYdOnz7dGZ82bVoi5pua\nbmxsrH5hNaIzKBERiZIalIiIREkNSkREoqQGJSIiUcp1SMJn4cKFzrjr9h3r1q1z5i5ZsiTTNdWj\n/v37J2IjRoxw5q5Zs8YZ9130rvfnGjU1NSVivuduDR061BlfsGBBInbcccc5cxsaGpy59TIk4RqA\nqpTr/xNX3Toz13uob4hn4sSJtV5OxXQGJSIiUVKDEhGRKKlBiYhIlNSgREQkSmpQIiISpXYfWJjj\nWgot7QML81hLPajkgXC1Xks9UD2zpXpmz1XTNhuUiIhIKPqIT0REoqQGJSIiUVKDEhGRKKlBiYhI\nlNSgREQkSmpQIiISJTUoERGJkhqUiIhESQ1KRESiFKxBkZxH8q7yf59PcmPK70ud25montlSPbOn\nmmarM9QzijMoM/ulmZ1WTS7JN0he6Msn2Z/kAZJ7SO4t/2/yMZN1pJb1LOd0Jzmb5HskW0j+dweX\nHLUaH5/XtDou95DcVz5ez85i7bHK4Ri9guT/kfyA5KskGzu65pjlUM/rSb5ePkafJNm3o2tOI4oG\nlQMDcIyZHW1mPcxsWugFFdxPAfQEMAhAbwDfDbuc4jKzR1sdlz0ATAKw2cx+FXptRUXyBAD/AeA7\nZnYMgO8DeJTksWFXVkwkLwAwDcDXUfp5/wOAhXm8dm4NiuTZJF8u/4tmEYBurb42iuRbrf7+VZKv\nlHOXkFzU6lT281ySjwA4CcAT5c7+Pd/Lo86acah6khwE4DIAE81sl5UU/s008PHZ2ngAj2S6c4EE\nrOmJAFrM7FkAMLMnAewDMLBmO5uDgPW8FMBSM/udmX0K4G4AI0l+qYa7CyCnN22ShwNYDmA+Sh14\nKYCxh6RZq9zHAMwt5y4EMMaVa2bXAngTwGXlf4He51mCAfgDyTdJziXZp+N7FU7geg4HsBXAXeWP\n+H5D8vJMdiyQCI7Pg+voD+DvUQcNKnBNXwKwkeRlJA8jORrAfgD/m8W+hRDLMVp2sG+cUfmeVCav\ns4pzAHQ1s5lm9pmZLQPQ7Mk9F0AXM3ugnLscwPp2tt/Ws1neBzAMQH8AQwAcDWBBZcuPTsh6nghg\nMIAWAH0B3AJgfvnMqqhC1rO1awH8wsy2psyPWbCamtkBlD7iWwjgTwD+E8ANZvZxxXsRj5DH6NMA\nxpE8g2R3AHcAOADgCxXuQ8XyalAnANh2SMz3Q9jXkfuWKzENM9tnZq+Y2QEzew/AzQC+RvLIarcZ\ngWD1BPAxgD8DmGpmn5rZiwCeB/C1DmwztJD1bO1bAB7OaFuhBaspyYsA/CuAkWZ2OIALAPw7ya9U\nu80IhHwPXQXgTpTOyraU/+wF8Ha120wrrwb1DoB+h8ROqiD3i21su5onLhqKfU0qZD0PfkzS+l9c\nRX/qZfDjk+R5KL2xLEuTXwAha3omgBcOXhs1s5cA/A+Ai9r5vpgFPUbNbI6ZnWJmfVFqVF0BvNre\n93VUXm+NfZF6AAAGZUlEQVTSawF8SvIWkl3L1yyGt5H7GcmbSHZhaTzUlwsAfwQwwPdFksNJnsKS\nPgB+AuB5M9tb5b7EIFg9AbyI0mfW/1ze3nko/Qv1mYr3Ih4h63nQeADLzGxfRSuPV8iaNgM4n+SZ\nQGm4AMD5KPA1KIR9Dz2C5Onl/z4JwEMA7jezD6rakwrk0qDM7BMAlwNoArATwDh4/qXYKvd6lK5z\nXAPgCZQ+S3b5MYApJHeRnOz4+gCUPkPdg9IBur+8zcIKWc/yFE8jSpM9uwE8COBbZvb7juxTSIGP\nT5A8AsA3UD8f74U+Rl8E8CMA/0XyA5QGCqaZ2c87tFMBBT5Gu6E0pr8XwDoAq1G6DlVzNIv/0xmS\n6wDMMbP5oddSD1TPbKme2VNNs1XUekZ5HYbkSJLHl09Px6M0NfZ06HUVleqZLdUze6pptuqlnl1D\nL8BjEIAlKI0xbgEw1szeDbukQlM9s6V6Zk81zVZd1LMQH/GJiEjn0+YZFEl1r5TMrN1fxlQ900tT\nT0A1TUv1zJbqmT1XTdv9iC/LM6zmZvcvPt9zzz3O+I4dOxKxNWvWVPSaLS0tznjPnj0r2k5byLQ3\nCsi2npWaPXt2Inb77bc7c7dv3+6Md+vWzRnPUiX1BPKp6f79+xOxuXPnOnN9NW1qakrEpk+f3rGF\npRBjPW+99dZEbPhw9yT0zJkznfFLLrkkEfPVPksx1nPVqlXO+A033JCIrVy50pk7aFC4m8H4ahrl\nkISIiIgalIiIREkNSkREopTrmPmcOXOc8ccff9wZ79WrVyI2a9YsZ25DQ4MznuW1pqJ77rnnErHe\nvXs7c/O41hSjbdsOvcdmyRVXXJGIbdzofmq2r6YrVqxIxPK4BhUj18/2+vXuG24fd9xxzviMGTMS\nsZtvvtmZW+/vAwsWuB/QsHnz5kTsoYcecubGeCzqDEpERKKkBiUiIlFSgxIRkSipQYmISJTUoERE\nJEq5TvENHTrUGX/xxRed8ZEjRyZiEyZMcOZ21qkzF98kmmtactGiRbVeTqH47qBxzjnnJGKrV692\n5rrukgAAW7ZsqX5hdWbcuHGJ2L333uvMHTDA/Sw91yRgvU/r+VTy3uqafgSAKVOmOOMha6ozKBER\niZIalIiIREkNSkREoqQGJSIiUWrzgYUkLctbxbse9wAAN910U+ptDBw40BnftGlTVWvKAsnUz4MK\neev9iy66KBHL43EklUpbz3JuLjV1DZ74BiouvvhiZ9z1uI08LkzHWE/X40u6d+/uzJ08ebIzPm3a\ntEQsr8fBFKGeAHD11Ven3oZr6ATwP1YmS76a6gxKRESipAYlIiJRUoMSEZEoqUGJiEiU1KBERCRK\nuU7x+SZNfLeMcXFNogFAHpMyPrFN8S1evNgZv+qqq1JvY8SIEc74/fff74wPGzYs9bbbE+OUFJlq\nORUbPXq0M758+fLMXiPGeo4ZMyYR27FjhzPXN0U2aNCgTNeUVoz1zILvNnJ33323M96vX7/MXltT\nfCIiUihqUCIiEiU1KBERiZIalIiIREkNSkREopTrFF+lmpubE7Hhw4c7c99++21nPMtJE5/Ypvh6\n9+7tjLvuu+eb0PF5+OGHnfEs74UYckrKN2nqmiR77rnnnLkbNmxwxl33lGtsbHTm5jEh5ckNNsW3\ncOFCZ67vfnJZTjpWIsZ6ZsH1fgsAc+bMccazvEefpvhERKRQ1KBERCRKalAiIhIlNSgREYmSGpSI\niESpa54v5puQ8k09uZ5M6rtHXB7TekXhq+eoUaNSb+Pmm292xn1PgN29e3cidtRRR6Fr11wPsQ7z\nPZF10qRJidjmzZudub57yrm2Ue98P/MDBgxInes7nuUvfLXbunVr6m1s2bLFGZ83b54zPmPGjEQs\n6595nUGJiEiU1KBERCRKalAiIhIlNSgREYlSrlewfRfsXMMQgPvWPCtXrsx0TfXINzAybdq0ROzG\nG2905vqGIZqampzxnj17plxd/XAdnwBwySWX5LySePmGTly1Gzp0qDPXdwsk+YsVK1Y441k8pNT3\nM+/6//aww7I959EZlIiIREkNSkREoqQGJSIiUVKDEhGRKKlBiYhIlNp9YGGOaym0tA8szGMt9aCS\nB8LVei31QPXMluqZPVdN22xQIiIioegjPhERiZIalIiIREkNSkREoqQGJSIiUVKDEhGRKP0/v/ve\n0G9cRSgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1140eb668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pylab as pl\n",
    "\n",
    "num_rows = 4\n",
    "num_cols = 5\n",
    "\n",
    "fig, ax = plt.subplots(nrows=num_rows, ncols=num_cols, sharex=True, sharey=True)\n",
    "ax = ax.flatten()\n",
    "for index in range(num_rows*num_cols):\n",
    "    img = digits.images[index]\n",
    "    label = digits.target[index]\n",
    "    ax[index].imshow(img, cmap='Greys', interpolation='nearest')\n",
    "    ax[index].set_title('digit ' + str(label))\n",
    "\n",
    "ax[0].set_xticks([])\n",
    "ax[0].set_yticks([])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Date Preprocessing\n",
    "Hint: How you divide training and test data set? And apply other techinques we have learned if needed.\n",
    "You could take a look at the Iris data set case in the textbook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "scikit-learn version: 0.18\n",
      "1. Complete removing the mean and scaling to unit variance.\n",
      "2. Complete splitting with 1257(70%) training data and 540(30%) test data.\n"
     ]
    }
   ],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "if Version(sklearn_version) < '0.18':\n",
    "    from sklearn.cross_validation import train_test_split\n",
    "else:\n",
    "    from sklearn.model_selection import train_test_split\n",
    "print ('scikit-learn version: ' + str(Version(sklearn_version)))\n",
    "\n",
    "# 1. Standardize features by removing the mean and scaling to unit variance\n",
    "X_std = StandardScaler().fit_transform(X) # fit_transform(X) will fit to data, then transform it.\n",
    "print ('1. Complete removing the mean and scaling to unit variance.')\n",
    "\n",
    "# 2. splitting data into 70% training and 30% test data: \n",
    "split_ratio = 0.3\n",
    "X_train, X_test, y_train, y_test = train_test_split(X_std, y, test_size=split_ratio, random_state=0)\n",
    "print('2. Complete splitting with ' + str(y_train.shape[0]) + \\\n",
    "      '(' + str(int((1-split_ratio)*100)) +'%) training data and ' + \\\n",
    "      str(y_test.shape[0]) + '(' + str(int(split_ratio*100)) +'%) test data.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifier #1 Perceptron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Misclassified samples: 39 out of 540\n",
      "Accuracy: 0.928\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import Perceptron\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "# Training\n",
    "ppn = Perceptron(n_iter=800, eta0=0.1, random_state=0)\n",
    "ppn.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = ppn.predict(X_test)\n",
    "# Results\n",
    "print('Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifier #2 Logistic Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Misclassified samples: 25 out of 540\n",
      "Accuracy: 0.954\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# Training\n",
    "lr = LogisticRegression(C=1.0, random_state=0) # we observe that changing C from 0.0001 to 1000 has ignorable effect\n",
    "lr.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = lr.predict(X_test)\n",
    "# Results\n",
    "print('Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifier #3 SVM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. Using linear kernel:\n",
      "   Misclassified samples: 14 out of 540\n",
      "   Accuracy: 0.974\n",
      "2. Using rbf kernel:\n",
      "   Misclassified samples: 8 out of 540\n",
      "   Accuracy: 0.985\n"
     ]
    }
   ],
   "source": [
    "from sklearn.svm import SVC\n",
    "\n",
    "# 1. Using linear kernel\n",
    "\n",
    "# Training\n",
    "svm = SVC(kernel='linear', C=1.0, random_state=0)\n",
    "svm.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = svm.predict(X_test)\n",
    "# Results\n",
    "print('1. Using linear kernel:')\n",
    "print('   Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('   Accuracy: %.3f' % accuracy_score(y_test, y_pred))\n",
    "\n",
    "\n",
    "# 2. Using rbf kernel\n",
    "\n",
    "# Training\n",
    "svm = SVC(kernel='rbf', C=1.0, random_state=0)\n",
    "svm.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = svm.predict(X_test)\n",
    "# Results\n",
    "print('2. Using rbf kernel:')\n",
    "print('   Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('   Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifier #4 Decision Tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. Using entropy criterion:\n",
      "   Misclassified samples: 71 out of 540\n",
      "   Accuracy: 0.869\n",
      "2. Using Gini criterion:\n",
      "   Misclassified samples: 77 out of 540\n",
      "   Accuracy: 0.857\n"
     ]
    }
   ],
   "source": [
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "# 1. Using entropy criterion\n",
    "\n",
    "# Training\n",
    "tree = DecisionTreeClassifier(criterion='entropy', random_state=0)\n",
    "tree.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = tree.predict(X_test)\n",
    "# Results\n",
    "print('1. Using entropy criterion:')\n",
    "print('   Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('   Accuracy: %.3f' % accuracy_score(y_test, y_pred))\n",
    "\n",
    "# 2. Using Gini criterion\n",
    "\n",
    "# Training\n",
    "tree = DecisionTreeClassifier(criterion='gini', random_state=0)\n",
    "tree.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = tree.predict(X_test)\n",
    "# Results\n",
    "print('2. Using Gini criterion:')\n",
    "print('   Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('   Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifer #5 Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1. Using entropy criterion:\n",
      "   Misclassified samples: 32 out of 540\n",
      "   Accuracy: 0.941\n",
      "2. Using Gini criterion:\n",
      "   Misclassified samples: 27 out of 540\n",
      "   Accuracy: 0.950\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# 1. Using entropy criterion\n",
    "\n",
    "# Training\n",
    "forest = RandomForestClassifier(criterion='entropy', n_estimators=10, random_state=1, n_jobs=2)\n",
    "forest.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = forest.predict(X_test)\n",
    "# Results\n",
    "print('1. Using entropy criterion:')\n",
    "print('   Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('   Accuracy: %.3f' % accuracy_score(y_test, y_pred))\n",
    "\n",
    "# 2. Using Gini criterion\n",
    "\n",
    "# Training\n",
    "forest = RandomForestClassifier(criterion='gini', n_estimators=10, random_state=1, n_jobs=2)\n",
    "forest.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = forest.predict(X_test)\n",
    "# Results\n",
    "print('2. Using Gini criterion:')\n",
    "print('   Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('   Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifier #6 KNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Misclassified samples: 15 out of 540\n",
      "Accuracy: 0.972\n"
     ]
    }
   ],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "\n",
    "# Training\n",
    "knn = KNeighborsClassifier(n_neighbors=5, p=2, metric='minkowski')\n",
    "knn.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = knn.predict(X_test)\n",
    "# Results\n",
    "print('Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Classifier #7 Naive Bayes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Misclassified samples: 123 out of 540\n",
      "Accuracy: 0.772\n"
     ]
    }
   ],
   "source": [
    "from sklearn.naive_bayes import GaussianNB\n",
    "\n",
    "# Training\n",
    "gnb = GaussianNB()\n",
    "gnb.fit(X_train, y_train)\n",
    "# Testing\n",
    "y_pred = gnb.predict(X_test)\n",
    "# Results\n",
    "print('Misclassified samples: %d out of %d' % ((y_test != y_pred).sum(), y_test.shape[0]))\n",
    "print('Accuracy: %.3f' % accuracy_score(y_test, y_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence in this example, the accuracy of the predictions is ranked as (in descending order):\n",
    "* SVM (rbf kernel) - 0.985 <br>\n",
    "* SVM (linear kernel) - 0.974 <br>\n",
    "* KNN - 0.972 <br>\n",
    "* Logistic Regression - 0.954 <br>\n",
    "* Random forest (Gini criterion) - 0.950 <br>\n",
    "* Random forest (entropy criterion) - 0.941 <br>\n",
    "* Perceptron - 0.928 <br>\n",
    "* Decision Tree (entropy criterion) - 0.869 <br>\n",
    "* Decision Tree (Gini criterion) - 0.857 <br>\n",
    "* Naive Bayes - 0.772 <br>\n",
    "\n",
    "The best is SVM (using rbf kernel). Because SVM maximize margins to nearest samples (called support vectors), which is considered as an effective way of classifying spacially separated samples. Also, SVM is more robust against outliers and offers slack variables as a solution to not linearly-separable samples. Moreover, using rbf kernel, the dimension is mapped to infinity and therefore samples are highly likely to be separated by some hyperplanes.\n",
    "\n",
    "The worst is Naive Bayes. Because it assumes that all data are independent, which leads to high bias and low variance. When the input samples are not generally independent, the assumption fails and therefore the accuracy is low. Also, Naive Bayes cannot deal with outliers or noise, therefore unideal samples may not be correctly classified."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}