В статистике оценка плотности ядра (KDE) - это применение сглаживания ядра для оценки плотности вероятности, т.е. непараметрический метод оценки функции плотности вероятности случайной величины на основе ядер в качестве весов. KDE - это фундаментальная задача сглаживания данных, в которой выводы о совокупности делаются на основе конечной выборки данных. В некоторых областях, таких как обработка сигналов и эконометрика, его также называют оконным методом Парзена-Розенблатта в честь Парзена и Розенблатта,которым обычно приписывают независимое создание его в его нынешнем виде.
Учитывая выборку независимых, одинаково распределенных (i.i.d) наблюдений (x1,x2,....,xn) случайной величины из неизвестного исходного распределения, оценка плотности ядра определяется как:
где К(а) является функцией ядра и h — параметр сглаживания, также называемый пропускной способностью
Оценка плотности ядра с использованием Python
Хотя в Python существует несколько способов вычисления оценки плотности ядра, мы будем использовать для этой цели популярную библиотеку машинного обучения scikit-learn.
(.env) boris@boris-All-Series:~/KDEPYTHON$ cat kdePython.py
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
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 = 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()
x_test = np.linspace(-1, 7, 2000)[:, np.newaxis]
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 позволяет настраивать параметр пропускной способности с помощью перекрестной проверки и возвращает значение параметра, которое максимизирует логарифмическую вероятность данных. Для этого мы можем использовать функцию GridSearchCV(),
для которой требуются разные значения параметра пропускной способности.
"""
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')
print("optimal bandwidth: " + "{:.2f}".format(kde.bandwidth))
plt.show()
"""
Простой способ понять, как работают эти ядра, — нарисовать их. Это означает построениемодели с использованием выборки только одного значения, например, 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()
"""
Эксперименты с разными ядрами
"""
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()
"""
The Final Optimized Model
"""
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()
No comments:
Post a Comment