jupyter_russian/topic07_unsupervised/topic7_pca_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "</center>\n",
    "Авторы материала: исследователь СколТеха Сергей Королев и программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Тема 7. Обучение без учителя: PCA и кластеризация\n",
    "\n",
    "Привет всем и добро пожаловать в седьмую статью нашего открытого курса машинного обучения\n",
    "\n",
    "Данное занятие мы посвятим методам обучения без учителя (unsupervised learning), в частности методу главных компонент (PCA — principal component analysis) и кластеризации. Вы узнаете, зачем снижать размерность в данных, как это делать и какие есть способы группирования схожих наблюдений в данных.\n",
    "\n",
    "#### Список статей серии\n",
    "1. [Первичный анализ данных с Pandas](https://habrahabr.ru/company/ods/blog/322626/)\n",
    "2. [Визуальный анализ данных c Python](https://habrahabr.ru/company/ods/blog/323210/)\n",
    "3. [Классификация, деревья решений и метод ближайших соседей](https://habrahabr.ru/company/ods/blog/322534/)\n",
    "4. [Линейные модели классификации и регрессии](https://habrahabr.ru/company/ods/blog/323890/)\n",
    "5. [Композиции: бэггинг, случайный лес. Кривые валидации и обучения](https://habrahabr.ru/company/ods/blog/324402/)\n",
    "6. [Построение и отбор признаков](https://habrahabr.ru/company/ods/blog/325422/)\n",
    "7. [Обучение без учителя: PCA и кластеризация](https://habrahabr.ru/company/ods/blog/325654/)\n",
    "\n",
    "\n",
    "### План этой статьи\n",
    "\n",
    "1. Метод главных компонент (PCA)\n",
    " - Интуиция, теория и особенности применения\n",
    " - Примеры использования\n",
    "2. Кластеризация\n",
    " - K-means\n",
    " - Affinity Propagation\n",
    " - Спектральная кластеризация\n",
    " - Агломеративная кластеризация\n",
    " - Метрики качества кластеризации\n",
    "4. Полезные источники\n",
    "\n",
    "## 0. Введение\n",
    "\n",
    "Основное отличие методов обучения без учителя от привычных классификаций и регрессий машинного обучения в том, что разметки для данных в этом случае нет. От этого образуются сразу несколько особенностей — во-первых это возможность использования несопоставимо больших объёмов данных, поскольку их не нужно будет размечать руками для обучения, а во-вторых это неясность измерения качества методов, из-за отсутствия таких же прямолинейных и интуитивно понятных метрик, как в задачах обучения с учителем.\n",
    "\n",
    "Одной из самых очевидных задач, которые возникают в голове в отсутствие явной разметки, является задача снижения размерности данных. С одной стороны её можно рассматривать как помощь в визуализации данных, для этого часто используется метод t-SNE, который мы рассмотрели во второй статье курса. С другой стороны подобное снижение размерности может убрать лишние сильно скоррелированные признаки у наблюдений и подготовить данные для дальнейшей обработки в режиме обучения с учителем, например сделать входные данные более \"перевариваемыми\" для деревьев решений.\n",
    "\n",
    "## 1. Метод главных компонент (PCA)\n",
    "\n",
    "### Интуиция, теория и особенности применения\n",
    "\n",
    "Метод главных компонент (Principal Component Analysis) — один из самых интуитивно простых и часто используемых методов для снижения размерности данных и проекции их на ортогональное подпространство признаков.\n",
    "\n",
    "<img align=\"right\" src=\"https://habrastorage.org/getpro/habr/post_images/bb6/fe7/f06/bb6fe7f06e114bcc9c354a1cb025b966.png\" width=\"400\">\n",
    "\n",
    "В совсем общем виде это можно представить как предположение о том, что все наши наблюдения скорее всего выглядят как некий эллипсоид в подпространстве нашего исходного пространства и наш новый базис в этом пространстве совпадает с осями этого эллипсоида. Это предположение позволяет нам одновременно избавиться от сильно скоррелированных признаков, так как вектора базиса пространства, на которое мы проецируем, будут ортогональными.\n",
    "\n",
    "В общем случае размерность этого эллипсоида будет равна размерности исходного пространства, но наше предположение о том, что данные лежат в подпространстве меньшей размерности, позволяет нам отбросить \"лишнее\" подпространство в новой проекции, а именно то подпространство, вдоль осей которого эллипсоид будет наименее растянут. Мы будем это делать \"жадно\", выбирая по-очереди в качестве нового элемента базиса нашего нового подпространства последовательно ось эллипсоида из оставшихся, вдоль которой дисперсия будет максимальной.\n",
    "\n",
    "> \"To deal with hyper-planes in a 14 dimensional space, visualize a 3D space and say 'fourteen' very loudly. Everyone does it.\" - Geoffrey Hinton\n",
    "\n",
    "Рассмотрим как это делается математически:\n",
    "\n",
    "Чтобы снизить размерность наших данных из $n$ в $k, k \\leq n$, нам нужно выбрать топ-$k$ осей такого эллипсоида, отсортированные по убыванию по дисперсии вдоль осей.\n",
    "\n",
    "Начнём с того, что посчитаем дисперсии и ковариации исходных признаков. Это делается просто с помощью матрицы ковариации. По определению ковариации, для двух признаков $X_i$ и $X_j$ их ковариация будет \n",
    "$$cov(X_i, X_j) = E[(X_i - \\mu_i) (X_j - \\mu_j)] = E[X_i X_j] - \\mu_i \\mu_j$$\n",
    "где $\\mu_i$ — матожидание $i$-ого признака.\n",
    "При этом отметим, что ковариация симметрична и ковариация вектора с самим собой будет равна его дисперсии.\n",
    "\n",
    "Таким образом матрица ковариации представляет собой симметричную матрицу, где на диагонали лежат дисперсии соответствующих признаков, а вне диагонали — ковариации соответствующих пар признаков. В матричном виде, где $\\mathbf{X}$ это матрица наблюдений, наша матрица ковариации будет выглядеть как\n",
    "$$\\Sigma = E[(\\mathbf{X} - E[\\mathbf{X}]) (\\mathbf{X} - E[\\mathbf{X}])^{T}]$$\n",
    "\n",
    "Чтобы освежить память — у матриц как у линейных операторов есть такое интересное свойство как собственные значения и собственные вектора (eigenvalues и eigenvectors). Эти штуки замечательны тем, что когда мы нашей матрицей действуем на соответствующее линейное пространство, собственные вектора остаются на месте и лишь умножаются на соответствующие им собственные значения. То есть определяют подпространство, которое при действии этой матрицей как линейным оператором, остаётся на месте или \"переходит в себя\". Формально собственный вектор $w_i$ с собственным значением $\\lambda_i$ для матрицы $M$ определяется просто как $M w_i = \\lambda_i w_i$.\n",
    "\n",
    "Матрицу ковариации для нашей выборки $\\mathbf{X}$ можно представить в виде произведения $\\mathbf{X}^{T} \\mathbf{X}$. Из [отношения Релея](https://ru.wikipedia.org/wiki/%D0%9E%D1%82%D0%BD%D0%BE%D1%88%D0%B5%D0%BD%D0%B8%D0%B5_%D0%A0%D1%8D%D0%BB%D0%B5%D1%8F) вытекает, что максимальная вариация нашего набора данных будет достигаться вдоль собственного вектора этой матрицы, соответствующего максимальному собственному значению. Таким образом главные компоненты, на которые мы бы хотели спроецировать наши данные, являются просто собственными векторами соответствующих топ-$k$ штук собственных значений этой матрицы.\n",
    "\n",
    "Дальнейшие шаги просты до безобразия — надо просто умножить нашу матрицу данных на эти компоненты и мы получим проекцию наших данных в ортогональном базисе этих компонент. Теперь если мы транспонируем нашу матрицу данных и матрицу векторов главных компонент, мы восстановим исходную выборку в том пространстве, из которого мы делали проекцию на компоненты. Если количество компонент было меньше размерности исходного пространства, мы потеряем часть информации при таком преобразовании."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Примеры использования\n",
    "### Набор данных по цветкам ириса\n",
    "Начнём с того, что загрузим все необходимые модули и покрутим привычный датасет с ирисами по примеру из документации пакета scikit-learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "sns.set(style=\"white\")\n",
    "%matplotlib inline\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from sklearn import datasets, decomposition\n",
    "\n",
    "# Загрузим наши ириски\n",
    "iris = datasets.load_iris()\n",
    "X = iris.data\n",
    "y = iris.target\n",
    "\n",
    "# Заведём красивую трёхмерную картинку\n",
    "fig = plt.figure(1, figsize=(6, 5))\n",
    "plt.clf()\n",
    "ax = Axes3D(fig, rect=[0, 0, 0.95, 1], elev=48, azim=134)\n",
    "\n",
    "plt.cla()\n",
    "\n",
    "for name, label in [(\"Setosa\", 0), (\"Versicolour\", 1), (\"Virginica\", 2)]:\n",
    "    ax.text3D(\n",
    "        X[y == label, 0].mean(),\n",
    "        X[y == label, 1].mean() + 1.5,\n",
    "        X[y == label, 2].mean(),\n",
    "        name,\n",
    "        horizontalalignment=\"center\",\n",
    "        bbox=dict(alpha=0.5, edgecolor=\"w\", facecolor=\"w\"),\n",
    "    )\n",
    "# Поменяем порядок цветов меток, чтобы они соответствовали правильному\n",
    "y_clr = np.choose(y, [1, 2, 0]).astype(np.float)\n",
    "ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y_clr, cmap=plt.cm.nipy_spectral)\n",
    "\n",
    "ax.w_xaxis.set_ticklabels([])\n",
    "ax.w_yaxis.set_ticklabels([])\n",
    "ax.w_zaxis.set_ticklabels([]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь посмотрим, насколько PCA улучшит результаты для модели, которая в данном случае плохо справится с классификацией из-за того, что у неё не хватит сложности для описания данных:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score, roc_auc_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "# Выделим из наших данных валидационную выборку\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.3, stratify=y, random_state=42\n",
    ")\n",
    "\n",
    "# Для примера возьмём неглубокое дерево решений\n",
    "clf = DecisionTreeClassifier(max_depth=2, random_state=42)\n",
    "clf.fit(X_train, y_train)\n",
    "preds = clf.predict_proba(X_test)\n",
    "print(\"Accuracy: {:.5f}\".format(accuracy_score(y_test, preds.argmax(axis=1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Теперь попробуем сделать тоже самое, но с данными, для которых мы снизили размерность до 2D:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Прогоним встроенный в sklearn PCA\n",
    "pca = decomposition.PCA(n_components=2)\n",
    "X_centered = X - X.mean(axis=0)\n",
    "pca.fit(X_centered)\n",
    "X_pca = pca.transform(X_centered)\n",
    "\n",
    "# И нарисуем получившиеся точки в нашем новом пространстве\n",
    "plt.plot(X_pca[y == 0, 0], X_pca[y == 0, 1], \"bo\", label=\"Setosa\")\n",
    "plt.plot(X_pca[y == 1, 0], X_pca[y == 1, 1], \"go\", label=\"Versicolour\")\n",
    "plt.plot(X_pca[y == 2, 0], X_pca[y == 2, 1], \"ro\", label=\"Virginica\")\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Повторим то же самое разбиение на валидацию и тренировочную выборку.\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X_pca, y, test_size=0.3, stratify=y, random_state=42\n",
    ")\n",
    "\n",
    "clf = DecisionTreeClassifier(max_depth=2, random_state=42)\n",
    "clf.fit(X_train, y_train)\n",
    "preds = clf.predict_proba(X_test)\n",
    "print(\"Accuracy: {:.5f}\".format(accuracy_score(y_test, preds.argmax(axis=1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видно, что качество возросло незначительно, но для более сложных данных более высокой размерности, где данные не разбиваются тривиально вдоль одного признака, применение PCA может достаточно сильно улучшить качество работы деревьев решений и ансамблей на их основе."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Посмотрим на 2 главные компоненты в последнем PCA-представлении данных и на тот процент исходной дисперсии в даных, который они \"объясняют\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, component in enumerate(pca.components_):\n",
    "    print(\n",
    "        \"{} component: {}% of initial variance\".format(\n",
    "            i + 1, round(100 * pca.explained_variance_ratio_[i], 2)\n",
    "        )\n",
    "    )\n",
    "    print(\n",
    "        \" + \".join(\n",
    "            \"%.3f x %s\" % (value, name)\n",
    "            for value, name in zip(component, iris.feature_names)\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Набор данных по рукописным цифрам\n",
    "Теперь возьмем набор данных по рукописным цифрам. Мы с ним уже работали в 3 [статье](https://habrahabr.ru/company/ods/blog/322534/#derevya-resheniy-i-metod-blizhayshih-sosedey-v-zadache-raspoznavaniya-rukopisnyh-cifr-mnist) про деревья решений и метод ближайших соседей."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "digits = datasets.load_digits()\n",
    "X = digits.data\n",
    "y = digits.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Вспомним, как выглядят эти цифры – посмотрим на первые десять. Картинки здесь представляются матрицей 8 x 8 (интенсивности белого цвета для каждого пикселя). Далее эта матрица \"разворачивается\" в вектор длины 64, получается признаковое описание объекта."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# f, axes = plt.subplots(5, 2, sharey=True, figsize=(16,6))\n",
    "plt.figure(figsize=(16, 6))\n",
    "for i in range(10):\n",
    "    plt.subplot(2, 5, i + 1)\n",
    "    plt.imshow(X[i, :].reshape([8, 8]));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Получается, размерность признакового пространства здесь – 64. Но давайте снизим размерность всего до 2 и увидим, что даже на глаз рукописные цифры неплохо разделяются на кластеры."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = decomposition.PCA(n_components=2)\n",
    "X_reduced = pca.fit_transform(X)\n",
    "\n",
    "print(\"Projecting %d-dimensional data to 2D\" % X.shape[1])\n",
    "\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.scatter(\n",
    "    X_reduced[:, 0],\n",
    "    X_reduced[:, 1],\n",
    "    c=y,\n",
    "    edgecolor=\"none\",\n",
    "    alpha=0.7,\n",
    "    s=40,\n",
    "    cmap=plt.cm.get_cmap(\"nipy_spectral\", 10),\n",
    ")\n",
    "plt.colorbar()\n",
    "plt.title(\"MNIST. PCA projection\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ну, правда, с t-SNE картинка получается еще лучше, поскольку у PCA ограничение – он находит только линейные комбинации исходных признаков. Зато даже на этом относительно небольшом наборе данных можно заметить, насколько t-SNE дольше работает."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "from sklearn.manifold import TSNE\n",
    "\n",
    "tsne = TSNE(random_state=17)\n",
    "\n",
    "X_tsne = tsne.fit_transform(X)\n",
    "\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.scatter(\n",
    "    X_tsne[:, 0],\n",
    "    X_tsne[:, 1],\n",
    "    c=y,\n",
    "    edgecolor=\"none\",\n",
    "    alpha=0.7,\n",
    "    s=40,\n",
    "    cmap=plt.cm.get_cmap(\"nipy_spectral\", 10),\n",
    ")\n",
    "plt.colorbar()\n",
    "plt.title(\"MNIST. t-SNE projection\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На практике, как правило, выбирают столько главных компонент, чтобы оставить 90% дисперсии исходных данных. В данном случае для этого достаточно выделить 21 главную компоненту, то есть снизить размерность с 64 признаков до 21."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = decomposition.PCA().fit(X)\n",
    "\n",
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot(np.cumsum(pca.explained_variance_ratio_), color=\"k\", lw=2)\n",
    "plt.xlabel(\"Number of components\")\n",
    "plt.ylabel(\"Total explained variance\")\n",
    "plt.xlim(0, 63)\n",
    "plt.yticks(np.arange(0, 1.1, 0.1))\n",
    "plt.axvline(21, c=\"b\")\n",
    "plt.axhline(0.9, c=\"r\")\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 2. Кластеризация\n",
    "\n",
    "Интуитивная постановка задачи кластеризации довольно проста и представляет из себя наше желание сказать: \"Вот тут у меня насыпаны точки. Я вижу, что они сваливаются в какие-то кучки вместе. Было бы круто иметь возможность эти точки относить к кучкам и в случае появления новой точки на плоскости говорить, в какую кучку она падает.\" Из такой постановки видно, что пространства для фантазии получается много, и от этого возникает соответствующее множество алгоритмов решения этой задачи. Перечисленные алгоритмы ни в коем случае не описывают данное множество полностью, но являются примерами самых популярных методов решения задачи кластеризации.\n",
    "\n",
    "<figure><img align=\"center\" src=\"https://habrastorage.org/getpro/habr/post_images/8b9/ae5/586/8b9ae55861f22a2809e8b3a00ef815ad.png\"><figcaption>Примеры работы алгоритмов кластеризации из документации пакета scikit-learn</figcaption></figure>\n",
    "\n",
    "### K-means\n",
    "\n",
    "Алгоритм К-средних, наверное, самый популярный и простой алгоритм кластеризации и очень легко представляется в виде простого псевдокода:\n",
    "1. Выбрать количество кластеров $inline$k$inline$, которое нам кажется оптимальным для наших данных.\n",
    "2. Высыпать случайным образом в пространство наших данных $inline$k$inline$ точек (центроидов).\n",
    "3. Для каждой точки нашего набора данных посчитать, к какому центроиду она ближе.\n",
    "4. Переместить каждый центроид в центр выборки, которую мы отнесли к этому центроиду.\n",
    "5. Повторять последние два шага фиксированное число раз, либо до тех пор пока центроиды не \"сойдутся\" (обычно это значит, что их смещение относительно предыдущего положения не превышает какого-то заранее заданного небольшого значения).\n",
    "\n",
    "В случае обычной евклидовой метрики для точек лежащих на плоскости, этот алгоритм очень просто расписывается аналитически и рисуется. Давайте посмотрим соответствующий пример:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Начнём с того, что насыпем на плоскость три кластера точек\n",
    "X = np.zeros((150, 2))\n",
    "\n",
    "np.random.seed(seed=42)\n",
    "X[:50, 0] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "X[:50, 1] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "\n",
    "X[50:100, 0] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "X[50:100, 1] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "\n",
    "X[100:150, 0] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "X[100:150, 1] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "\n",
    "plt.figure(figsize=(5, 5))\n",
    "plt.plot(X[:, 0], X[:, 1], \"bo\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# В scipy есть замечательная функция, которая считает расстояния\n",
    "# между парами точек из двух массивов, подающихся ей на вход\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "# Прибьём рандомность и насыпем три случайные центроиды для начала\n",
    "np.random.seed(seed=42)\n",
    "centroids = np.random.normal(loc=0.0, scale=1.0, size=6)\n",
    "centroids = centroids.reshape((3, 2))\n",
    "\n",
    "cent_history = []\n",
    "cent_history.append(centroids)\n",
    "\n",
    "for i in range(3):\n",
    "    # Считаем расстояния от наблюдений до центроид\n",
    "    distances = cdist(X, centroids)\n",
    "    # Смотрим, до какой центроиде каждой точке ближе всего\n",
    "    labels = distances.argmin(axis=1)\n",
    "\n",
    "    # Положим в каждую новую центроиду геометрический центр её точек\n",
    "    centroids = centroids.copy()\n",
    "    centroids[0, :] = np.mean(X[labels == 0, :], axis=0)\n",
    "    centroids[1, :] = np.mean(X[labels == 1, :], axis=0)\n",
    "    centroids[2, :] = np.mean(X[labels == 2, :], axis=0)\n",
    "\n",
    "    cent_history.append(centroids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# А теперь нарисуем всю эту красоту\n",
    "plt.figure(figsize=(8, 8))\n",
    "for i in range(4):\n",
    "    distances = cdist(X, cent_history[i])\n",
    "    labels = distances.argmin(axis=1)\n",
    "\n",
    "    plt.subplot(2, 2, i + 1)\n",
    "    plt.plot(X[labels == 0, 0], X[labels == 0, 1], \"bo\", label=\"cluster #1\")\n",
    "    plt.plot(X[labels == 1, 0], X[labels == 1, 1], \"co\", label=\"cluster #2\")\n",
    "    plt.plot(X[labels == 2, 0], X[labels == 2, 1], \"mo\", label=\"cluster #3\")\n",
    "    plt.plot(cent_history[i][:, 0], cent_history[i][:, 1], \"rX\")\n",
    "    plt.legend(loc=0)\n",
    "    plt.title(\"Step {:}\".format(i + 1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Также стоит заметить, что хоть мы и рассматривали евклидово расстояние, алгоритм будет сходиться и в случае любой другой метрики, поэтому для различных задач кластеризации в зависимости от данных можно экспериментировать не только с количеством шагов или критерием сходимости, но и с метрикой, по которой мы считаем расстояния между точками и центроидами кластеров.\n",
    "\n",
    "Другой особенностью этого алгоритма является то, что он чувствителен к исходному положению центроид кластеров в пространстве. В такой ситуации спасает несколько последовательных запусков алгоритма с последующим усреднением полученных кластеров."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Выбор числа кластеров для kMeans\n",
    "\n",
    "В отличие от задачи классификации или регресии, в случае кластеризации сложнее выбрать критерий, с помощью которого было бы просто представить задачу кластеризации как задачу оптимизации.\n",
    "В случае kMeans распространен вот такой критерий – сумма квадратов расстояний от точек до центроидов кластеров, к которым они относятся.\n",
    "$$ J(C) = \\sum_{k=1}^K\\sum_{i~\\in~C_k} ||x_i - \\mu_k|| \\rightarrow \\min\\limits_C,$$\n",
    "\n",
    "здесь $C$ – множество кластеров мощности $K$, $\\mu_k$ – центроид кластера $C_k$.\n",
    "\n",
    "Понятно, что здравый смысл в этом есть: мы хотим, чтобы точки распологались кучно возле центров своих кластеров. Но вот незадача: минимум такого фнукционала будет достигаться тогда, когда кластеров столько же, сколько и точек (то есть каждая точка – это кластер из одного элемента).\n",
    "Для решения этого вопроса (выбора числа кластеров) часто пользуются такой эвристикой: выбирают то число кластеров, начиная с которого описанный функционал $ J(C) $ падает \"уже не так быстро\". Или более формально: \n",
    "$$ D(k) = \\frac{|J(C_k) - J(C_{k+1})|}{|J(C_{k-1}) - J(C_k)|}  \\rightarrow \\min\\limits_k $$\n",
    "\n",
    "Рассмотрим пример."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inertia = []\n",
    "for k in range(1, 8):\n",
    "    kmeans = KMeans(n_clusters=k, random_state=1).fit(X)\n",
    "    inertia.append(np.sqrt(kmeans.inertia_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(range(1, 8), inertia, marker=\"s\")\n",
    "plt.xlabel(\"$k$\")\n",
    "plt.ylabel(\"$J(C_k)$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что $J(C_k)$ падает сильно при увеличении числа кластеров с 1 до 2 и с 2 до 3 и уже не так сильно – при изменении $k$ с 3 до 4. Значит, в данной задаче оптимально задать 3 кластера.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Сложности\n",
    "\n",
    "Само по себе решение задачи K-means NP-трудное (NP-Hard, [статья](https://habrahabr.ru/post/132127/) \"Еще немного про P и NP\" на Хабре), и для размерности $d$, числа кластеров $k$ и числа точек $n$ решается за $O(n^{d k+1})$. Для решения такой боли часто используются эвристики, например MiniBatch K-means, который для обучения использует не весь датасет целиком, а лишь маленькие его порции (batch) и обновляет центроиды используя среднее за всю историю обновлений центроида от всех относящихся к нему точек. Сравнение обычного K-means и его MiniBatch имплементации можно посмотреть в [документации scikit-learn](http://scikit-learn.org/stable/auto_examples/cluster/plot_mini_batch_kmeans.html).\n",
    "\n",
    "[Реализация](http://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html) алгоритма в scikit-learn обладает массой удобных плюшек, таких как возможность задать количество запусков через параметр `n_init`, что даст более устойчивые центроиды для кластеров в случае скошенных данных. К тому же эти запуски можно делать параллельно, не жертвуя временем вычисления."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Affinity Propagation\n",
    "\n",
    "Ещё один пример алгоритма кластеризации. В отличие от алгоритма К-средних, данный подход не требует заранее определять число кластеров на которое мы хотим разбить наши данные. Основная идея алгоритма заключается в том, что нам хотелось бы, чтобы наши наблюдения кластеризовались в группы на основе того, как они \"общаются\" или насколько они похожи друг на друга.\n",
    "\n",
    "Заведём для этого какую-нибудь метрику \"похожести\", определяющуюся тем, что $s(x_i, x_j) > s(x_i, x_k)$ если наблюдение $x_i$ больше похоже на наблюдение $x_j$, чем на $x_k$. Простым примером такой похожести будет отрицательный квадрат расстояния $s(x_i, x_j) = - ||x_i - x_j||^{2}$.\n",
    "\n",
    "Теперь опишем сам процесс \"общения\" для этого заведём две матрицы, инициализируемые нулями, одна из которых $r_{i,k}$ будет описывать насколько хорошо $k$-тое наблюдение подходит для того, чтобы быть \"примером для подражания\" для $i$-того наблюдения, относительно всех остальных потенциальных \"примеров\", а вторая — $a_{i,k}$ будет описывать насколько правильным было бы для $i$-того наблюдения выбрать $k$-тое в качестве такого \"примера\". Звучит немного запутанно, но чуть дальше увидим пример \"на пальцах\".\n",
    "\n",
    "После этого данные матрицы обновляются по очереди по правилам:\n",
    "\n",
    "$$r_{i,k} \\leftarrow s_(x_i, x_k) - \\max_{k' \\neq k} \\left\\{ a_{i,k'} + s(x_i, x_k') \\right\\}$$\n",
    "\n",
    "$$a_{i,k} \\leftarrow \\min \\left( 0, r_{k,k} + \\sum_{i' \\not\\in \\{i,k\\}} \\max(0, r_{i',k}) \\right), \\ \\ \\  i \\neq k$$\n",
    "\n",
    "$$a_{k,k} \\leftarrow \\sum_{i' \\neq k} \\max(0, r_{i',k})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Спектральная кластеризация\n",
    "Спектральная кластеризация объединяет несколько описанных выше подходов, чтобы получить максимальное количество профита от сложных многообразий размерности меньшей исходного пространства. \n",
    "\n",
    "Для работы этого алгоритма нам потребуется определить матрицу похожести наблюдений (adjacency matrix). Можно это сделать таким же образом, как и для Affinity Propagation выше: $A_{i, j} = - ||x_i - x_j||^{2}$. Эта матрица также описывает полный граф с вершинами в наших наблюдениях и рёбрами между каждой парой наблюдений с весом, соответствующим степени похожести этих вершин. Для нашей выше выбранной метрики и точек, лежащих на плоскости, эта штука будет интуитивной и простой — две точки более похожи, если ребро между ними короче. Теперь нам бы хотелось разделить наш получившийся граф на две части так, чтобы получившиеся точки в двух графах были в общем больше похожи на другие точки внутри получившейся \"своей\" половины графа, чем на точки в \"другой\" половине. Формальное название такой задачи называется Normalized cuts problem и подробнее про это можно почитать [тут](http://people.eecs.berkeley.edu/~malik/papers/SM-ncut.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Алгомеративная кластеризация\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Наверное самый простой и понятный алгоритм кластеризации без фиксированного числа кластеров — агломеративная кластеризация. Интуиция у алгоритма очень простая: \n",
    "1. Начинаем с того, что высыпаем на каждую точку свой кластер\n",
    "2. Сортируем попарные расстояния между центрами кластеров по возрастанию\n",
    "3. Берём пару ближайших кластеров, склеиваем их в один и пересчитываем центр кластера\n",
    "4. Повторяем п. 2 и 3 до тех пор, пока все данные не склеятся в один кластер\n",
    "\n",
    "Сам процесс поиска ближайших кластеров может происходить с использованием разных методов объединения точек:\n",
    "1. Single linkage — минимум попарных расстояний между точками из двух кластеров\n",
    "$d(C_i, C_j) = min_{x_i \\in C_i, x_j \\in C_j} ||x_i - x_j||$\n",
    "2. Complete linkage — максимум попарных расстояний между точками из двух кластеров\n",
    "$d(C_i, C_j) = max_{x_i \\in C_i, x_j \\in C_j} ||x_i - x_j||$\n",
    "3. Average linkage — среднее попарных расстояний между точками из двух кластеров\n",
    "$d(C_i, C_j) = \\frac{1}{n_i n_j} \\sum_{x_i \\in C_i} \\sum_{x_j \\in C_j} ||x_i - x_j||$\n",
    "4. Centroid linkage — расстояние между центроидами двух кластеров\n",
    "$d(C_i, C_j) = ||\\mu_i - \\mu_j||$\n",
    "\n",
    "Профит первых трёх подходов по сравнению с четвёртым в том, что для них не нужно будет пересчитывать расстояния каждый раз после склеивания, что сильно снижает вычислительную сложность алгоритма.\n",
    "\n",
    "По итогам выполнения такого алгоритма можно также построить замечательное дерево склеивания кластеров и глядя на него определить, на каком этапе нам было бы оптимальнее всего остановить алгоритм. Либо воспользоваться тем же правилом локтя, что и в k-means.\n",
    "\n",
    "К счастью для нас в питоне уже есть замечательные инструменты для построения таких дендрограмм для агломеративной кластеризации. Рассмотрим на примере наших кластеров из K-means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster import hierarchy\n",
    "from scipy.spatial.distance import pdist\n",
    "\n",
    "X = np.zeros((150, 2))\n",
    "\n",
    "np.random.seed(seed=42)\n",
    "X[:50, 0] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "X[:50, 1] = np.random.normal(loc=0.0, scale=0.3, size=50)\n",
    "\n",
    "X[50:100, 0] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "X[50:100, 1] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "\n",
    "X[100:150, 0] = np.random.normal(loc=-1.0, scale=0.2, size=50)\n",
    "X[100:150, 1] = np.random.normal(loc=2.0, scale=0.5, size=50)\n",
    "\n",
    "distance_mat = pdist(\n",
    "    X\n",
    ")  # pdist посчитает нам верхний треугольник матрицы попарных расстояний\n",
    "\n",
    "Z = hierarchy.linkage(\n",
    "    distance_mat, \"single\"\n",
    ")  # linkage — реализация агломеративного алгоритма\n",
    "plt.figure(figsize=(10, 5))\n",
    "dn = hierarchy.dendrogram(Z, color_threshold=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Метрики качества кластеризации\n",
    "Задача оценки качества кластеризации является более сложной по сравнению с оценкой качества классификации. Во-первых, такие оценки не должны зависеть от самих значений меток, а только от самого разбиения выборки. Во-вторых, не всегда известны истинные метки объектов, поэтому также нужны оценки, позволяющие оценить качество кластеризации, используя только неразмеченную выборку.\n",
    "\n",
    "Выделяют *внешние* и *внутренние* метрики качества. Внешние используют информацию об истинном разбиении на кластеры, в то время как внутренние метрики не используют никакой внешней информации и оценивают качество кластеризации, основываясь только на наборе данных. Оптимальное число кластеров обычно определяют с использованием внутренних метрик.\n",
    "\n",
    "Все указанные ниже метрики реализованы в `sklearn.metrics`.\n",
    "\n",
    "**Adjusted Rand Index (ARI)**\n",
    "\n",
    "Предполагается, что известны истинные метки объектов. Данная мера не зависит от самих значений меток, а только от разбиения выборки на кластеры. Пусть $N$ - число объектов в выборке. Обозначим через $a$ - число пар объектов, имеющих одинаковые метки и находящихся в одном кластере, через $b$ - число пар объектов, имеющих различные метки и находящихся в разных кластерах. Тогда Rand Index это $$\\text{RI} = \\frac{2(a + b)}{n(n-1)}.$$ То есть это доля объектов, для которых эти разбиения (исходное и полученное в результате кластеризации) \"согласованы\". Rand Index (RI) выражает схожесть двух разных кластеризаций одной и той же выборки. Чтобы этот индекс давал значения близкие к нулю для случайных кластеризаций при любом $N$ и числе кластеров, необходимо нормировать его. Так определяется Adjusted Rand Index: $$\\text{ARI} = \\frac{\\text{RI} - E[\\text{RI}]}{\\max(\\text{RI}) - E[\\text{RI}]}.$$\n",
    "\n",
    "Эта мера симметрична, не зависит от значений и перестановок меток. Таким образом, данный индекс является мерой расстояния между различными разбиениями выборки. $\\text{ARI}$ принимает значения в диапазоне $[-1, 1]$. Отрицательные значения соответствуют \"независимым\" разбиениям на кластеры, значения, близкие к нулю, - случайным разбиениям, и положительные значения говорят о том, что два разбиения схожи (совпадают при $\\text{ARI} = 1$).\n",
    "\n",
    "**Adjusted Mutual Information (AMI)**\n",
    "\n",
    "Данная мера очень похожа на $\\text{ARI}$. Она также симетрична, не зависит от значений и перестановок меток. Определяется с использованием функции [энтропии](https://en.wikipedia.org/wiki/Entropy_(information_theory), интерпретируя разбиения выборки, как дискретные распределения (вероятность отнесения к кластеру равна доле объектов в нём). Индекс $MI$ определяется как  [взаимная информация](https://en.wikipedia.org/wiki/Mutual_information) для двух распределений, соответствующих разбиениям выборки на кластеры. Интуитивно, взаимная информация измеряет долю информации, общей для обоих разбиений: насколько информация об одном из них уменьшает неопределенность относительно другого.\n",
    "\n",
    "Аналогично $\\text{ARI}$ определяется индекс $\\text{AMI}$, позволяющий избавиться от роста индекса $MI$ с увеличением числа классов. Он принимает значения в диапазоне $[0, 1]$. Значения, близкие к нулю, говорят о независимости разбиений, а близкие к единице - об их схожести (совпадении при $\\text{AMI} = 1$).\n",
    "\n",
    "**Гомогенность, полнота, V-мера**\n",
    "\n",
    "Формально данные меры также определяются с использованием функций энтропии и условной энтропии, рассматривая разбиения выборки как дискретные распределения:$$h = 1 - \\frac{H(C\\mid K)}{H(C)}, c = 1 - \\frac{H(K\\mid C)}{H(K)},$$\n",
    "здесь $K$ - результат кластеризации, $C$ - истинное разбиение выборки на классы. Таким образом, $h$ измеряет, насколько каждый кластер состоит из объектов одного класса, а $c$ - насколько объекты одного класса относятся к одному кластеру. Эти меры не являются симметричными. Обе величины принимают значения в диапазоне $[0, 1]$, и большие значения соответствуют более точной кластеризации. Эти меры не являются нормализованными, как $\\text{ARI}$ или $\\text{AMI}$, и поэтому зависят от числа кластеров. Случайная кластеризация не будет давать нулевые показатели при большом числе классов и малом числе объектов. В этих случаях предпочтительнее использовать $\\text{ARI}$. Однако при числе объектов более 1000 и числе кластеров менее 10 данная проблема не так явно выражена и может быть проигнорирована.\n",
    "\n",
    "Для учёта обеих величин $h$ и $c$ одновременно вводится $V$-мера, как их среднее гармоническое:\n",
    "$$v = 2\\frac{hc}{h+c}.$$\n",
    "Она является симметричной и показывает, насколько две кластеризации схожи между собой.\n",
    "\n",
    "**Силуэт**\n",
    "\n",
    "В отличие от описанных выше метрик, данный коэффициент не предполагает знания истинных меток объектов, и позволяет оценить качество кластеризации, используя только саму (неразмеченную) выборку и результат кластеризации. Сначала силуэт определяется отдельно для каждого объекта. Обозначим через $a$ - среднее расстояние от данного объекта до объектов из того же кластера, через $b$ - среднее расстояние от данного объекта до объектов из ближайшего кластера (отличного от того, в котором лежит сам объект). Тогда силуэтом данного объекта называется величина: $$s = \\frac{b - a}{\\max(a, b)}.$$\n",
    "Силуэтом выборки называется средняя величина силуэта объектов данной выборки. Таким образом, силуэт показывает, насколько среднее расстояние до объектов своего кластера отличается от среднего расстояния до объектов других кластеров. Данная величина лежит в диапазоне $[-1, 1]$. Значения, близкие к -1, соответствуют плохим (разрозненным) кластеризациям, значения, близкие к нулю, говорят о том, что кластеры пересекаются и накладываются друг на друга, значения, близкие к 1, соответствуют \"плотным\" четко выделенным кластерам. Таким образом, чем больше силуэт, тем более четко выделены кластеры, и они представляют собой компактные, плотно сгруппированные облака точек.\n",
    "\n",
    "С помощью силуэта можно выбирать оптимальное число кластеров $k$ (если оно заранее неизвестно) - выбирается число кластеров, максимизирующее значение силуэта. В отличие от предыдущих метрик, силуэт зависит от формы кластеров, и достигает больших значений на более выпуклых кластерах, получаемых с помощью алгоритмов, основанных на восстановлении плотности распределения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "И напоследок давайте посмотрим на эти метрики для наших алгоритмов, запущенных на данных рукописных цифр MNIST:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn import datasets, metrics\n",
    "from sklearn.cluster import (\n",
    "    AffinityPropagation,\n",
    "    AgglomerativeClustering,\n",
    "    KMeans,\n",
    "    SpectralClustering,\n",
    ")\n",
    "\n",
    "data = datasets.load_digits()\n",
    "X, y = data.data, data.target\n",
    "\n",
    "algorithms = []\n",
    "algorithms.append(KMeans(n_clusters=10, random_state=1))\n",
    "algorithms.append(AffinityPropagation())\n",
    "algorithms.append(\n",
    "    SpectralClustering(n_clusters=10, random_state=1, affinity=\"nearest_neighbors\")\n",
    ")\n",
    "algorithms.append(AgglomerativeClustering(n_clusters=10))\n",
    "\n",
    "data = []\n",
    "for algo in algorithms:\n",
    "    algo.fit(X)\n",
    "    data.append(\n",
    "        (\n",
    "            {\n",
    "                \"ARI\": metrics.adjusted_rand_score(y, algo.labels_),\n",
    "                \"AMI\": metrics.adjusted_mutual_info_score(y, algo.labels_),\n",
    "                \"Homogenity\": metrics.homogeneity_score(y, algo.labels_),\n",
    "                \"Completeness\": metrics.completeness_score(y, algo.labels_),\n",
    "                \"V-measure\": metrics.v_measure_score(y, algo.labels_),\n",
    "                \"Silhouette\": metrics.silhouette_score(X, algo.labels_),\n",
    "            }\n",
    "        )\n",
    "    )\n",
    "\n",
    "results = pd.DataFrame(\n",
    "    data=data,\n",
    "    columns=[\"ARI\", \"AMI\", \"Homogenity\", \"Completeness\", \"V-measure\", \"Silhouette\"],\n",
    "    index=[\"K-means\", \"Affinity\", \"Spectral\", \"Agglomerative\"],\n",
    ")\n",
    "\n",
    "results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Полезные источники\n",
    "- [статья](https://habrahabr.ru/company/io/blog/265089/) \"Как подобрать платье с помощью метода главных компонент\"\n",
    "- [статья](https://habrahabr.ru/users/lexxo/topics/) \"Как работает метод главных компонент (PCA) на простом примере\"\n",
    "- [статья](https://habrahabr.ru/post/321216/) \"Интересные алгоритмы кластеризации, часть первая: Affinity propagation\"\n",
    "- [статья](https://habrahabr.ru/post/322034/) \"Интересные алгоритмы кластеризации, часть вторая: DBSCAN\"\n",
    "- [обзор](http://scikit-learn.org/stable/modules/clustering.html) методов кластеризации в пакете scikit-learn (англ.)\n",
    "- [Q&A](http://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues) Разбор PCA с интуицией и примерами (англ.)\n",
    "- [тетрадка](https://github.com/diefimov/MTH594_MachineLearning/blob/master/ipython/Lecture10.ipynb) про k-Means и EM-алгоритм в курсе Дмитрия Ефимова (англ.)\n",
    "- [конспект](https://github.com/esokolov/ml-course-hse/blob/master/2016-fall/lecture-notes/lecture12-unsupervised.pdf) \"Обучение без учителя\" в курсе Евгения Соколова"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}