smd_unsupervised.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:26.174916Z",
     "start_time": "2018-11-27T15:05:26.151988Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:90% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:90% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unsupervised Machine Learning  and Clustering\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 id=\"tocheading\">Table of Contents</h2>\n",
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "% vectors and matrices\n",
    "\\def\\v#1{\\boldsymbol{#1}}\n",
    "\\def\\m#1{\\boldsymbol{#1}}\n",
    "\\def\\p#1{\\mathtt{#1}}\n",
    "% confusion matrix stuff\n",
    "\\def\\TP\\mathit{TP}\n",
    "\\def\\TN\\mathit{TN}\n",
    "\\def\\FP\\mathit{FP}\n",
    "\\def\\FN\\mathit{FN}\n",
    "% statistics \n",
    "\\DeclareMathOperator\\mse{mse}\n",
    "\\DeclareMathOperator\\E{E}\n",
    "\\DeclareMathOperator\\Var{Var}\n",
    "\\DeclareMathOperator\\Cov{Cov}\n",
    "\\DeclareMathOperator\\Bias{Bias}\n",
    "\\DeclareMathOperator*{\\argmax}{arg\\,max}\n",
    "\\DeclareMathOperator*{\\argmin}{arg\\,min}\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "The story so far:\n",
    "\n",
    "- Linear Discriminant Analysis (LDA) and Fisher's linear discriminant\n",
    "- Principal Component Analysis (PCA)\n",
    "- Feature Selection\n",
    "- Supervised Learning\n",
    "\n",
    "\n",
    "Notations:\n",
    "\n",
    "* Uppercase letters such as $X$ or $Y$ denote generic aspects of a variable (i.e. the actual random variable)\n",
    "* Observed values are written in lowercase. The ith observed value of $X$ is written as $x_i$\n",
    "* Matrices are written in bold uppercase letters as in $\\m{X}$\n",
    "* Observations map as *rows* in the matrix while the observed variables are the *columns*.\n",
    "\n",
    "So if I measure two observables $p = 2$ the size and weight of $N = 100$ people, I get a $N \\times p$ matrix $\\v{X}$.\n",
    "One observation in that matrix is denoted as $x_i = [ \\mathrm{size}, \\mathrm{weight} ]$ while all observations of the variable weight are denoted by $\\v{x}_{\\bullet 1}$, analogous to the numpy indexing `X[:, 1]`\n",
    "\n",
    "So far we have been occupied with \n",
    "predicting the values of one or more outputs or response variables $Y = (Y_1, \\ldots, Y_m)$ for a given set of input or predictor variables $X = (X_1, \\ldots , X_p)$. \n",
    "\n",
    "One possible definition of supervised machine learning\n",
    "\n",
    "> Given a $N \\times p$ matrix $\\v{X}$ and some associated output vector $\\v{Y} \\in \\mathbb{R}^N$,\n",
    " find a function $f(X) = \\hat{Y}$ that takes a vector $X \\in \\mathbb{R}^p$ and returns a prediction for $\\hat{Y}$\n",
    " where some \"loss function\" $L(Y, f(X))$ is minimized for all $X$.\n",
    " \n",
    "We've tried to split points into different classes by finding a decision function which seperates the points in an \"optimal\" way. This process is often called __classification__\n",
    "\n",
    "In __regression__ the dependent variable is not discrete but a continous value. The problem remains similar however. From *known* input data we try to find a function which accurately predicts $Y \\in \\mathbb{R}$.\n",
    "\n",
    "In __unsupervised__ machine learning problems the output vector $\\v{y}$ is *unknown*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:27.643365Z",
     "start_time": "2018-11-27T15:05:26.177785Z"
    }
   },
   "outputs": [],
   "source": [
    "from ml import plots\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib.colors import ListedColormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:27.643365Z",
     "start_time": "2018-11-27T15:05:26.177785Z"
    }
   },
   "outputs": [],
   "source": [
    "pd.options.display.max_rows = 10\n",
    "plots.set_plot_style()\n",
    "\n",
    "\n",
    "def create_discrete_cmap(k, under='k'):\n",
    "    cmap = ListedColormap([f'C{i}' for i in range(k)])\n",
    "    cmap.set_under('k')\n",
    "    return cmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification Example (Iris)\n",
    "\n",
    "The sterotypical machine learning data set. Given some measurements from flower petals is it possible to deduce the species for a new sample of measurements?\n",
    "\n",
    "Originally published in 1936 by Ronald Fisher in his paper about LDA. \n",
    "\n",
    "<table style=\"width:100%\">\n",
    "  <tr>\n",
    "    <td><img src=\"./ml/images/setosa.jpg\" alt=\"Drawing\" style=\"width: 200px;\"/></td>\n",
    "    <td><img src=\"./ml/images/versicolor.jpg\" alt=\"Drawing\" style=\"width: 200px;\"/></td> \n",
    "    <td><img src=\"./ml/images/virginica.jpg\" alt=\"Drawing\" style=\"width: 200px;\"/></td>\n",
    "  </tr>\n",
    "</table>\n",
    "\n",
    "<p style=\"color:gray\"> Iris flower data set. (2017, October 24). In Wikipedia, The Free Encyclopedia.  <p>\n",
    "\n",
    "Length and width of the petals and sepals making the dataset four dimensional.  Measurements where taken in Canada by Botanist Edgar Anderson:\n",
    "\n",
    "> all from the same pasture, and picked on the same day and measured at the same time by the same person with the same apparatus\n",
    "\n",
    "The label vector $\\v{y}$ is known because the botanist categorized the flowers into species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:28.673241Z",
     "start_time": "2018-11-27T15:05:27.646648Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b8d1998ab5a44d9cacb35b24e9da6e3d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from itertools import combinations\n",
    "\n",
    "\n",
    "iris = load_iris()\n",
    "X = iris.data  # we only take the first two features.\n",
    "y = iris.target\n",
    "\n",
    "cmap = create_discrete_cmap(len(iris.target_names))\n",
    "\n",
    "f, axs = plt.subplots(3, 2, constrained_layout=True)\n",
    "axs = axs.flatten()\n",
    "\n",
    "for i, (col_x, col_y) in enumerate(combinations(range(X.shape[1]), 2)):\n",
    "\n",
    "    plot = axs[i].scatter(\n",
    "        X[:, col_x], X[:, col_y],\n",
    "        c=y,\n",
    "        cmap=cmap,\n",
    "        label='species',\n",
    "        vmin=-0.5, vmax=len(iris.target_names) - 0.5,\n",
    "    )\n",
    "    axs[i].set_xlabel(iris.feature_names[col_x].replace('(cm)', '/ cm'))\n",
    "    axs[i].set_ylabel(iris.feature_names[col_y].replace('(cm)', '/ cm'))\n",
    "    \n",
    "bar = f.colorbar(plot)\n",
    "bar.set_ticks(np.arange(len(iris.target_names)))\n",
    "bar.set_ticklabels(iris.target_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The question this data set can help us answer is:\n",
    "\n",
    "> What species does this flower likely belong to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised Learning\n",
    "\n",
    "In unsupervised learning we have no given $Y$.\n",
    "\n",
    "These methods try to find the underlying (joint) probability density $P(X)$ so that we might learn some properties about it.\n",
    "\n",
    "One common question is whether $X$ is created by a mixture of two or more underlying random variables.\n",
    "\n",
    "In simpler terms:\n",
    "\n",
    " - Given some points, can we infer something about the underlying distributions or labels\n",
    "\n",
    "\n",
    "### Clustering\n",
    "\n",
    "Clustering algorithms try to find modes of $P(X)$ based on densities, neighborhood relations or any other measure of 'similarity'  between points.\n",
    "\n",
    "More generally speaking to quote wikipedia again:\n",
    "\n",
    "> Cluster analysis or clustering is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar (in some sense or another) to each other than to those in other groups (clusters).\n",
    "\n",
    "\n",
    "In the example below an unknown number of gaussian distributions are sampled to create a scatter plot.\n",
    "\n",
    "Can one infer $P(X)$ by looking at $X$, i.e. the blue dots?\n",
    "\n",
    "In this case we know that this distribution of blue dots, $X$, was created by sampling $k$ two-dimensional gaussians with known standard deviation.\n",
    "This is just what \n",
    "\n",
    "    X, y = make_blobs(n_samples=300, centers=k, center_box=(-2, 2), cluster_std=0.5)\n",
    " \n",
    "does.\n",
    "We even know in what region of space we have to look for to find the centroids centroids of these blobs. \n",
    "\n",
    "Questions: \n",
    "\n",
    "   - How many blobs do we have?\n",
    "   - Which points belong to which?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:28.984232Z",
     "start_time": "2018-11-27T15:05:28.675757Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a959ddfdf7604758aa7aea4c8176cba5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "rng = np.random.default_rng(1234)\n",
    "\n",
    "\n",
    "# choose a random number between 1 and 3\n",
    "k = rng.integers(1, 4)\n",
    "\n",
    "# create k blobs \n",
    "X, y = make_blobs(n_samples=300, centers=k, center_box=(-2, 2), cluster_std=0.5)\n",
    "\n",
    "# plot points without any color coding.\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_aspect('equal')\n",
    "ax.scatter(X[:, 0], X[:, 1])\n",
    "ax.set_axis_off()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "If you guessed k=2, you're wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:29.586188Z",
     "start_time": "2018-11-27T15:05:28.987673Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b0c3cc1ee304c75827458fa6122a9b2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmap = create_discrete_cmap(k)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_aspect('equal')\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap)\n",
    "ax.set_axis_off()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### K-Means Algorithm\n",
    "\n",
    "The k-Means algorithm tries to find a partition of the data into $k$ clusters $S = \\{S_1, \\ldots, S_k\\}$  which minimize the variance within those clusters. The number of clusters $k$ has to specified by the user.\n",
    "\n",
    "Formally speaking the algorithm solves\n",
    "$$\n",
    "{\\underset {S}{\\operatorname {arg\\,min} }}\\sum _{S_i \\in S}\\sum _{x \\in S_{i}}\\left\\|x -{\\overline{x}}_{S_i}\\right\\|^{2}.\n",
    "$$\n",
    "\n",
    "\n",
    "This problem is similar in nature to that of feature selection. A minimization over all possible subsets of the points $S$. Solving the general optimization problem is NP-Hard and therefore intractable. \n",
    "\n",
    "In simpler terms, it would take too long to find the optimal solution in the general case. \n",
    "Again there is a popular greedy heuristic which is usualy used to solve the problem.\n",
    "\n",
    "#### Loyds Algorithm\n",
    "\n",
    "Loyds Algorithm (sometimes also simply called *the* k-means algorithm) finds a local optimum using a greedy heuristic.\n",
    "\n",
    "It does so iterativly according to the following steps \n",
    "\n",
    "1. Pick some initial cluster means (or centroids) $\\{m_1, \\ldots, m_k \\}$ either randomly or according to some heuristic.\n",
    "\n",
    "2. Create a partition $S$  by assigning each point $x \\in X$ to the cluster $S_i$ where the distance to $m_i$ is the smallest.\n",
    "\n",
    "3. Update the cluster means by calculating the means within the assigned clusters. \n",
    "\n",
    "4. Repeat steps 2 and 3 until convergence.\n",
    "\n",
    "One popular convergence is the total change in cluster centroids per iteration. If the cluster centroids do not change from one iteration to the next, return the cluster centroids.\n",
    "\n",
    "\n",
    "The example below shows a k-mean clustering on $k=3$ random blobs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:30.581347Z",
     "start_time": "2018-11-27T15:05:29.591438Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ddfbcbea240444d99d83c99c38cb0d4f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "\n",
    "k = 3\n",
    "cmap = create_discrete_cmap(3)\n",
    "\n",
    "# create three random blobs\n",
    "X, y = make_blobs(n_samples=300, centers=k, center_box=(-2, 2), cluster_std=0.5, random_state=1234)\n",
    "\n",
    "# use KMeans to predict clusters for each sample in X\n",
    "kmeans = KMeans(n_clusters=3).fit(X)\n",
    "prediction = kmeans.predict(X)\n",
    "\n",
    "\n",
    "# match colors\n",
    "true_colors = [np.argmax(np.bincount(y[prediction == i])) for i in range(k)]\n",
    "\n",
    "# shift labels to get the right colors\n",
    "prediction = np.choose(prediction, true_colors)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_aspect(1)\n",
    "# plot shades with predicted clusters\n",
    "ax.scatter(X[:, 0], X[:, 1], c=prediction, cmap=cmap, alpha=0.3, s=250, lw=0, label='prediction')\n",
    "\n",
    "# plot points with true cluster associations\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap, label='truth')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set_axis_off()\n",
    "\n",
    "for center in kmeans.cluster_centers_:\n",
    "    ax.scatter(center[0], center[1], c=\"k\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below you can see an animation of k-Means as it converges. The black hexagons indicate the cluster centroids in the corresponding iteration. \n",
    "\n",
    "Create the video by running \n",
    "\n",
    "    python ani_kmeans.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video width=\"1280\" height=\"720\" style=\"max-width: 100%;\" controls>\n",
       "  <source src=\"kmeans_clustering.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%HTML\n",
    "<video width=\"1280\" height=\"720\" style=\"max-width: 100%;\" controls>\n",
    "  <source src=\"kmeans_clustering.mp4\" type=\"video/mp4\">\n",
    "</video>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the number of clusters in the center is not known, the algorithm does not produce meaningful results.\n",
    "This is a serious limitation as the number of clusters in the data is rarely known.  \n",
    "\n",
    "\n",
    "In the example below we sample two uniform 2D distributions and use k-means to cluster them into $k=4$ regions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:58.745695Z",
     "start_time": "2018-11-27T15:05:58.417540Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "041c92a98c4a44fb93a48669a9882bbe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "u1 = rng.uniform(-1, 0, size=(500, 2))\n",
    "u2 = rng.uniform(0, 1, size=(500, 2))\n",
    "X = np.append(u1, u2, axis=0)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(X[:, 0], X[:, 1], color='gray')\n",
    "ax.set_axis_off()\n",
    "ax.set_aspect(1)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:59.247113Z",
     "start_time": "2018-11-27T15:05:58.748135Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "85b6041c2aaf4f9db366a67449514ea2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use KMeans to predict clusters for each sample in X using the deliberatley wrong value of k=4\n",
    "k = 4\n",
    "prediction = KMeans(n_clusters=k).fit_predict(X)\n",
    "cmap = create_discrete_cmap(k)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(X[:, 0], X[:, 1], c=prediction, cmap=cmap)\n",
    "ax.set_axis_off()\n",
    "ax.set_aspect(1)\n",
    "ax.set_title(f'k = {k}')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The k-Means algorithm works well on convex clusters with similar standard deviations. But it fails on elongated or concave shapes.\n",
    "\n",
    "Like all cluster algorithms k-Means has a number of advantages and disadvantages\n",
    "\n",
    "Pros:\n",
    "\n",
    " - Relatively fast\n",
    " - Predictable results when inital cluster centroids are fixed.\n",
    " - works well on convex and 'blob' like shapes\n",
    "\n",
    "Cons:\n",
    "\n",
    " - Number of clusters has to be known beforehand\n",
    " - Gets worse and slower in high dimensions\n",
    " \n",
    "The follwing snippet transforms two random blobs to be elongated along the x axis.\n",
    "In this case k-Means fails as the neighborhood relation does not yield meaningful results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:00.102111Z",
     "start_time": "2018-11-27T15:05:59.663065Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "36b5a59ff0c4452e9856be0720ab9de2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "k = 2\n",
    "X, y = make_blobs(n_samples=1300, centers=k, random_state=2)\n",
    "\n",
    "transformation = [[20, 0], [0, 1]]\n",
    "X = np.dot(X, transformation)\n",
    "prediction = KMeans(n_clusters=2,).fit_predict(X)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.scatter(X[:, 0], X[:, 1], alpha=0.5, lw=0, c=prediction, cmap=create_discrete_cmap(k))\n",
    "ax.set_aspect('equal', 'datalim')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Scaling\n",
    "\n",
    "k-Means and many other algorithms, especially the ones using kernel functions, are sensitive to data scaling and centering.\n",
    "\n",
    "The scikit-learn toolkit includes a  preprocessing module with  methods to scale and transform data.\n",
    "\n",
    "##### Standard Scaling\n",
    "\n",
    "The simplest scaler is the `StandardScaler` scaling data to unit variance and zero mean. It makes it 'gaussian like'.\n",
    "\n",
    "##### MinMax Scaling\n",
    "\n",
    "Forces the values of each attribute in the data to be within the given feature range. When transforming  to values between 0 and 1 the applied transformation is\n",
    "\n",
    "$$\n",
    "\\mathbf{x}_j^{\\prime} = \\frac{(\\mathbf{x}_j - \\min(\\mathbf{x}_j))}{(\\max(\\mathbf{x}_j) - \\min(\\mathbf{x}_j))}\n",
    "$$\n",
    "\n",
    "\n",
    "##### Quantile Scaling\n",
    "\n",
    "Transform the given data to follow a normal, uniform or any other distribution with a CDF that can be inverted.\n",
    "Use inverse transform sampling to transform the data. I.e. apply the inverse CDF to each column in the data.\n",
    "\n",
    "https://en.wikipedia.org/wiki/Inverse_transform_sampling\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:01.555566Z",
     "start_time": "2018-11-27T15:06:00.104848Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "865f281ae62242048f5ffd5c815b09a1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn import preprocessing\n",
    "\n",
    "cmap = create_discrete_cmap(k)\n",
    "\n",
    "fig, axs = plt.subplots(4, 1, figsize=(9, 20), sharex=True)\n",
    "\n",
    "scalers=[\n",
    "        preprocessing.StandardScaler(),\n",
    "        preprocessing.MinMaxScaler(feature_range=(-1, 1)),\n",
    "        preprocessing.MaxAbsScaler(),\n",
    "        preprocessing.QuantileTransformer()\n",
    "]\n",
    "\n",
    "for scaler, ax in zip(scalers, axs):\n",
    "    X_prime = scaler.fit_transform(X)\n",
    "    prediction = KMeans(n_clusters=2, random_state=0).fit_predict(X_prime)\n",
    "    ax.set_title(scaler.__class__.__name__)\n",
    "    ax.scatter(X_prime[:, 0], X_prime[:, 1], c=prediction, cmap=cmap, alpha=0.5, lw=0)\n",
    "    ax.set_aspect(1)\n",
    "    \n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image Compression Using k-Means\n",
    "\n",
    "For a typical color, or RGB,  image the value for each color component must be stored. It is basically a list of color vectors \n",
    "\n",
    "$$\n",
    "c_i = \\begin{pmatrix} R, G, B \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "For a $420 \\times 640$ pixel image that results in  $420 \\cdot 640 \\cdot 3 \\; \\text{byte} = 793.8 \\;  \\text{kilobyte}$ of memory.\n",
    "This is how a typical Bitmap Image works.\n",
    "\n",
    "The GIF format allows the compression of color information.  \n",
    "\n",
    "In the example below we use k-Means to compress the color information of an image. We cluster the space of color components (either RGB, HSV or similar) into $k$ cluster centroids. \n",
    "Then only the association of each pixel to the cluster centroids and the centroids themselves have to be stored. \n",
    "\n",
    "When clustering the colors of an image into $k=50$ centroids for example  we only need to store $420 \\cdot 640 + k \\cdot 3  \\; \\text{byte} = 269.1 \\; \\text{kilobyte}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:03.176783Z",
     "start_time": "2018-11-27T15:06:01.557800Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e4fe55213691491083739be10d915c85",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import load_sample_image\n",
    "from skimage import color\n",
    "\n",
    "#load an example image as a 3D array\n",
    "image_rgb = load_sample_image(\"flower.jpg\")\n",
    "image_rgb = np.array(image_rgb, dtype=np.float64) / 255\n",
    "\n",
    "#store width length and number of colors\n",
    "width, length, channels = image_rgb.shape  \n",
    "\n",
    "# show image\n",
    "plt.figure()\n",
    "plt.imshow(image_rgb)\n",
    "\n",
    "#convert image to hsv for nicer plots\n",
    "image_hsv = color.rgb2hsv(image_rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:28.343954Z",
     "start_time": "2018-11-27T15:06:03.179701Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46bd448fc4e2415aa90c84fd1ad9c510",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot H, V and S values of each pixel into a 3D coordinate system.\n",
    "\n",
    "fig = plt.figure(figsize=(14, 12))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "\n",
    "ax.scatter3D(\n",
    "    image_hsv[:, :, 0],\n",
    "    image_hsv[:, :, 1],\n",
    "    image_hsv[:, :, 2],\n",
    "    c=image_rgb.reshape(-1, 3),\n",
    "    alpha=0.5,\n",
    "    lw=0,\n",
    ")\n",
    "\n",
    "ax.set_xlabel('hue')\n",
    "ax.set_ylabel('saturation')\n",
    "ax.set_zlabel('lightness value')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:53.066282Z",
     "start_time": "2018-11-27T15:06:28.346351Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b03cf29eba584f7b9ec2b569e59de6b3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(14, 12))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "\n",
    "ax.scatter3D(\n",
    "    image_rgb[:, :, 0],\n",
    "    image_rgb[:, :, 1],\n",
    "    image_rgb[:, :, 2],\n",
    "    c=image_rgb.reshape(-1, 3),\n",
    "    alpha=0.5,\n",
    "    lw=0,\n",
    ")\n",
    "ax.set_xlabel('red')\n",
    "ax.set_ylabel('green')\n",
    "ax.set_zlabel('blue')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:53.403058Z",
     "start_time": "2018-11-27T15:06:53.068733Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.utils import resample\n",
    "\n",
    "# perform k-means on a small random sample of pixels\n",
    "# we could use all pixels here but it would take much too long.\n",
    "# we reshape the image into a 2D array for resampling and k-means fitting\n",
    "flattened_image = image_hsv.reshape(-1, 3)\n",
    "\n",
    "sample = resample(flattened_image, n_samples=10000, replace=False)\n",
    "# sample = flattened_image\n",
    "\n",
    "\n",
    "# get the desired number of cluster centers\n",
    "kmeans = KMeans(n_clusters=50).fit(sample)\n",
    "centroids = kmeans.cluster_centers_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:53.840907Z",
     "start_time": "2018-11-27T15:06:53.405643Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cb80daac91a044e38604adf059f5989c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(10, 2))\n",
    "\n",
    "centroids_rgb = color.hsv2rgb(centroids.reshape(1, -1, 3)).reshape(-1, 3)\n",
    "#centroids_rgb = centroids.reshape(1, -1, 3).reshape(-1, 3)\n",
    "\n",
    "ax.scatter(\n",
    "    np.linspace(0, 1, len(centroids_rgb)),\n",
    "    np.ones(len(centroids_rgb)),\n",
    "    c=centroids_rgb,\n",
    "    s=10000.\n",
    ")\n",
    "ax.margins(0.1)\n",
    "ax.set_axis_off()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:54.792692Z",
     "start_time": "2018-11-27T15:06:53.843275Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2a2a1e2d44064bcc9de749fee02b97f7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# associate each pixel with the number of a cluster (i.e. the nearest centroid/color)\n",
    "labels = kmeans.predict(flattened_image)\n",
    "\n",
    "# get the actual value for each cluster centroid and reshape into a 3D image\n",
    "reconstructed  = centroids[labels].reshape(width, length, channels)\n",
    "\n",
    "# convert to RGB and plot again.\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))\n",
    "\n",
    "ax1.set_title('Compressed Using KMeans')\n",
    "ax1.imshow(color.hsv2rgb(reconstructed))\n",
    "\n",
    "ax2.set_title('Original')\n",
    "ax2.imshow(image_rgb)\n",
    "\n",
    "for ax in (ax1, ax2):\n",
    "    ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sorting colors using PCA\n",
    "\n",
    "When plotting the image centroids there are various ways to arrange the colors. There is no natural order \n",
    "to colors. Just as there is no natural order to most things with more than 1 dimension. \n",
    "\n",
    "We can try a few different things and try to get a nice output.\n",
    "\n",
    "1. Order by hue, Value, Saturation in HSV colorspace\n",
    "2. Order by red, green, blue in RGB colorspace\n",
    "3. Order by *length* of vector in e.g. eucledian norm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:55.375841Z",
     "start_time": "2018-11-27T15:06:54.795259Z"
    }
   },
   "outputs": [],
   "source": [
    "flattened_image = image_hsv.reshape(-1, 3)\n",
    "sample = resample(flattened_image, n_samples=10000, replace=False)\n",
    "\n",
    "# get the desired number of cluster centers\n",
    "kmeans = KMeans(n_clusters = 100).fit(sample)\n",
    "\n",
    "hsv = kmeans.cluster_centers_\n",
    "rgb = color.hsv2rgb(hsv.reshape(1, -1, 3)).reshape(-1, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:56.425540Z",
     "start_time": "2018-11-27T15:06:55.378062Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d485bd94d8a6455ca3541b8ba2c71b6f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig, axs = plt.subplots(6, 1, figsize=(12, 12))\n",
    "\n",
    "\n",
    "for i, label, ax in zip(range(3), 'rgb', axs[:3]):\n",
    "    \n",
    "    idx = np.argsort(rgb[:, i])\n",
    "    \n",
    "    ax.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "    ax.set_axis_off()\n",
    "    ax.set_title(f'Order by {label}')\n",
    "    ax.margins(0.1)\n",
    "    \n",
    "\n",
    "for i, label, ax in zip(range(3), 'hsv', axs[3:]):\n",
    "    \n",
    "    idx = np.argsort(hsv[:, i])\n",
    "    \n",
    "    ax.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "    ax.set_axis_off()\n",
    "    ax.set_title(f'Order by {label}')\n",
    "    ax.margins(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:57.328236Z",
     "start_time": "2018-11-27T15:06:56.937404Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ce09e20e17b74860afc444d7ffcd88bb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 6))\n",
    "\n",
    "\n",
    "norm = np.linalg.norm(rgb, axis=1)\n",
    "idx = np.argsort(norm)\n",
    "\n",
    "ax1.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax1.axis('off')\n",
    "ax1.set_title('Order by norm RGB')\n",
    "ax1.margins(0.1)\n",
    "\n",
    "norm = np.linalg.norm(centroids, axis=1)\n",
    "idx = np.argsort(norm)\n",
    "ax2.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax2.axis('off')\n",
    "ax2.set_title('Order by norm HSV')\n",
    "ax2.margins(0.1)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need some other way to reduce the information down to 1D so we can sort it. \n",
    "Previously we used the PCA to do just that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:57.333572Z",
     "start_time": "2018-11-27T15:06:57.330440Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:57.768544Z",
     "start_time": "2018-11-27T15:06:57.335763Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a573a7bc9f8d4a05aaa741de0c405fe2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2,) = plt.subplots(2, 1, figsize=(12, 6))\n",
    "\n",
    "c = PCA(n_components=1).fit_transform(rgb).ravel()\n",
    "idx = np.argsort(c)\n",
    "\n",
    "ax1.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax1.axis('off')\n",
    "ax1.set_title('Order by PCA transformed RGB.')\n",
    "ax1.margins(0.1)\n",
    "\n",
    "c = PCA(n_components=1).fit_transform(hsv).ravel()\n",
    "idx = np.argsort(c)\n",
    "\n",
    "ax2.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax2.axis('off')\n",
    "ax2.set_title('Order by PCA transformed HSV')\n",
    "ax2.margins(0.1)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density Based Clustering\n",
    "\n",
    "K-Means only takes distances between neighbours into acount. This works well for convex distributions that are approximately round in shape. Below you'll find another example where k-Means fails miserably. The moons are neither round nor convex yet easily separated into clusters by eye. A simple transformation will not help in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:58.165627Z",
     "start_time": "2018-11-27T15:06:57.771441Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4456aeca348a4f16b8dcb558c567f5a5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_moons\n",
    "\n",
    "X, y = make_moons(n_samples=400, noise=0.1, random_state=1234)\n",
    "\n",
    "kmeans = KMeans(n_clusters=2, random_state=0).fit(X)\n",
    "prediction = kmeans.predict(X)\n",
    "\n",
    "cmap = create_discrete_cmap(2)\n",
    "\n",
    "plt.figure()\n",
    "plt.axes(aspect=1)\n",
    "plt.scatter(X[:, 0], X[:, 1], c=prediction, s=200, lw=0, alpha=0.5, cmap=cmap)\n",
    "plt.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap)\n",
    "plt.axis('off')\n",
    "for center in kmeans.cluster_centers_:\n",
    "    plt.scatter(center[0], center[1], c=\"k\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DBSCAN\n",
    "\n",
    "Density-based spatial clustering of applications with noise (DBSCAN) is a clustering method based on the density of different regions within the parameter space. It clusters together points within regions of high density and labels points as outliers that do not lie within a dense region.\n",
    "\n",
    "1. Start by selecting *core points*. These are all points that have at least $m$ points in their neighbourhood region of radius $\\epsilon$.\n",
    "\n",
    "2. For each core point find connected points to build a cluster $C$. A point is connected to $C$ if its within the $\\epsilon$  neighbourhood of any point in $C$ and is also a core point.\n",
    "\n",
    "3. Assign each point which is not a core point to the nearest cluster with distance being at most $\\epsilon$ \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:00.159918Z",
     "start_time": "2018-11-27T15:06:58.167799Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "808fb826ea9b413ab05d8b5534bff346",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_moons, make_checkerboard, make_circles\n",
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "cmap = ListedColormap([f'C{i}' for i in range(2)])\n",
    "cmap.set_under('k')\n",
    "\n",
    "f, [top, center, bottom] = plt.subplots(3, 4, sharex=True, sharey=True)\n",
    "\n",
    "X_moon, y_moon = make_moons(noise=0.05, random_state=0)\n",
    "X_circle, y_circle = make_circles(noise=0.05, factor=0.4, random_state=0, n_samples=200)\n",
    "X_blobs, y_blobs = make_blobs(centers=2, center_box=(-0.5, 0.5), cluster_std=0.4, random_state=0)\n",
    "X_long, y_long = make_blobs(centers=2, center_box=(-2.1, 2.1), cluster_std=0.1, random_state=0)\n",
    "\n",
    "data = [(X_moon, y_moon), (X_long, y_long), (X_circle, y_circle), (X_blobs, y_blobs)]\n",
    "\n",
    "for ax, (X, y) in zip(top, data):\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=y, s=10, cmap=cmap, vmin=-0.5, vmax=1.5)\n",
    "    ax.axis('off')\n",
    "    ax.set_aspect(1)\n",
    "    \n",
    "for ax, (X, y) in zip(center, data):\n",
    "    prediction = KMeans(n_clusters=2, random_state=0).fit_predict(X)\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=prediction, s=10, cmap=cmap, vmin=-0.5, vmax=1.5)\n",
    "    ax.axis('off')\n",
    "    ax.set_aspect(1)\n",
    "    \n",
    "for ax, (X, y) in zip(bottom, data):\n",
    "    prediction = DBSCAN(eps=0.339).fit_predict(X)\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=prediction, s=10, cmap=cmap, vmin=-0.5, vmax=1.5)\n",
    "    ax.axis('off')\n",
    "    ax.set_aspect(1)\n",
    "    \n",
    "top[1].set_title('True Clusters')\n",
    "center[1].set_title('K - Mean Clusters')\n",
    "bottom[1].set_title('DBSCAN Clusters')\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluating Cluster Performance\n",
    "\n",
    "In the case of unsupervised learning there is no ground truth to which the cluster structure can be compared. \n",
    "\n",
    "Some heuristic has to be applied to measure how well a clustering performed. \n",
    "\n",
    "#### Silhuette Coefficent \n",
    "\n",
    "This evaluation ciriterion assumes a clustering is 'good' if the clusters are dense instead of sparse. \n",
    "\n",
    "Define $a$ as the mean distance between a single point $x_0$ and all other points in the cluster $S_p$ with $x_0 \\in S_p$.\n",
    "\n",
    "$$\n",
    "a(x_0) = \\frac{1}{|S_p|} \\sum_{x_i \\in S_p} \\left\\|x_i - x_0 \\right\\|\n",
    "$$\n",
    "\n",
    "and $b$ as the mean distance between $x_0$ and all the points in the *nearest cluster* $S_p^\\prime$ of which $x_0$ is not a member\n",
    "\n",
    "$$\n",
    "b(x_0) = \\frac{1}{|S_p^\\prime|} \\sum_{x_i \\in S_p^\\prime} \\left\\|x_i - x_0 \\right\\|\n",
    "$$\n",
    "\n",
    "The Silhuette Coefficent is then defined as \n",
    "\n",
    "$$\n",
    "s = \\frac{b - a}{\\text{max}(a, b)}\n",
    "$$\n",
    "\n",
    "\n",
    "The coefficent takes a value close to +1 for dense clustering and -1 for sparse clusters. \n",
    "\n",
    "Unfortunately it doesn't work very reliable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:00.582110Z",
     "start_time": "2018-11-27T15:07:00.162684Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c618a2003ec64502b304dc8a10461249",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.metrics import silhouette_score\n",
    "\n",
    "X, y = make_moons(n_samples=360, noise=0.1, random_state=172)\n",
    "\n",
    "km = KMeans(n_clusters=2)\n",
    "prediction_kmeans = km.fit_predict(X)\n",
    "score_kmeans = silhouette_score(X, km.labels_) \n",
    "\n",
    "db = DBSCAN(eps=0.14)\n",
    "prediction_db = db.fit_predict(X)\n",
    "score_db = silhouette_score(X, db.labels_) \n",
    "\n",
    "fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(10, 5))\n",
    "ax1.set_title('k-Means clustering score: {:0.3f}'.format(score_kmeans))\n",
    "ax1.scatter(X[:, 0], X[:, 1], c=prediction_kmeans, cmap=cmap, vmin=-0.5, vmax=1.5)\n",
    "ax1.axis('off')\n",
    "\n",
    "ax2.set_title('DBSCAN clustering: {:0.3f}'.format(score_db))\n",
    "ax2.scatter(X[:, 0], X[:, 1], c=prediction_db, cmap=cmap, vmin=-0.5, vmax=1.5)\n",
    "ax2.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomly Sampling Cluster Centers\n",
    "\n",
    "If we have the correct number of clusters and there are clear clusters, the kMeans algorithm should be stable. \n",
    "\n",
    "We can do many kMeans iterations with different initial centers and check if this is the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ffc34c80a5b34bfa9f8a4e4703950337",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/600 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from tqdm.auto import tqdm\n",
    "\n",
    "X, y = make_blobs(n_samples=300, centers=3, cluster_std=2, random_state=1234)\n",
    "\n",
    "\n",
    "centers = {}\n",
    "\n",
    "ks = range(1, 7)\n",
    "n_iter = 100\n",
    "bar = tqdm(total=len(ks) * n_iter)\n",
    "\n",
    "for k in ks:\n",
    "    centers[k] = []\n",
    "    for i in range(100):\n",
    "        bar.update(1)\n",
    "        kmeans = KMeans(k, random_state=i).fit(X)\n",
    "        centers[k].append(kmeans.cluster_centers_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7c771137e9ba41a9b8a0770a30696606",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(2, 3)\n",
    "\n",
    "\n",
    "\n",
    "for ax, (k, c),  in zip(axs.flat, centers.items()):\n",
    "    ax.scatter(X[:, 0], X[:, 1], s=5)\n",
    "    ax.set_aspect(1)\n",
    "    c = np.array(c)\n",
    "    ax.scatter(c[:, :, 0].ravel(), c[:, :, 1].ravel())\n",
    "    ax.set_title(f'k = {k}')"
   ]
  }
 ],
 "metadata": {
  "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}