<jupyter_start><jupyter_text>Открытый курс по машинному обучениюАвтор материала: <NAME> (@elenadevyataykina)  Оценка плотности, используя Gaussian Mixture Models В данном курсе были рассмаотрены следующие алгоритмы кластеризации: - K-means- Affinity Propagation- Спектральная кластеризация- Агломеративная кластеризацияЗдесь будет рассмотрен EM-алгоритм, реализованный в модуле `sklearn.mixture.GMM`.**Что это такое?**Это алгоритм, используемый для нахождения оценок максимального правдоподобия параметров вероятностных моделей, в случае, когда модель зависит от некоторых скрытых переменных.**Как он работает?**Оценивая максимальное правдоподобие, EM-алгоритм создает модель, которая назначает метки класса точкам данных.EM-алгоритм начинает с того, что пытается сделать вывод на основании параметров модели.Затем следует итерационный трехшаговый процесс:1. E-шаг: На этом шаге на основании параметров модели вычисляются вероятность принадлежности каждой точки данных к кластеру.2. М-шаг: Обновляет параметры модели в соответствии с кластерным распределением, проведенным на шаге E.3. Предыдущие два шага повторяются до тех пор, пока параметры модели и кластерное распределение не уравняются.Часто данный алгоритм используют для разделения смеси гаусиан (GMM). То есть в результате GMM-EM получаем полноценную смесь Гауссианов полученную в локальном максимуме правдоподобия.Далее рассмаотрим GMM более подробно.**GMM** - один из алгоритмов кластеризации наряду с иерархической кластеризацией (кластеризация на основе связности), k-means (кластеризация на основе центров), DBSCAN (кластеризация на основе плотности). В основе данного алгоритма лежат распределения. **Гауссова смесь распределений (GMM — Gaussian mixture models)** - статистическая модель для представления нормально распределенных подвыборок внутри общей выборки. Гауссова смесь распределений параметризируется двумя типами значений — смесь весов компонентов и средних компонентов или ковариаций (для многомерного случая) и др. Если количество компонентов известно, техника, чаще всего используемая для оценки параметров смеси распределений — ЕМ-алгоритм.Метод GMM непосредственно вытекает из теоремы, гласящей, что любая функция плотности вероятности может быть представлена как взвешенная сумма нормальных распределений$$p = \sum_{j=1}^k w_i \varphi (x; \theta_j),$$$$\sum_{j=1}^k w_j = 1,$$где $\varphi (x; \theta_j)$ - функция распределения многомерного аргумента $x$ с параметрами $\theta_j,$ $$\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$$$w_j$ - ее вес, $k$ - количество компонент в смеси. Здесь $n$ - размерность пространства признаков, $\mu_j \in \mathbb R^n$ - вектор математического ожидания $j$-й компоненты смеси, $R_j \in \mathbb R^{n*n}$ - ковариационная матрица.<jupyter_code>import warnings
warnings.simplefilter('ignore')

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

import seaborn as sns; sns.set()

from sklearn.mixture import GMM
from sklearn.neighbors import KernelDensity<jupyter_output><empty_output><jupyter_text>Введение в GMM Создадим одномерный набор данных с азаднным распределением и построим гистограмму<jupyter_code>np.random.seed(2)
x = np.concatenate([np.random.normal(0, 2, 2000),
                    np.random.normal(5, 5, 2000),
                    np.random.normal(3, 0.5, 600)])
plt.hist(x, 80, normed=True)
plt.xlim(-10, 20);<jupyter_output><empty_output><jupyter_text>С помощью GMM отобразим плотность распределения<jupyter_code>X = x[:, np.newaxis]
clf = GMM(4, n_iter=500, random_state=3).fit(X)
xpdf = np.linspace(-10, 20, 1000)
density = np.exp(clf.score(xpdf[:, np.newaxis]))

plt.hist(x, 80, normed=True, alpha=0.5)
plt.plot(xpdf, density, '-r')
plt.xlim(-10, 20);<jupyter_output><empty_output><jupyter_text>Посмотрим на полученные атрибуты GMM.Среднее значение каждой компоненты:<jupyter_code>clf.means_<jupyter_output><empty_output><jupyter_text>Ковариация каждой компоненты<jupyter_code>clf.covars_<jupyter_output><empty_output><jupyter_text>Вес каждой компоненты<jupyter_code>clf.weights_<jupyter_output><empty_output><jupyter_text>Так как в GMM задано 4 компоненты, построим все плотности распределения<jupyter_code>plt.hist(x, 80, normed=True, alpha=0.3)
plt.plot(xpdf, density, '-r')

for i in range(clf.n_components):
    pdf = clf.weights_[i] * stats.norm(clf.means_[i, 0],
                                       np.sqrt(clf.covars_[i, 0])).pdf(xpdf)
    plt.fill(xpdf, pdf, facecolor='gray',
             edgecolor='none', alpha=0.3)
plt.xlim(-10, 20);<jupyter_output><empty_output><jupyter_text>Для построенной модели можно использовать одно из нескольких средних, чтобы оценить, насколько хорошо он подходит данным.Обычно для этого используется **Информационный критерий Акаике (AIC)** или **Байесовский информационный критерий (BIC)**.**Информационный критерий Акаике** - критерий, применяющийся исключительно для выбора из нескольких статистических моделей.$$aic = 2k - 2l,$$где $k$ - число параметров в статистической модели, $l$ - значение логарифмической функции правдоподобия построенной модели.Критерий не только вознаграждает за качество приближения, но и штрафует за использование излишнего количества параметров модели. Считается, что наилучшей будет модель с наименьшим значением критерия aic.**Байесовский информационный критерий** - критерий, разработанный исходя из баесовского подхода, модификация aic.$$bic = k ln(n) - 2l$$Данный критерий налагает больший штраф на увеличение количества параметров по сравнению с aic. Посмотрим на значение aic, bic на заданной выборке<jupyter_code>print(clf.bic(X))
print(clf.aic(X))<jupyter_output>25911.1937804
25840.421853<jupyter_text>Давайте рассмотрим их как функцию числа гауссиан:<jupyter_code>n_estimators = np.arange(1, 10)
clfs = [GMM(n, n_iter=1000).fit(X) for n in n_estimators]
bics = [clf.bic(X) for clf in clfs]
aics = [clf.aic(X) for clf in clfs]

plt.plot(n_estimators, bics, label='BIC')
plt.plot(n_estimators, aics, label='AIC')
plt.xlabel('Number of components')
plt.legend();<jupyter_output><empty_output><jupyter_text>Судя по графику можно сказать, что для aic предпочтительнее 4 компоненты, а для bic - 2 Так как GMM является Generative Model(вероятностная модель для генерации набора данных). С помощью даннной модели можно обнаружить выбросы (Рассчитывается вероятность каждой точки под моделью и чем ниже $p(x)$, тем вероятнее, что данная точка является выбросом. Есть подробная статья, связанная с этой темой - https://arxiv.org/pdf/1706.02690.pdf)<jupyter_code>np.random.seed(0)

# Добавим к выборке 20 выбросов
true_outliers = np.sort(np.random.randint(0, len(x), 20))
y = x.copy()
y[true_outliers] += 50 * np.random.randn(20)
clf = GMM(4, n_iter=500, random_state=0).fit(y[:, np.newaxis])
xpdf = np.linspace(-10, 20, 1000)
density_noise = np.exp(clf.score(xpdf[:, np.newaxis]))

plt.hist(y, 80, normed=True, alpha=0.5)
plt.plot(xpdf, density_noise, '-r')
plt.xlim(-15, 30);<jupyter_output><empty_output><jupyter_text>Далее оценим логарифмическую вероятность каждой точки под моделью<jupyter_code>log_likelihood = clf.score_samples(y[:, np.newaxis])[0]
plt.plot(y, log_likelihood, '.k');
detected_outliers = np.where(log_likelihood < -9)[0]

print("Все выбросы:")
print(true_outliers)
print("\nОбнаруженные выбросы:")
print(detected_outliers)<jupyter_output>Все выбросы:
[  99  537  705 1033 1653 1701 1871 2046 2135 2163 2222 2496 2599 2607 2732
 2893 2897 3264 3468 4373]

Обнаруженные выбросы:
[  99  537  705 1653 2046 2135 2163 2222 2496 2732 2893 2897 3067 3173 3253
 3468 3483 4373]<jupyter_text>Как можно заметить, алгоритм пропустил несколько точек (те, которые находятся в середине распределения).Вот пропущенные выбросы:<jupyter_code>set(true_outliers) - set(detected_outliers)<jupyter_output><empty_output><jupyter_text>И вот значения, которые были неверно определены как выбросы<jupyter_code>set(detected_outliers) - set(true_outliers)<jupyter_output><empty_output><jupyter_text>Необходимо отметить, что здесь используются одномерные данные, но GMM также может делать обобщения на несколько измерений.  Другие оценки плотности Важным средством оценки плотности является оценка плотности ядра, реализованная в `sklearn.neighbors.KernelDensity`. Можно это рассматривать как обобщение GMM.<jupyter_code>kde = KernelDensity(0.15).fit(x[:, None])
density_kde = np.exp(kde.score_samples(xpdf[:, None]))

plt.hist(x, 80, normed=True, alpha=0.5)
plt.plot(xpdf, density, '-b', label='GMM')
plt.plot(xpdf, density_kde, '-r', label='KDE')
plt.xlim(-10, 20)
plt.legend();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1046189376443418, 'fraction_numerical': 0.03533487297921478, 'mean_word_length': 6.187551867219917, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 56, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '46751', 'license': 'MIT', 'n_tokens_mistral': 3768, 'n_tokens_neox': 4039, 'n_words': 823, 'path': 'jupyter/tutorials/gmm_tutorial_edevyataykina.ipynb', 'repo_name': 'ivan-magda/mlcourse_open_homeworks', 'size': '8671'}