07_Clasificacion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OvXvAf62ddu2"
   },
   "source": [
    "## Clasificación - Parte 1. Funciones Discriminantes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "pqoMFY44ddu6"
   },
   "source": [
    "Primero, corremos la celda de preparación."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "p6Qox9gcddu7"
   },
   "outputs": [],
   "source": [
    "# To support both python 2 and python 3\n",
    "from __future__ import division, print_function, unicode_literals\n",
    "\n",
    "# Common imports\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# to make this notebook's output stable across runs\n",
    "np.random.seed(42)\n",
    "\n",
    "# To plot pretty figures\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n",
    "\n",
    "# Where to save the figures\n",
    "PROJECT_ROOT_DIR = \".\"\n",
    "CHAPTER_ID = \"06_Regularizacion\"\n",
    "IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, \"plots\", CHAPTER_ID)\n",
    "\n",
    "def save_fig(fig_id, tight_layout=True, fig_extension=\"png\", resolution=300):\n",
    "    os.makedirs(IMAGES_PATH, exist_ok=True)\n",
    "    path = os.path.join(IMAGES_PATH, fig_id + \".\" + fig_extension)\n",
    "    print(\"Saving figure\", fig_id)\n",
    "    if tight_layout:\n",
    "        plt.tight_layout()\n",
    "    plt.savefig(path, format=fig_extension, dpi=resolution)\n",
    "\n",
    "# Ignore useless warnings (see SciPy issue #5998)\n",
    "# import warnings\n",
    "# warnings.filterwarnings(action=\"ignore\", message=\"^internal gelsd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discriminante lineal de Fischer\n",
    "\n",
    "Veamos cómo el criterio de Fischer nos permite encontrar la dirección óptima para realizar la proyección de las instancias (datos).\n",
    "\n",
    "Creemos un set aleatorio en dos dimensiones, compuesto de dos clases, cuya distribución es multinormal\n",
    "\n",
    "$$\n",
    "p(x | C_k) = \\mathcal{N}(\\mu_k, \\Sigma_k)\\;\\;.\n",
    "$$\n",
    "\n",
    "Por supuesto, si supieramos a priori que esta es la distribución, podríamos usar las técnicas que vamos a ver la próxima, para modelos discriminativos.\n",
    "\n",
    "Elijamos parámetros a ojo, y muestremos elementos de cada clase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from scipy.stats import multivariate_normal\n",
    "from numpy.random import multivariate_normal\n",
    "\n",
    "size1 = 250\n",
    "mu1 = [2, 1]\n",
    "cov1 = [[1, 0.95],[0.95, 1]]\n",
    "\n",
    "size2 = 200\n",
    "mu2 = [1, 3]\n",
    "cov2 = [[1, 0.2],[0.2, 1]]\n",
    "\n",
    "# Sample classes\n",
    "xc1 = multivariate_normal(mean=mu1, cov=cov1, size=size1).T\n",
    "\n",
    "xc2 = multivariate_normal(mean=mu2, cov=cov2, size=size2).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Veamos cómo se ven\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(*xc1, 'ob', mfc='None', label='C1')\n",
    "ax.plot(*xc2, 'or', mfc='None', label='C2')\n",
    "\n",
    "ax.set_xlabel('$x_1$')\n",
    "ax.set_ylabel('$x_2$')\n",
    "ax.legend(loc='lower right', fontsize=16)\n",
    "ax.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Claramente, si proyectamos, digamos sobre el eje de $x_2$, las clases no se separan muy bien."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Veamos cómo se ven\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "nbins = 20\n",
    "ax.hist(xc1[1], nbins, histtype='step', color='b')\n",
    "ax.hist(xc2[1], nbins, histtype='step', color='r')\n",
    "\n",
    "ax.set_xlabel('$y$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1S5Jyl0Cddu4"
   },
   "source": [
    "Veamos si podemos mejorar esto usando el criterio de Fischer. Para eso, tenemos que calcular la matriz de covarianza intra-clase\n",
    "\n",
    "$$\n",
    "\\mathbf{S}_I = \\underbrace{\\sum_{n\\,\\in\\,\\mathcal{C}_1} \\left(\\mathbf{x}_n - \\mathbf{m}_1\\right)\\left(\\mathbf{x}_n - \\mathbf{m}_1\\right)^\\mathrm{T}}_{\\text{Matriz de covarianza intra-clase 1}}\n",
    "+ \\underbrace{\\sum_{n\\,\\in\\,\\mathcal{C}_2} \\left(\\mathbf{x}_n - \\mathbf{m}_2\\right)\\left(\\mathbf{x}_n - \\mathbf{m}_2\\right)^\\mathrm{T}}_{\\text{Matriz de covarianza intra-clase 2}}\\;\\;,\n",
    "$$\n",
    "donde\n",
    "$$\n",
    "\\mathbf{m}_i = \\frac{1}{N_i}\\sum_{n\\,\\in\\,\\mathcal{C}_i} \\mathbf{x}_n\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = np.zeros([2,2])\n",
    "Si = np.zeros([2, 2, 2])\n",
    "for i, x in enumerate([xc1, xc2]):\n",
    "    # Calcula media empírica\n",
    "    m[i] = np.mean(x, axis=1)\n",
    "    xm = (x.T - m[i]).T\n",
    "    Si[i] = np.matmul(xm,xm.T)/(x.shape[1]-1)\n",
    "\n",
    "    print('Covariance Matrix for class {}\\n'.format(i+1), Si[i])\n",
    "    \n",
    "print(m[0], m[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sumemos ahora para obetener la matriz total\n",
    "S = np.sum(Si, axis=0)\n",
    "print('Covariance Matrix Total\\n', S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora resolvemos el problema lineal\n",
    "\n",
    "$$\n",
    "S_I \\mathbf{w} = \\mathbf{m}_2 - \\mathbf{m}_1\n",
    "$$\n",
    "\n",
    "para obtener $\\mathbf{w}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.linalg.solve(S, (m[1] - m[0])[:, np.newaxis])\n",
    "\n",
    "# Normalizamos, por amor al arte\n",
    "w /= np.linalg.norm(w)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Agregemos el vector al plot\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "ax.plot(*xc1, 'ob', mfc='None', label='C1')\n",
    "ax.plot(*xc2, 'or', mfc='None', label='C2')\n",
    "\n",
    "# Ploteo vector de pesos\n",
    "ax.quiver([0], [0], w[0], w[1], color='green', scale=10)\n",
    "\n",
    "# ploteo plano perpendicular\n",
    "xp = np.array([-1.5, 4])\n",
    "yp = -w[0] * xp / w[1]\n",
    "\n",
    "plt.plot(xp, yp, 'o:k')\n",
    "\n",
    "ax.set_xlabel('$x_1$')\n",
    "ax.set_ylabel('$x_2$')\n",
    "ax.legend(loc='lower right', fontsize=16)\n",
    "ax.set_aspect('equal')\n",
    "\n",
    "ax.spines['left'].set_position('zero')\n",
    "ax.spines['bottom'].set_position('zero')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a ver si me sale proyectar en la dirección perpendicular a $\\mathbf{w}$. Para eso, obtengo:\n",
    "\n",
    "$$\n",
    "\\mathbf{y} = \\mathbf{w}^\\mathrm{T} \\mathbf{x}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Encuentro proyección de cada clase\n",
    "yc1 = np.dot(w.T, xc1)\n",
    "yc2 = np.dot(w.T, xc2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hagamos nuevamente el histograma\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "nbins = 20\n",
    "ax.hist(yc1[0], nbins, histtype='step', color='b')\n",
    "ax.hist(yc2[0], nbins, histtype='step', color='r')\n",
    "\n",
    "ax.set_xlabel('$y$ óptimo')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xn = np.array([[4,],[-1,]])\n",
    "yn = np.dot(w.T, xn)\n",
    "\n",
    "if yn < 0:\n",
    "    print('xn pertenece a la Clase 1')\n",
    "else:\n",
    "    print('xn pertenece a la Clase 2')    \n",
    "    \n",
    "print(yn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1S5Jyl0Cddu4"
   },
   "source": [
    "## Perceptron\n",
    "\n",
    "Vimos cómo funciona el perceptron. Ahora vamos a usarlo para clasificar datos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Empecemos con los datos que generamos recién. Hagamos un set de entrenamiento y uno de test, juntando los elementos de ambas clases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.hstack([xc1, xc2])\n",
    "tc1 = np.ones(xc1.shape[1])\n",
    "tc2 = -np.ones(xc2.shape[1])\n",
    "t = np.hstack([tc1, tc2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "*Tarea*: Encontrar la superficie de decisión para los datos que generamos arriba usando el modelo del perceptron\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Datos de MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a usar un set de datos clásico de Machine Learning, los números de MNIST. Se trata de 70 000 imágenes pequeñas de dígitos escritos a mano. El \"target\" de cada uno de estos dígitos es el número que representan.\n",
    "\n",
    "Este set de datos es tan común, que en <tt>sklearn</tt> hay una función que permite bajarlos directamente. Dependiendo de la versión que <tt>sklearn</tt> que estén usando, la función relevante del paquete <tt>datasets</tt> cambia. Además, cambia la forma en la que devuelven los datos (antes estaban ordenados por valor del *target*, ahora viene así nomás. Para que el resultado sea idéntico con ambas versiones, usamos este código, que nos presta amablemente Géron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sort_by_target(mnist):\n",
    "    reorder_train = np.array(sorted([(target, i) for i, target in enumerate(mnist.target[:60000])]))[:, 1]\n",
    "    reorder_test = np.array(sorted([(target, i) for i, target in enumerate(mnist.target[60000:])]))[:, 1]\n",
    "    mnist.data[:60000] = mnist.data[reorder_train]\n",
    "    mnist.target[:60000] = mnist.target[reorder_train]\n",
    "    mnist.data[60000:] = mnist.data[reorder_test + 60000]\n",
    "    mnist.target[60000:] = mnist.target[reorder_test + 60000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    from sklearn.datasets import fetch_openml\n",
    "    mnist = fetch_openml('mnist_784', version=1, cache=True)\n",
    "    mnist.target = mnist.target.astype(np.int8) # fetch_openml() returns targets as strings\n",
    "    sort_by_target(mnist) # fetch_openml() returns an unsorted dataset\n",
    "except ImportError:\n",
    "    from sklearn.datasets import fetch_mldata\n",
    "    mnist = fetch_mldata('MNIST original')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El resultado es un tipo que todavía no habíamos visto. No vamos a entrar en detalles, pero digamos que tiene los datos en el atributo data y los valores de los labels en el atributo target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mnist.data.shape)\n",
    "print(mnist.target.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los datos tienen 784 *features*, que corresponden a cada uno de los píxeles de las imágenes de 28 x 28. El valor oscila entre 0 (blanco) y 255 (negro). Veamos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mnist.data[500])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Separemos ahora datos de labels. Seguimos usando nuestra notación, en la que los labels se llaman *t*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "X, t = mnist[\"data\"], mnist[\"target\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Agarremos un dígito cualquiera\n",
    "un_digito = X[36000]\n",
    "\n",
    "# Lo ponemos en forma de imagen y lo vemos.\n",
    "un_digito_image = un_digito.reshape(28, 28)\n",
    "plt.imshow(un_digito_image, cmap = mpl.cm.binary,\n",
    "           interpolation=\"None\")\n",
    "plt.axis(\"off\")\n",
    "\n",
    "#save_fig(\"some_digit_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Es un cinco (creo). Confirmemos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t[36000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lo que hicimos recién para plotear el número está piola. Vamos a convertirlo en una función para tener a mano."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_digit(data):\n",
    "    image = data.reshape(28, 28)\n",
    "    plt.imshow(image, cmap = mpl.cm.binary,\n",
    "               interpolation=\"nearest\")\n",
    "    plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora veamos varios números. Usamos otro código que nos vuelve a prestar nuestro amigo Aurélien."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# EXTRA\n",
    "def plot_digits(instances, images_per_row=10, **options):\n",
    "    size = 28\n",
    "    images_per_row = min(len(instances), images_per_row)\n",
    "    images = [instance.reshape(size,size) for instance in instances]\n",
    "    n_rows = (len(instances) - 1) // images_per_row + 1\n",
    "    row_images = []\n",
    "    n_empty = n_rows * images_per_row - len(instances)\n",
    "    images.append(np.zeros((size, size * n_empty)))\n",
    "    for row in range(n_rows):\n",
    "        rimages = images[row * images_per_row : (row + 1) * images_per_row]\n",
    "        row_images.append(np.concatenate(rimages, axis=1))\n",
    "    image = np.concatenate(row_images, axis=0)\n",
    "    plt.imshow(image, cmap = mpl.cm.binary, **options)\n",
    "    plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9,9))\n",
    "example_images = np.r_[X[:12000:600], X[13000:30600:600], X[30600:60000:590]]\n",
    "plot_digits(example_images, images_per_row=10)\n",
    "#save_fig(\"more_digits_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Antes de seguir mirando, tenemos que separar el conjunto en entrenamiento y testeo. Por suerte, los datos MNIST ya vienen separados, de forma de tener buena representación de cada clase. Las primeras 60000 instancias son de entrenamiento y las últimas 10000 de testeo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, t_train, t_test = X[:60000], X[60000:], t[:60000], t[60000:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a barajar el conjunto de entrenamiento. Esto es para asegurarse de que todos los dígitos estarán bien representandos en distintos *folds* de validación cruzada. (Ver clase anterior y [notebook](06_Regularización.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "shuffle_index = np.random.permutation(60000)\n",
    "X_train, t_train = X_train[shuffle_index], t_train[shuffle_index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clasificación binaria\n",
    "\n",
    "Vamos a empezar haciendo las cosas fáciles y dividir el problema. Vamos a intentar detectar \"cincos\". Para eso, generamos un nuevo label, que sea 1 cuando el número es cinco, y cero cuando no lo sea. Obviamente, hacemos lo mismo para el test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_train_5 = (t_train == 5)\n",
    "t_test_5 = (t_test == 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<tt>sklearn</tt> tiene una clase <tt>Perceptron</tt>, que es la que vamos a usar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Perceptron\n",
    "\n",
    "perce = Perceptron()\n",
    "\n",
    "# Hagamos un fit usando como features los valores de los píxeles directamente (linear regression)\n",
    "# En ese caso, la matriz de diseño es simplemente, Xtrain\n",
    "phi = X_train.copy()\n",
    "\n",
    "# Pero podríamos probar otras cosas, como esto que está comentado más abajo.\n",
    "#phi = np.log(1 + X_train.copy()**2)\n",
    "\n",
    "perce = perce.fit(phi, t_train_5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(perce.coef_[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Perceptron says', perce.predict([un_digito]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bueno, parece que ese cinco raro lo identifica correctamente. Pero claro, estaba dentro del conjunto de entrenamiento. Esto no quiere decir absolutamente nada.\n",
    "\n",
    "***\n",
    "**Pregunta**: ¿O sí?\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "#cross_val_score(sgd_clf, X_train, y_train_5, cv=3, scoring=\"accuracy\")\n",
    "scores = cross_val_score(perce, X_train, t_train_5, cv=5, scoring=\"accuracy\")\n",
    "print(scores)\n",
    "print(scores.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esto parece increíble. Más del 96% de \"accuracy\", con un modelo lineal muy simple. ¿Es posible?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "**Pregunta**: Pensemos un poco más en detenimiento, considerando la naturaleza del dataset. ¿Es realmente un valor tan alto? ¿Qué pasaría si hicieramos un clasificador que dijera que el número nunca es cinco?\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En efecto, en este tipo de datasets, donde los datos están muy desbalanceados, no es muy útil la \"accuracy\". Es mejor usar la matriz de confusión.\n",
    "\n",
    "Para crearla, neceistamos calcular predicciones en cada uno de nuestras instancias de entrenamiento, si queremos evitar usar el conjunto de test. En ese caso, podemos hacer CV y calcular las predicciones en cada uno de los folds. Eso lo hace una función de <tt>sklearn</tt> (gracias por tanto!), aunque la implementación no sería complicada."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_predict\n",
    "\n",
    "t_train_pred = cross_val_predict(perce, phi, t_train_5, cv=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "cm = confusion_matrix(t_train_5, t_train_pred)\n",
    "print(cm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con estos valores en mano, podemos calcular el exhaustividad (*recall*) y la precisión.\n",
    "\n",
    "Recordemos:\n",
    "\n",
    "$$\n",
    "\\mathrm{precision} = \\frac{TP}{TP + FP}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathrm{recall} = \\frac{TP}{TP + FN}\n",
    "$$\n",
    "\n",
    "donde TP son los *true positives* (es decir, la cantidad de casos relevantes recuperados correctamente), FP son los *falsos positivos* (es decir, la cantidad de casos no relevante recuperados incorrectamente), y FN son los *falsos negativos* (es decir, la cantidad de casos relevantes **no** recuperados).\n",
    "\n",
    "Calculemos estas cosas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tp = cm[1, 1]\n",
    "fp = cm[0, 1]\n",
    "fn = cm[1, 0]\n",
    "\n",
    "print('Precision: ', tp/(tp + fp))\n",
    "print('Recall: ',  tp/(tp + fn))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estos resultados hay que interpretarlos de la siguiente manera: la precisión nos dice en qué fracción de los casos en los que perceptrón dice que tiene un \"5\", realmente lo tiene. Además, el *recall* nos dice que fracción de todos los \"5\" encontró."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esto todavía no lo vimos, pero uno puede, en principio, ajustar estos valores mirando el nivel del umbral que usa para clasificar una instancia como verdadera o no. Como vimos, por defecto esto es $y(\\mathbf{x}) = 0$ para el perceptron, pero podemos cambiarlo. Para eso, necesitamos el perceptron nos diga el valor de $y(\\mathbf{x_k})$, para cada imagen $x_k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab_type": "text",
    "id": "3HNFr3DPddvC"
   },
   "outputs": [],
   "source": [
    "# Podríamos hacer así\n",
    "y = perce.decision_function(X_train)\n",
    "\n",
    "# Pero mejor sería hacerlo con CV, usando 5 folds\n",
    "y = cross_val_predict(perce, X_train, t_train_5, cv=3, method=\"decision_function\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Veamos cómo se ven\n",
    "A = plt.hist(y, 100)\n",
    "plt.xlabel('y(x)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora hagamos un código que vaya variando el umbral y nos calcule la precisión y el recall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "umb = np.linspace(y.min()-1, y.max()+1, 5000)\n",
    "\n",
    "recall = np.zeros_like(umb)\n",
    "precision = np.zeros_like(umb)\n",
    "fpr = np.zeros_like(umb)\n",
    "\n",
    "for i, u in enumerate(umb):\n",
    "    \n",
    "    # Calcula los índices con detecciones para este umbral\n",
    "    det = y > u\n",
    "    \n",
    "    # Compara esto con los verdaderos casos en esos índices\n",
    "\n",
    "    tp = np.sum(t_train_5[det] == True)\n",
    "    # Falsos positivos\n",
    "    fp = np.sum(t_train_5[det] == False)\n",
    "\n",
    "    # Verdaderos y falsos negativos\n",
    "    tn = np.sum(t_train_5[~det] == False)\n",
    "    fn = np.sum(t_train_5[~det] == True)\n",
    "    \n",
    "    recall[i] = tp/(tp + fn)\n",
    "    precision[i] = tp/(tp + fp)\n",
    "    fpr[i] = fp/(fp + tn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 5))\n",
    "plt.plot(umb, precision, label='Precisión', lw=2)\n",
    "plt.plot(umb, recall, label='Recall', lw=2)\n",
    "plt.xlabel(\"Umbral\", fontsize=16)\n",
    "plt.legend(loc=0, fontsize=16)\n",
    "plt.xlim(-umb.max(), umb.max())\n",
    "plt.axvline(0.0, ls=':', color='0.5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se puede plotear directo uno vs. el otro:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(7, 7))\n",
    "plt.plot(precision, recall)\n",
    "plt.axvline(0.1, ls=':', color='0.5')\n",
    "plt.xlabel('Precisión')\n",
    "plt.ylabel('Recall')\n",
    "plt.title('Curva PR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otra forma de ver esto es con la curva ROC (Receiver Operating Characteristic, o Característica Operativa del Receptor), que es muy similar, salvo que plotea la tasa de verdaderos positivos (es decir, el recall, en función de la tasa de falsos positivos). Cuanto más alejado de la recta unidad esté el sistema mejor, pero por supuesto esto depende del problema a resolver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(7, 7))\n",
    "plt.plot(fpr, recall)\n",
    "plt.plot([0, 1], [0, 1], color='0.5', ls=':')\n",
    "plt.xlabel('Tasa de falsos positivos (FPR)')\n",
    "plt.ylabel('Tasa de verdaderos positivos (TPR) / Recall')\n",
    "plt.title('Curva ROC')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otra característica interesante es el area bajo la curva. Podemos hacer una estimación veloz, sumando:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(recall[1:] * np.diff(fpr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hay código de <tt>sklearn</tt> para todo esto."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import precision_recall_curve, roc_curve, roc_auc_score\n",
    "precisions, recalls, thresholds = precision_recall_curve(t_train_5, y)\n",
    "fpr, tpr, thresholds = roc_curve(t_train_5, y)\n",
    "print(roc_auc_score(t_train_5, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extras\n",
    "\n",
    "¿Cómo se verá el cinco promedio?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cinco = X_train[t_train==5].mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_digit(cinco)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¿Qué distribución tienen los pesos en la imagen?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_digit(perce.coef_[0])"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "06_Regularización.ipynb",
   "provenance": []
  },
  "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.6.4"
  },
  "nav_menu": {
   "height": "279px",
   "width": "309px"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}