{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "

\n", "## Открытый курс по машинному обучению\n", "Автор материала: Елена Девятайкина (@elenadevyataykina)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Оценка плотности, используя Gaussian Mixture Models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "В данном курсе были рассмаотрены следующие алгоритмы кластеризации: \n", "- K-means\n", "\n", "- Affinity Propagation\n", "\n", "- Спектральная кластеризация\n", "\n", "- Агломеративная кластеризация\n", "\n", "Здесь будет рассмотрен EM-алгоритм, реализованный в модуле `sklearn.mixture.GMM`.\n", "\n", "**Что это такое?**\n", "\n", "Это алгоритм, используемый для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных.\n", "\n", "**Как он работает?**\n", "\n", "Оценивая максимальное правдоподобие, EM-алгоритм создает модель, которая назначает метки класса точкам данных.\n", "\n", "EM-алгоритм начинает с того, что пытается сделать вывод на основании параметров модели.\n", "\n", "Затем следует итерационный трехшаговый процесс:\n", "\n", "1. E-шаг: На этом шаге на основании параметров модели вычисляются вероятность принадлежности каждой точки данных к кластеру.\n", "2. М-шаг: Обновляет параметры модели в соответствии с кластерным распределением, проведенным на шаге E.\n", "3. Предыдущие два шага повторяются до тех пор, пока параметры модели и кластерное распределение не уравняются.\n", "\n", "Часто данный алгоритм используют для разделения смеси гаусиан (GMM). То есть в результате GMM-EM получаем полноценную смесь Гауссианов полученную в локальном максимуме правдоподобия.\n", "\n", "Далее рассмаотрим GMM более подробно.\n", "\n", "\n", "**GMM** - один из алгоритмов кластеризации наряду с иерархической кластеризацией (кластеризация на основе связности), k-means (кластеризация на основе центров), DBSCAN (кластеризация на основе плотности). В основе данного алгоритма лежат распределения. \n", "\n", "**Гауссова смесь распределений (GMM — Gaussian mixture models)** - статистическая модель для представления нормально распределенных подвыборок внутри общей выборки. Гауссова смесь распределений параметризируется двумя типами значений — смесь весов компонентов и средних компонентов или ковариаций (для многомерного случая) и др. Если количество компонентов известно, техника, чаще всего используемая для оценки параметров смеси распределений — ЕМ-алгоритм.\n", "\n", "Метод GMM непосредственно вытекает из теоремы, гласящей, что любая функция плотности вероятности может быть представлена как взвешенная сумма нормальных распределений\n", "\n", "$$\n", "p = \\sum_{j=1}^k w_i \\varphi (x; \\theta_j),$$\n", "\n", "$$\n", "\\sum_{j=1}^k w_j = 1,$$\n", "где $\\varphi (x; \\theta_j)$ - функция распределения многомерного аргумента $x$ с параметрами $\\theta_j,$ \n", "\n", "$$\\varphi (x; \\theta_j) = p(x | \\mu_j, R_j) = \\frac{1}{{2\\pi}^{n/2}{|R_j|}^{1/2}}e^{-\\frac{1}{2}{(x-\\mu_j)}^TR_j^{-1}(x-\\mu_j)}, x \\in \\mathbb R^n$$\n", "\n", "$w_j$ - ее вес, $k$ - количество компонент в смеси. Здесь $n$ - размерность пространства признаков, $\\mu_j \\in \\mathbb R^n$ - вектор математического ожидания $j$-й компоненты смеси, $R_j \\in \\mathbb R^{n*n}$ - ковариационная матрица.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "\n", "warnings.simplefilter('ignore')\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy import stats\n", "from sklearn.mixture import GMM\n", "from sklearn.neighbors import KernelDensity\n", "\n", "import seaborn as sns; sns.set()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Введение в GMM" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Создадим одномерный набор данных с азаднным распределением и построим гистограмму" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(2)\n", "x = np.concatenate([np.random.normal(0, 2, 2000),\n", " np.random.normal(5, 5, 2000),\n", " np.random.normal(3, 0.5, 600)])\n", "plt.hist(x, 80, normed=True)\n", "plt.xlim(-10, 20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "С помощью GMM отобразим плотность распределения" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "X = x[:, np.newaxis]\n", "clf = GMM(4, n_iter=500, random_state=3).fit(X)\n", "xpdf = np.linspace(-10, 20, 1000)\n", "density = np.exp(clf.score(xpdf[:, np.newaxis]))\n", "\n", "plt.hist(x, 80, normed=True, alpha=0.5)\n", "plt.plot(xpdf, density, '-r')\n", "plt.xlim(-10, 20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Посмотрим на полученные атрибуты GMM.\n", "\n", "Среднее значение каждой компоненты:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "clf.means_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ковариация каждой компоненты" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "clf.covars_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Вес каждой компоненты" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "clf.weights_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Так как в GMM задано 4 компоненты, построим все плотности распределения" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.hist(x, 80, normed=True, alpha=0.3)\n", "plt.plot(xpdf, density, '-r')\n", "\n", "for i in range(clf.n_components):\n", " pdf = clf.weights_[i] * stats.norm(clf.means_[i, 0],\n", " np.sqrt(clf.covars_[i, 0])).pdf(xpdf)\n", " plt.fill(xpdf, pdf, facecolor='gray',\n", " edgecolor='none', alpha=0.3)\n", "plt.xlim(-10, 20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Для построенной модели можно использовать одно из нескольких средних, чтобы оценить, насколько хорошо он подходит данным.\n", "Обычно для этого используется **Информационный критерий Акаике (AIC)** или **Байесовский информационный критерий (BIC)**.\n", "\n", "**Информационный критерий Акаике** - критерий, применяющийся исключительно для выбора из нескольких статистических моделей.\n", "\n", "$$aic = 2k - 2l,$$\n", "\n", "где $k$ - число параметров в статистической модели, $l$ - значение логарифмической функции правдоподобия построенной модели.\n", "Критерий не только вознаграждает за качество приближения, но и штрафует за использование излишнего количества параметров модели. Считается, что наилучшей будет модель с наименьшим значением критерия aic.\n", "\n", "**Байесовский информационный критерий** - критерий, разработанный исходя из баесовского подхода, модификация aic.\n", "\n", "$$bic = k ln(n) - 2l$$\n", "\n", "Данный критерий налагает больший штраф на увеличение количества параметров по сравнению с aic." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Посмотрим на значение aic, bic на заданной выборке" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(clf.bic(X))\n", "print(clf.aic(X))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Давайте рассмотрим их как функцию числа гауссиан:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_estimators = np.arange(1, 10)\n", "clfs = [GMM(n, n_iter=1000).fit(X) for n in n_estimators]\n", "bics = [clf.bic(X) for clf in clfs]\n", "aics = [clf.aic(X) for clf in clfs]\n", "\n", "plt.plot(n_estimators, bics, label='BIC')\n", "plt.plot(n_estimators, aics, label='AIC')\n", "plt.xlabel('Number of components')\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Судя по графику можно сказать, что для aic предпочтительнее 4 компоненты, а для bic - 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Так как GMM является Generative Model(вероятностная модель для генерации набора данных). С помощью даннной модели можно обнаружить выбросы (Рассчитывается вероятность каждой точки под моделью и чем ниже $p(x)$, тем вероятнее, что данная точка является выбросом. Есть подробная статья, связанная с этой темой - https://arxiv.org/pdf/1706.02690.pdf)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(0)\n", "\n", "# Добавим к выборке 20 выбросов\n", "true_outliers = np.sort(np.random.randint(0, len(x), 20))\n", "y = x.copy()\n", "y[true_outliers] += 50 * np.random.randn(20)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "clf = GMM(4, n_iter=500, random_state=0).fit(y[:, np.newaxis])\n", "xpdf = np.linspace(-10, 20, 1000)\n", "density_noise = np.exp(clf.score(xpdf[:, np.newaxis]))\n", "\n", "plt.hist(y, 80, normed=True, alpha=0.5)\n", "plt.plot(xpdf, density_noise, '-r')\n", "plt.xlim(-15, 30);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Далее оценим логарифмическую вероятность каждой точки под моделью" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "log_likelihood = clf.score_samples(y[:, np.newaxis])[0]\n", "plt.plot(y, log_likelihood, '.k');" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "detected_outliers = np.where(log_likelihood < -9)[0]\n", "\n", "print(\"Все выбросы:\")\n", "print(true_outliers)\n", "print(\"\\nОбнаруженные выбросы:\")\n", "print(detected_outliers)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Как можно заметить, алгоритм пропустил несколько точек (те, которые находятся в середине распределения).\n", "\n", "Вот пропущенные выбросы:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "set(true_outliers) - set(detected_outliers)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "И вот значения, которые были неверно определены как выбросы" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "set(detected_outliers) - set(true_outliers)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Необходимо отметить, что здесь используются одномерные данные, но GMM также может делать обобщения на несколько измерений." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Другие оценки плотности" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Важным средством оценки плотности является оценка плотности ядра, реализованная в `sklearn.neighbors.KernelDensity`. Можно это рассматривать как обобщение GMM. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kde = KernelDensity(0.15).fit(x[:, None])\n", "density_kde = np.exp(kde.score_samples(xpdf[:, None]))\n", "\n", "plt.hist(x, 80, normed=True, alpha=0.5)\n", "plt.plot(xpdf, density, '-b', label='GMM')\n", "plt.plot(xpdf, density_kde, '-r', label='KDE')\n", "plt.xlim(-10, 20)\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "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.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }