Оценка плотности ядра в Python с помощью Scikit-Learn

Вступление

Эта статья представляет собой введение в оценку плотности ядра с использованием библиотеки машинного обучения Python scikit-learn .

Оценка плотности ядра (KDE) – это непараметрический метод оценки функции плотности вероятности заданной случайной величины. Он также упоминается под своим традиционным названием “Окно Парцена-Розенблатта” в честь его первооткрывателей.

Учитывая выборку независимых, тождественно распределенных (i.i.d) наблюдений $(x_1,x_2,\ldots,x_n)$ случайной величины из неизвестного исходного распределения, оценка плотности ядра задается:

$$ p(x) = \frac{1}{nh}}^{n}K(\frac{x-x_j}{h}) $$

где $K(a)$ – функция ядра, а $h$ – параметр сглаживания, также называемый полосой пропускания. Различные ядра обсуждаются далее в этой статье, но просто чтобы понять математику, давайте рассмотрим простой пример.

Пример Вычисления

Предположим, у нас есть точки выборки [-2,-1,0,1,2] , с линейным ядром, заданным:-\frac{|a|}{h}\) и\).

Подключите вышеизложенное к формуле для $p(x)$:

$$ p(0) = \frac{1}{(5)(10)} ( 0.8+0.9+1+0.9+0.8.088 $$

Оценка Плотности Ядра С Помощью Python

Хотя существует несколько способов вычисления оценки плотности ядра в Python, мы будем использовать для этой цели популярную библиотеку машинного обучения scikit-learn . Импортируйте в свой код следующие библиотеки:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV

Синтетические Данные

Чтобы продемонстрировать оценку плотности ядра, синтетические данные генерируются из двух различных типов распределений. Одно из них-асимметричное логарифмически нормальное распределение, а другое – гауссово распределение. Следующая функция возвращает 2000 точек данных:

def generate_data(seed=17):
    # Fix the seed to reproduce the results
    rand = np.random.RandomState(seed)
    x = []
    dat = rand.lognormal(0, 0.3, 1000)
    x = np.concatenate((x, dat))
    dat = rand.normal(3, 1, 1000)
    x = np.concatenate((x, dat))
    return x

Приведенный ниже код хранит точки в x_train . Мы можем либо построить диаграмму рассеяния этих точек вдоль оси y, либо сгенерировать гистограмму этих точек.

x_train = generate_data()[:, np.newaxis]
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))
plt.subplot(121)
plt.scatter(np.arange(len(x_train)), x_train, c='red')
plt.xlabel('Sample no.')
plt.ylabel('Value')
plt.title('Scatter plot')
plt.subplot(122)
plt.hist(x_train, bins=50)
plt.title('Histogram')
fig.subplots_adjust(wspace=.3)
plt.show()

Использование плотности ядра Scikit-Learn

Чтобы найти форму оцениваемой функции плотности, мы можем сгенерировать набор точек, равноудаленных друг от друга, и оценить плотность ядра в каждой точке. Тестовые баллы задаются:

x_test = np.linspace(-1, 7, 2000)[:, np.newaxis]

Теперь мы создадим объект Kernel Density и с помощью метода fit() найдем оценку каждого образца, как показано в приведенном ниже коде. Метод Kernel Density() использует два параметра по умолчанию, т. е. kernel=gaussian и bandwidth=1 .

model = KernelDensity()
model.fit(x_train)
log_dens = model.score_samples(x_test)

Форму распределения можно увидеть, построив график оценки плотности для каждой точки, как показано ниже:

plt.fill(x_test, np.exp(log_dens), c='cyan')
plt.show()

Понимание параметра Пропускной способности

Предыдущий пример представляет собой не очень впечатляющую оценку функции плотности, приписываемую в основном параметрам по умолчанию. Давайте поэкспериментируем с различными значениями пропускной способности, чтобы увидеть, как это влияет на оценку плотности.

bandwidths = [0.01, 0.05, 0.1, 0.5, 1, 4]
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(10, 7))
plt_ind = np.arange(6) + 231

for b, ind in zip(bandwidths, plt_ind):
    kde_model = KernelDensity(kernel='gaussian', bandwidth=b)
    kde_model.fit(x_train)
    score = kde_model.score_samples(x_test)
    plt.subplot(ind)
    plt.fill(x_test, np.exp(score), c='cyan')
    plt.title("h="+str(b))

fig.subplots_adjust(hspace=0.5, wspace=.3)
plt.show()

Мы можем ясно видеть, что увеличение полосы пропускания приводит к более гладкой оценке. Очень малые значения полосы пропускания приводят к колючим и нервным кривым, в то время как очень высокие значения приводят к очень обобщенной гладкой кривой, которая пропускает важные детали. Важно выбрать сбалансированное значение для этого параметра.

Настройка параметра Пропускной способности

Библиотека scikit-learn позволяет настраивать параметр bandwidth с помощью перекрестной проверки и возвращает значение параметра, которое максимизирует логарифмическую вероятность данных. Функция , которую мы можем использовать для достижения этой цели, – это GridSearchCV () , которая требует различных значений параметра bandwidth .

bandwidth = np.arange(0.05, 2, .05)
kde = KernelDensity(kernel='gaussian')
grid = GridSearchCV(kde, {'bandwidth': bandwidth})
grid.fit(x_train)

Лучшая модель может быть получена с помощью поля best_estimator_ объекта GridSearchCV .

Давайте рассмотрим оптимальную оценку плотности ядра с помощью гауссовского ядра и выведем также значение пропускной способности:

kde = grid.best_estimator_
log_dens = kde.score_samples(x_test)
plt.fill(x_test, np.exp(log_dens), c='green')
plt.title('Optimal estimate with Gaussian kernel')
plt.show()
print("optimal bandwidth: " + "{:.2f}".format(kde.bandwidth))

optimal bandwidth: 0.15

Теперь эта оценка плотности, по-видимому, очень хорошо моделирует данные. Первая половина графика согласуется с логнормальным распределением, а вторая половина графика довольно хорошо моделирует нормальное распределение.

Различные ядра для оценки плотности

scikit-learn позволяет оценивать плотность ядра с использованием различных функций ядра:

ядро : $K(a;h) \propto \cos (\frac{\pi a}{2h}) \text { if } |a| < h $
ядро : $K(a;h) \propto 1 – \frac{a^2}{h^2}$
ядро : $K(a;h) \propto \exp (-\frac{|a|}{h})$
ядро : $K(a;h) \propto \exp(-\frac{a^2}{2h^2})$
ядро : $K(a;h) \propto 1 – \frac{|a|}{h} \text { if } |a| < h $
ядро : $K(a;h) \propto 1 \text { if } |a| < h $

Простой способ понять, как работают эти ядра, – построить их график. Это означает построение модели, используя выборку только одного значения, например 0. Затем оцените плотность всех точек вокруг нуля и постройте график плотности вдоль оси y. Приведенный ниже код показывает весь процесс:

kernels = ['cosine', 'epanechnikov', 'exponential', 'gaussian', 'linear', 'tophat']
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(10, 7))
plt_ind = np.arange(6) + 231

for k, ind in zip(kernels, plt_ind):
    kde_model = KernelDensity(kernel=k)
    kde_model.fit([[0]])
    score = kde_model.score_samples(np.arange(-2, 2, 0.1)[:, None])
    plt.subplot(ind)
    plt.fill(np.arange(-2, 2, 0.1)[:, None], np.exp(score), c='blue')
    plt.title(k)

fig.subplots_adjust(hspace=0.5, wspace=.3)
plt.show()

Эксперименты С Различными Ядрами

Давайте поэкспериментируем с различными ядрами и посмотрим, как они оценивают функцию плотности вероятности для наших синтетических данных.

Мы можем использовать GridSearchCV() , как и раньше, чтобы найти оптимальное значение bandwidth . Однако для cosine , linear и top hat kernels GridSearchCV() может выдавать предупреждение во время выполнения из-за некоторых оценок, приводящих к значениям - inf . Один из возможных способов решения этой проблемы-написать пользовательскую функцию подсчета очков для GridSearchCV() .

В приведенном ниже коде -in баллы для тестовых точек опущены в пользовательской функции my_scores() и возвращается среднее значение. Это не обязательно лучшая схема для обработки значений -in score, и в зависимости от рассматриваемых данных может быть принята какая-то другая стратегия.

def my_scores(estimator, X):
    scores = estimator.score_samples(X)
    # Remove -inf
    scores = scores[scores != float('-inf')]
    # Return the mean values
    return np.mean(scores)

kernels = ['cosine', 'epanechnikov', 'exponential', 'gaussian', 'linear', 'tophat']
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(10, 7))
plt_ind = np.arange(6) + 231
h_vals = np.arange(0.05, 1, .1)

for k, ind in zip(kernels, plt_ind):
    grid = GridSearchCV(KernelDensity(kernel=k),
                        {'bandwidth': h_vals},
                        scoring=my_scores)
    grid.fit(x_train)
    kde = grid.best_estimator_
    log_dens = kde.score_samples(x_test)
    plt.subplot(ind)
    plt.fill(x_test, np.exp(log_dens), c='cyan')
    plt.title(k + " h=" + "{:.2f}".format(kde.bandwidth))

fig.subplots_adjust(hspace=.5, wspace=.3)
plt.show()

Окончательная Оптимизированная Модель

Приведенный выше пример показывает, как различные ядра оценивают плотность по-разному. Последний шаг-настроить GridSearchCV() так, чтобы он не только обнаружил оптимальную пропускную способность, но и оптимальное ядро для наших данных примера. Вот окончательный код, который также выводит окончательную оценку плотности и ее настроенные параметры в заголовке графика:

grid = GridSearchCV(KernelDensity(),
                    {'bandwidth': h_vals, 'kernel': kernels},
                    scoring=my_scores)
grid.fit(x_train)
best_kde = grid.best_estimator_
log_dens = best_kde.score_samples(x_test)
plt.fill(x_test, np.exp(log_dens), c='green')
plt.title("Best Kernel: " + best_kde.kernel+" h="+"{:.2f}".format(best_kde.bandwidth))
plt.show()

Вывод

Оценка плотности ядра с использованием библиотеки scikit-learn ‘ s library sklearn.neighbors обсуждалась в этой статье. Примеры приведены для одномерных данных, однако они также могут быть применены к данным с несколькими измерениями.

Будучи интуитивно понятным и простым способом оценки плотности для распределений неизвестных источников, специалист по обработке данных должен использовать его с осторожностью, поскольку проклятие размерности может значительно замедлить его.