K-Means without Scikit-Learn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "9183169c-9c32-4159-a460-f77b09480d79"
    }
   },
   "source": [
    "****\n",
    "# K-Means without Scikit-Learn\n",
    "****\n",
    "<p style=\"text-align:right\"><i>Jesus Perez Colino<br>First version: November 2016</i></p>\n",
    "\n",
    "## About this notebook: \n",
    "****\n",
    "Notebook prepared by **Jesus Perez Colino** Version 0.2, First Released: 01/10/2016, Alpha (work-in-progress)\n",
    "\n",
    "- This work is licensed under a [Creative Commons Attribution-ShareAlike 3.0 Unported License](http://creativecommons.org/licenses/by-sa/3.0/deed.en_US). This work is offered for free, with the hope that it will be useful.\n",
    "\n",
    "\n",
    "- **Summary**: This Jupyter notebook is the simplest Python implementation of the **K-Means algorithm**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------- Reproducibility conditions for this notebook -------------------\n",
      "Python version:       3.5.3 |Anaconda 4.4.0 (x86_64)| (default, Mar  6 2017, 12:15:08) \n",
      "[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]\n",
      "IPython version:      5.3.0\n",
      "-------------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "import IPython\n",
    "from sys import version \n",
    "\n",
    "print (' Reproducibility conditions for this notebook '.center(85,'-'))\n",
    "print ('Python version:       ' + version)\n",
    "print ('IPython version:      ' + IPython.__version__)\n",
    "print ('-'*85)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-Means algorithm\n",
    "\n",
    "Given a set of observations $(x_1,x_2,\\ldots, x_n)$, where each observation is a $d$-dimensional real vector, **$k$-means clustering** look for a partition of the $n$ observations into $k\\leq n$ sets or clusters $S=\\{S_1,S_2,\\ldots,S_k\\}$, such that the partition  $S$ *minimize* the within-cluster sum of squares between observations and the mean, or **centroid**, of the cluster (or alternatively, minimize the variance). \n",
    "\n",
    "More formally, the objective is to find a particular partition of clusters $S$ that is the solution of:\n",
    "\n",
    "\n",
    "$$ \\underset{S} {\\operatorname{arg\\,min}}  \\sum_{i=1}^{k} \\sum_{ x \\in S_i} \\left\\|  x - \\mu_i \\right\\|^2 = \\underset{S} {\\operatorname{arg\\,min}}  \\sum_{i=1}^{k} |S_i| \\operatorname{Var} S_i $$\n",
    "\n",
    "\n",
    "where $\\mu_i$ is the mean or centroid of points in the cluster $S_i$. \n",
    "\n",
    "The most common **k-means algorithm** uses an iterative refinement technique. Given an initial set of $k$ means $\\mu_1,…,\\mu_k$ (see below), the algorithm proceeds by alternating between two steps:\n",
    "\n",
    "- **Assignment step**: Assign each observation to the cluster whose mean has the least squared *Euclidean distance*, this is intuitively the \"nearest\" mean. \n",
    "\n",
    "\\begin{equation}\n",
    "S_i^{(t)} = \\big \\{ x_p : \\big \\| x_p - \\mu^{(t)}_i \\big \\|^2 \\le \\big \\| x_p - \\mu^{(t)}_j \\big \\|^2 \\ \\forall j, 1 \\le j \\le k \\big\\}\n",
    "\\end{equation} \n",
    "\n",
    "where each $x_p$ is assigned to exactly one $S^{(t)}$, even if it could be assigned to two or more of them.\n",
    "\n",
    "- **Update step**: Calculate the new means to be the **centroids** of the observations in the new clusters:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu^{(t+1)}_i = \\frac{1}{|S^{(t)}_i|} \\sum_{x_j \\in S^{(t)}_i} x_j \n",
    "\\end{equation}\n",
    "\n",
    "The algorithm will converge as soon as the assignments no longer change.  However, there is no guarantee that the optimum is found using this algorithm.\n",
    "\n",
    "The algorithm is often presented as assigning objects to the nearest cluster by distance.  Using a different distance function other than (squared) Euclidean distance may stop the algorithm from converging. Various modifications of k-means such as spherical k-means and k-medoids have been proposed to allow using other distance measures.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbpresent": {
     "id": "7b3839ec-1359-48c2-8945-09e0e7e39efc"
    }
   },
   "outputs": [],
   "source": [
    "from random import sample\n",
    "from math import fsum, sqrt\n",
    "from collections import defaultdict\n",
    "from functools import partial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbpresent": {
     "id": "6c932994-378c-42d9-8cd6-ebba30c8c621"
    }
   },
   "outputs": [],
   "source": [
    "def mean(data):\n",
    "    'Accurate arithmetic mean'\n",
    "    if isinstance(data,list)==False:\n",
    "        data = list(data)\n",
    "    return fsum(data) / len(data)\n",
    "\n",
    "def transpose(matrix):\n",
    "    'Swap rows with columns for a 2-D array'\n",
    "    return zip(*matrix)\n",
    "\n",
    "def distance(p, q, sqrt=sqrt, fsum=fsum, zip=zip):\n",
    "    'Multi-dimensional euclidean distance between points p and q'\n",
    "    return sqrt(fsum((x1 - x2) ** 2.0 for x1, x2 in zip(p, q)))\n",
    "\n",
    "def assign_data(centroids, data):\n",
    "    'Assign data the closest centroid'\n",
    "    d = defaultdict(list)\n",
    "    for point in data:\n",
    "        centroid = min(centroids, key=partial(distance, point))\n",
    "        d[centroid].append(point)\n",
    "    return dict(d)\n",
    "\n",
    "def compute_centroids(groups):\n",
    "    'Compute the centroid of each group'\n",
    "    return [tuple(map(mean, transpose(group))) for group in groups]\n",
    "\n",
    "def k_means(data, k=2, iterations=10):\n",
    "    'Return k-centroids for the data'\n",
    "    data = list(data)\n",
    "    centroids = sample(data, k)\n",
    "    for i in range(iterations):\n",
    "        labeled = assign_data(centroids, data)\n",
    "        centroids = compute_centroids(labeled.values())\n",
    "    return centroids\n",
    "\n",
    "def quality(labeled):\n",
    "    'Mean value of squared distances from data to its assigned centroid'\n",
    "    return mean(distance(c, p) ** 2 for c, pts in labeled.items() for p in pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "### Simple example with six 3-D points clustered into two groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{(17.666666666666668, 35.666666666666664, 25.0): [(10, 41, 23), (22, 30, 29), (21, 36, 23)], (14.333333333333334, 38.0, 7.0): [(11, 42, 5), (20, 32, 4), (12, 40, 12)]}\n"
     ]
    }
   ],
   "source": [
    "points=[(10, 41, 23),\n",
    "        (22, 30, 29),\n",
    "        (11, 42, 5),\n",
    "        (20, 32, 4),\n",
    "        (12, 40, 12),\n",
    "        (21, 36, 23)]\n",
    "\n",
    "centroids = k_means(points, k=2)\n",
    "print(assign_data(centroids, points))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Example with a richer dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbpresent": {
     "id": "637be784-e046-4c94-bfd7-9b6003f55307"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k     quality\n",
      "-     -------\n",
      "1     5,583.5\n",
      "2     1,337.8\n",
      "3       851.2\n",
      "4       666.6\n",
      "5       434.9\n",
      "6       239.5\n",
      "7       386.0\n"
     ]
    }
   ],
   "source": [
    "data = [ (10, 30),\n",
    "         (12, 50),\n",
    "         (14, 70),\n",
    "         (9, 150),\n",
    "         (20, 175),\n",
    "         (8, 200),\n",
    "         (14, 240),\n",
    "         (50, 35),\n",
    "         (40, 50),\n",
    "         (45, 60),\n",
    "         (55, 45),\n",
    "         (60, 130),\n",
    "         (60, 220),\n",
    "         (70, 150),\n",
    "         (60, 190),\n",
    "         (90, 160)]\n",
    "\n",
    "print('k     quality')\n",
    "print('-     -------')\n",
    "for k in range(1, 8):\n",
    "    centroids = k_means(data, k, iterations=20)\n",
    "    d = assign_data(centroids, data)\n",
    "    print('{0}    {1:8,.1f}'.format(k, quality(d)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py35]",
   "language": "python",
   "name": "conda-env-py35-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}