Estimación de la densidad del kernel en Python usando Scikit-Learn

Introducción a la estimación de la densidad del núcleo utilizando scikit-learn. Ejemplos del uso de diferentes parámetros de kernel y ancho de banda para la optimización.

Introducción

Este artículo es una introducción a la estimación de la densidad del kernel utilizando la biblioteca de aprendizaje automático scikit-learn de Python.

La estimación de la densidad del kernel (KDE) es un método no paramétrico para estimar la función de densidad de probabilidad de una variable aleatoria determinada. También se le conoce por su nombre tradicional, el método Ventana de Parzen-Rosenblatt, en honor a sus descubridores.

Dada una muestra de observaciones independientes e idénticamente distribuidas (i.i.d) \((x_1,x_2,\ldots,x_n)\) de una variable aleatoria de una fuente de distribución desconocida, la estimación de la densidad del kernel viene dada por:

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

donde \(K(a)\) es la función kernel y \(h\) es el parámetro de suavizado, también llamado ancho de banda. Varios núcleos se discuten más adelante en este artículo, pero solo para entender las matemáticas, echemos un vistazo a un ejemplo simple.

Cálculo de ejemplo {#cálculo de ejemplo}

Supongamos que tenemos los puntos muestrales [-2,-1,0,1,2], con un kernel lineal dado por: \(K(a)= 1-\frac{|a| {h}\) y \(h=10\).

$$\begin{matriz} x_{j} & = & \lbrack & {- 2} & {- 1} & 0 & 1 & 2 & \rbrack \ {|0 - x_{j}|} & = & \lbrack & 2 & 1 & 0 & 1 & 2 & \rbrack \ {|\frac{0 - x_{j}}{h}|} & = & \lbrack & 0.2 & 0.1 & 0 & 0.1 & 0.2 & \rbrack \ {K(|\frac{0 - x_{j}}{h}|)} & = & \lbrack & 0.8 & 0.9 & 1 & 0.9 & 0.8 & \rbrack \ \end{matriz}$$

Introduce lo anterior en la fórmula para \(p(x)\):

$$
p(0) = \frac{1}{(5)(10)} ( 0,8+0,9+1+0,9+0,8 ) = 0,088
$$

Estimación de la densidad del kernel mediante Python

Si bien hay varias formas de calcular la estimación de la densidad del kernel en Python, usaremos la popular biblioteca de aprendizaje automático scikit-learn para este propósito. Importe las siguientes bibliotecas en su código:

1
2
3
4
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV

Datos sintéticos

Para demostrar la estimación de la densidad del núcleo, se generan datos sintéticos a partir de dos tipos diferentes de distribuciones. Una es una distribución logarítmica normal asimétrica y la otra es una distribución gaussiana. La siguiente función devuelve 2000 puntos de datos:

1
2
3
4
5
6
7
8
9
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

El siguiente código almacena los puntos en x_train. Podemos hacer un gráfico de dispersión de estos puntos a lo largo del eje y o podemos generar un histograma de estos puntos.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
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()

scatter plot and histogram estimation

Uso de KernelDensity {#usingscikitlearnskerneldensity} de Scikit-Learn

Para encontrar la forma de la función de densidad estimada, podemos generar un conjunto de puntos equidistantes entre sí y estimar la densidad del kernel en cada punto. Los puntos de prueba están dados por:

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

Ahora crearemos un objeto KernelDensity y usaremos el método fit() para encontrar el puntaje de cada muestra como se muestra en el código a continuación. El método KernelDensity() utiliza dos parámetros predeterminados, es decir, kernel=gaussian y bandwidth=1.

1
2
3
model = KernelDensity()
model.fit(x_train)
log_dens = model.score_samples(x_test)

La forma de la distribución se puede ver trazando el puntaje de densidad para cada punto, como se indica a continuación:

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

density score

Comprensión del parámetro de ancho de banda {#comprensión del parámetro de ancho de banda}

El ejemplo anterior no es una estimación muy impresionante de la función de densidad, atribuida principalmente a los parámetros predeterminados. Experimentemos con diferentes valores de ancho de banda para ver cómo afecta la estimación de la densidad.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
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()

different bandwidth estimation

Podemos ver claramente que aumentar el ancho de banda da como resultado una estimación más fluida. Los valores de ancho de banda muy pequeños dan como resultado curvas puntiagudas y nerviosas, mientras que los valores muy altos dan como resultado una curva suave muy generalizada que pierde detalles importantes. Es importante seleccionar un valor equilibrado para este parámetro.

Ajuste del parámetro de ancho de banda

La biblioteca scikit-learn permite ajustar el parámetro ancho de banda a través de la validación cruzada y devuelve el valor del parámetro que maximiza la probabilidad de registro de los datos. La función que podemos usar para lograr esto es GridSearchCV(), que requiere diferentes valores del parámetro bandwidth.

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

El mejor modelo se puede recuperar usando el campo best_estimator_ del objeto GridSearchCV.

Veamos la estimación óptima de la densidad del kernel utilizando el kernel gaussiano e imprimamos también el valor del ancho de banda:

1
2
3
4
5
6
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))

gaussian best estimator

1
optimal bandwidth: 0.15

Ahora, esta estimación de densidad parece modelar muy bien los datos. La primera mitad de la gráfica está de acuerdo con la distribución logarítmica normal y la segunda mitad de la gráfica modela bastante bien la distribución normal.

Kernels diferentes para la estimación de la densidad {# differentkernelsfordensityestimation}

scikit-learn permite la estimación de la densidad del kernel utilizando diferentes funciones del kernel:

  1. kernel ='coseno':      \(K(a;h) \propto \cos (\frac{\pi a}{2h}) \text { if } |a| < h \)
  2. núcleo = 'epanechnikov': \(K(a;h) \propto 1 - \frac{a^2}{h^2}\)
  3. `núcleo = ’exponencial’’: \(K(a;h) \propto \exp (-\frac{|a|}{h})\)
  4. núcleo = 'gaussiano':     \(K(a;h) \propto \exp(-\frac{a^2}{2h^2})\)
  5. núcleo = 'lineal':      \(K(a;h) \propto 1 - \frac{|a|}{h} \text { if } |a| < h \)
  6. kernel = 'tophat':      \(K(a;h) \propto 1 \text { if } |a| < h \)

Una forma sencilla de entender la forma en que funcionan estos núcleos es graficarlos. Esto significa construir un modelo usando una muestra de un solo valor, por ejemplo, 0. Luego, estime la densidad de todos los puntos alrededor de cero y trace la densidad a lo largo del eje y. El siguiente código muestra todo el proceso:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
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()

different kernels

Experimentando con diferentes núcleos {#experimentando con diferentes núcleos}

Experimentemos con diferentes núcleos y veamos cómo estiman la función de densidad de probabilidad para nuestros datos sintéticos.

Podemos usar GridSearchCV(), como antes, para encontrar el valor óptimo de ancho de banda. Sin embargo, para los núcleos coseno, lineal y tophat, GridSearchCV() puede generar una advertencia de tiempo de ejecución debido a que algunas puntuaciones dan como resultado valores -inf. Una forma posible de solucionar este problema es escribir una función de puntuación personalizada para GridSearchCV().

En el siguiente código, las puntuaciones -inf para los puntos de prueba se omiten en la función de puntuación personalizada my_scores() y se devuelve un valor medio. Este no es necesariamente el mejor esquema para manejar los valores de puntuación -inf y se puede adoptar alguna otra estrategia, dependiendo de los datos en cuestión.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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()

different kernels custom scoring

El modelo optimizado final

El ejemplo anterior muestra cómo diferentes núcleos estiman la densidad de diferentes maneras. Un paso final es configurar GridSearchCV() para que no solo descubra el ancho de banda óptimo, sino también el kernel óptimo para nuestros datos de ejemplo. Aquí está el código final que también traza la estimación de densidad final y sus parámetros ajustados en el título del gráfico:

1
2
3
4
5
6
7
8
9
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()

optimized kernel density model

Yendo más lejos: proyecto de extremo a extremo portátil

¿Tu naturaleza inquisitiva te hace querer ir más allá? Recomendamos consultar nuestro Proyecto guiado: ["Predicción práctica del precio de la vivienda: aprendizaje automático en Python"](https://wikihtp.com/courses/hands-on-house- precio-predicción-aprendizaje-máquina-en-python/#cta){target="_blank"}.

[](https://wikihtp.com/ cursos/predicción-de-precio-de-la-casa-práctica-aprendizaje-de-máquina-en-python/#cta)

En este proyecto guiado, aprenderá a crear potentes modelos tradicionales de aprendizaje automático, así como modelos de aprendizaje profundo, utilizar Ensemble Learning y capacitar a los meta-aprendices para predecir los precios de la vivienda a partir de una bolsa de modelos Scikit-Learn y Keras.

Con Keras, la API de aprendizaje profundo creada sobre Tensorflow, experimentaremos con arquitecturas, crearemos un conjunto de modelos apilados y entrenaremos una red neuronal meta-aprendizaje (modelo de nivel 1) para averiguar el precio de un casa.

El aprendizaje profundo es sorprendente, pero antes de recurrir a él, se recomienda intentar resolver el problema con técnicas más simples, como los algoritmos de aprendizaje superficial. Nuestro rendimiento de referencia se basará en un algoritmo de Regresión de bosque aleatorio. Además, exploraremos la creación de conjuntos de modelos a través de Scikit-Learn a través de técnicas como embalaje y votación.

Este es un proyecto integral y, como todos los proyectos de aprendizaje automático, comenzaremos con Análisis exploratorio de datos, seguido de Preprocesamiento de datos y, finalmente, Creación de modelos de aprendizaje superficial y Profundo para ajustarse a los datos que hemos explorado y limpiado previamente.

Conclusión

En este artículo se ha discutido la estimación de la densidad del kernel utilizando la biblioteca sklearn.neighbors de scikit-learn. Los ejemplos se dan para datos univariados, sin embargo, también se pueden aplicar a datos con múltiples dimensiones.

Si bien es una forma intuitiva y sencilla de estimar la densidad para distribuciones de fuentes desconocidas, un científico de datos debe usarla con precaución, ya que la maldición de la dimensionalidad puede ralentizarla considerablemente.