k-means.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# k-Means Algorithmus\n",
    "\n",
    "Die Aufgabe ist, gegebene Daten in möglichst ähnliche Cluster zu segmentieren.\n",
    "\n",
    "\"k\" bezieht sich dabei auf die Anzahl der Cluster, die man erhalten möchte."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Zunächst mal importieren wir ein paar Bibliotheken:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()  # for plot styling\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Grobe Beschreibung\n",
    "\n",
    "Der k-Means-Algorithmus sucht nach einer vorgegebenen Anzahl von Clustern innerhalb eines unmarkierten mehrdimensionalen Datensatzes.\n",
    "Dazu verwendet er eine einfache Vorstellung davon, wie das optimale Clustering aussieht:\n",
    "\n",
    "-  Das \"Clusterzentrum\" ist das arithmetische Mittel aller zum Cluster gehörenden Punkte.\n",
    "-  Jeder Punkt liegt näher an seinem eigenen Clusterzentrum als an anderen Clusterzentren.\n",
    "\n",
    "Diese beiden Annahmen bilden die Grundlage des k-Means-Algorithmus. Wir werden bald genau darauf\n",
    "eingehen, wie der Algorithmus diese Lösung erreicht, aber sehen wir uns zunächst einen einfachen\n",
    "Datensatz an und sehen uns das Ergebnis des k-Means an."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Unsere Daten\n",
    "Zunächst erzeugen wir uns einen Datensatz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "X, y_true = make_blobs(n_samples=300, centers=4,\n",
    "                       cluster_std=0.6, random_state=0)\n",
    "print(X.shape)\n",
    "fig, ax = plt.subplots(1, figsize=(10,10))\n",
    "ax.scatter(X[:, 0], X[:, 1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2,2, figsize=(10,10))\n",
    "axs[0,0].scatter(y_true, X[:,1]);\n",
    "axs[1,1].scatter(X[:,0], y_true);\n",
    "axs[0,1].scatter(X[:,0], X[:,1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "Man sieht schon mit bloßem Auge vier Cluster."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fertiger Algorithmus in Scikit-Learn\n",
    "\n",
    "Wir nutzen nun den k-Means Algorithmus von Scikit-learn. Das Muster ist dabei immer gleich:\n",
    "\n",
    "- zunächst wählt man ein Modell\n",
    "- und ruft dann die `.fit`-Methode, um das Modell zu trainieren\n",
    "- schließlich wendet man mit der `.predict`-Methode die Modell auf an\n",
    "\n",
    "Hier könnten wir auch die Kurzform `fit_predict` verwenden, wir arbeiten ja bei `fit` und `predict` mit dem gleichen Datensätzen.\n",
    "\n",
    "Hier ist die [Doku zu KMeans](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "kmeans = KMeans(n_clusters=4)\n",
    "kmeans.fit(X)\n",
    "y_kmeans = kmeans.predict(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Wir können das Ergebnis visualisieren - die großen roten Punkte zeigen die Zenter der Cluster an:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "centers = kmeans.cluster_centers_\n",
    "\n",
    "fig, ax = plt.subplots(1, figsize=(10,10))\n",
    "ax.scatter(X[:, 0], X[:, 1], c=y_kmeans, s=50, cmap='viridis')\n",
    "ax.scatter(centers[:, 0], centers[:, 1], c='red', s=100);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Der Algorithmus benötigte {} Iterationen\".format(kmeans.n_iter_))\n",
    "print(\"Inertia: {}\".format(kmeans.inertia_))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Implementierung\n",
    "\n",
    "Der Algorithmus ist einfach genug, dass wir ihn in wenigen Zeilen selbst schreiben können:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import pairwise_distances_argmin\n",
    "\n",
    "def find_clusters(X, n_clusters, rseed=2):\n",
    "    # 1. wähle zufällige Zentren:\n",
    "    rng = np.random.default_rng(rseed)\n",
    "    i = rng.permutation(X.shape[0])[:n_clusters]\n",
    "    centers = X[i]\n",
    "    n_iter = 0\n",
    "    while True:\n",
    "        # Färbe entsprechend dem nächtsgelegenen Zentrum ein\n",
    "        labels = pairwise_distances_argmin(X, centers)\n",
    "\n",
    "        # Ermittle neue Zentren als Mittelpunkt des Clusters\n",
    "        new_centers = np.array([X[labels == i].mean(0) for i in range(n_clusters)])\n",
    "\n",
    "        # Prüfe, ob fertig\n",
    "        if np.all(centers == new_centers):\n",
    "            break\n",
    "        centers = new_centers\n",
    "        n_iter = n_iter + 1\n",
    "\n",
    "    inertia = 0\n",
    "    for x in X:\n",
    "        closestDistance = float(\"inf\")\n",
    "        for y in centers:\n",
    "            d = (x-y).T.dot(x-y) # Vektor-Äquivalent von Euklidischer Norm zum Quadrat\n",
    "            if d < closestDistance:\n",
    "                closestDistance = d\n",
    "        inertia += closestDistance\n",
    "        \n",
    "    return centers, labels, n_iter, inertia\n",
    "\n",
    "centers, labels, n_iter, inertia = find_clusters(X, 4)\n",
    "\n",
    "print(\"Der Algorithmus benötigte {} Iterationen\".format(n_iter))\n",
    "print(\"Inertia: {}\".format(inertia))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Schauen wir uns das Ergebnis an:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, figsize=(10,10))\n",
    "ax.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis')\n",
    "ax.scatter(centers[:, 0], centers[:, 1], c='red', s=100);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Wir können uns das auch interaktiv anschauen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "from sklearn.metrics import pairwise_distances_argmin\n",
    "import ipywidgets as widgets\n",
    "\n",
    "def plot_kmeans_interactive(min_clusters=1, max_clusters=6):\n",
    "    def plot_points(ax, X, labels, n_clusters):\n",
    "        ax.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis', vmin=0, vmax=n_clusters - 1)\n",
    "\n",
    "    def plot_centers(ax, centers):\n",
    "        ax.scatter(centers[:, 0], centers[:, 1], marker='o',\n",
    "                    c=np.arange(centers.shape[0]),\n",
    "                    s=200, cmap='viridis')\n",
    "        ax.scatter(centers[:, 0], centers[:, 1], marker='o', c='red', s=100)\n",
    "\n",
    "    def _kmeans_step(frame=0, n_clusters=4, rseed=2):\n",
    "        labels = np.zeros(X.shape[0])\n",
    "        rng = np.random.default_rng(rseed)\n",
    "        i = rng.permutation(X.shape[0])[:n_clusters]\n",
    "        centers = X[i]\n",
    "\n",
    "        nsteps = frame // 3\n",
    "\n",
    "        for i in range(nsteps + 1):\n",
    "            old_centers = centers\n",
    "            if i < nsteps or frame % 3 > 0:\n",
    "                labels = pairwise_distances_argmin(X, centers)\n",
    "\n",
    "            if i < nsteps or frame % 3 > 1:\n",
    "                centers = np.array([X[labels == j].mean(0)\n",
    "                                    for j in range(n_clusters)])\n",
    "                nans = np.isnan(centers)\n",
    "                centers[nans] = old_centers[nans]\n",
    "\n",
    "        # plot the data and cluster centers\n",
    "        fig, ax = plt.subplots(1, figsize=(10,10))\n",
    "        ax.set_xlim(-6, 6)\n",
    "        ax.set_ylim(-2, 10)\n",
    "\n",
    "        plot_points(ax, X, labels, n_clusters)\n",
    "        plot_centers(ax, old_centers)\n",
    "\n",
    "        # plot new centers if third frame\n",
    "        if frame % 3 == 2:\n",
    "            for i in range(n_clusters):\n",
    "                ax.annotate('', centers[i], old_centers[i],\n",
    "                             arrowprops=dict(arrowstyle='->', linewidth=3, color='yellow'))\n",
    "            plot_centers(ax, centers)\n",
    "\n",
    "        if frame == 0:\n",
    "            ax.set_title(\"Initialisierung: Setze die Zentren zufällig\")\n",
    "        elif frame % 3 == 1:\n",
    "            ax.set_title(\"1. Färbe die Punkte entsprechend dem nächsten Zentrum ein\")\n",
    "        elif frame % 3 == 2:\n",
    "            ax.set_title(\"2. Setze die Zentren auf den Cluster-Durchschnitt\")\n",
    "        elif frame % 3 == 0 and frame > 0:\n",
    "            ax.set_title(\"Ergebnis: Neue Positionen der Zentren\")\n",
    "\n",
    "    return interact(_kmeans_step, frame=widgets.IntSlider(min=0, max=100, step=1, value=0),\n",
    "                    n_clusters=[1, 2, 3, 4, 5, 6, 15],\n",
    "                    rseed=[2,4])\n",
    "\n",
    "plot_kmeans_interactive();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Limitierungen des k-Means Algorithmus\n",
    "\n",
    "### Die Anzahl der Cluster muss im Vorfeld gewählt werden\n",
    "\n",
    "Wie man sieht ist die Wahl von k ganz entscheidend für eine saubere Trennung der Cluster.\n",
    "\n",
    "- Zu kleines k differenziert die Cluster ggf. nicht ausreichend.\n",
    "- Zu großes k trennt ggf. willkürlich Cluster auf.\n",
    "\n",
    "### Lösung ist nicht unbedingt die global beste Lösung\n",
    "\n",
    "Das Verfahren konvergiert zwar garantiert, allerdings nicht notwendigerweise zur global optimalen Lösung.\n",
    "\n",
    "Hier wirkt sich ganz entscheidend die anfängliche Anordnung der Zentren aus. Daher wird die Berechnung\n",
    "in der Praxis mehrere Male mit jeweils zufällig gewählten Anfangszentren wiederholt.\n",
    "\n",
    "Scikit-Learn macht das defaultmäßig 10 mal, gesteuert über den Parameter `n_init`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "centers, labels, n_iter, inertia = find_clusters(X, n_clusters=4, rseed=4)\n",
    "\n",
    "fig, ax = plt.subplots(1, figsize=(10,10))\n",
    "ax.scatter(X[:, 0], X[:, 1], c=labels, s=50, cmap='viridis')\n",
    "ax.scatter(centers[:, 0], centers[:, 1], c='red', s=100, alpha=0.5);\n",
    "print(\"Inertia: {}\".format(inertia))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wir sehen hier, dass für dieses suboptimale Ergebnis die Inertia wesentlich größer ist.\n",
    "\n",
    "Können wir das auch verwenden, um unterschiedliche Cluster-Anzahlen zu bewerten?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in np.append(np.arange(10),299):\n",
    "    k = n + 1\n",
    "    kmeans = KMeans(n_clusters=k)\n",
    "    kmeans.fit(X)\n",
    "    print('Anzahl Cluster: {} - Inertia: {}'.format(k, kmeans.inertia_))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nein. Das ist auch einleuchtend, wenn man sich überlegt, dass bei zunehmenden Clusteranzahlen die Wege\n",
    "der Punkte zu ihren Zentren natürlich abnimmt.\n",
    "\n",
    "Wenn jeder Punkt sein eigenes Cluster ist, dann ist Inertia = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Langsam, bei großen Datensets\n",
    "\n",
    "Da jede Iteration von k-Means auf jeden Punkt im Datensatz zugreifen muss, kann der Algorithmus\n",
    "relativ langsam sein, wenn die Anzahl der Stichproben wächst. Sie könnten sich fragen, ob diese\n",
    "Anforderung, alle Daten bei jeder Iteration zu verwenden, gelockert werden kann;\n",
    "beispielsweise könnten Sie einfach eine Teilmenge der Daten verwenden, um die Clusterzentren\n",
    "bei jedem Schritt zu aktualisieren.\n",
    "\n",
    "Dies ist die Idee hinter den batchbasierten k-Means-Algorithmen, von denen eine Form\n",
    "in `sklearn.cluster.MiniBatchKMeans` implementiert ist. Die Schnittstelle hierfür ist die\n",
    "gleiche wie für Standard-K-Means.\n",
    "\n",
    "### Limitiert auf lineare Cluster Grenzen\n",
    "\n",
    "Die Grenze zwischen zwei Clustern ist immer eine Gerade, daher hat der Algorithmus Probleme, wenn\n",
    "die Cluster kompliziertere Formen haben. Die Grenzen eines Clusters sind Linienabschnitte, sie bilden ein Voronoi:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.spatial import Voronoi, voronoi_plot_2d\n",
    "\n",
    "kmeans = KMeans(n_clusters=4)\n",
    "kmeans.fit(X)\n",
    "y_kmeans = kmeans.predict(X)\n",
    "\n",
    "centers = kmeans.cluster_centers_\n",
    "vor = Voronoi(centers)\n",
    "voronoi_plot_2d(vor, show_vertices=False, line_colors='orange', line_width=2, line_alpha=0.6, point_size=10);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Diese Form limitiert natürlich, was KMeans als Cluster differenzieren kann. Hier ein Beispiel, an dem KMeans scheitern muss:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_moons, make_circles\n",
    "\n",
    "X_1, y = make_moons(200, noise=.05, random_state=0)\n",
    "X_2, y = make_moons(200, noise=.05, random_state=0)\n",
    "X_strange = np.vstack([X_1, 2 + X_2])\n",
    "print(X_strange.shape)\n",
    "\n",
    "kmeans = KMeans(4, random_state=0)\n",
    "kmeans.fit(X_strange)\n",
    "y_means = kmeans.predict(X_strange)\n",
    "centers = kmeans.cluster_centers_\n",
    "print('Inertia: {}'.format(kmeans.inertia_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import DecisionPlot\n",
    "\n",
    "DecisionPlot.plot_decision_surface(X_strange, y_means, kmeans);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    },
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Alternative Algorithmen\n",
    "\n",
    "- Gaussian mixture models\n",
    "- Algorithmen, die Anzahl der Cluster selber wählen können:\n",
    "    - DBSCAN\n",
    "    - mean-shift\n",
    "    - affinity propagation\n",
    "    - ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}