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

Introducción

Este artículo es una introducción a la estimación de la densidad del kernel utilizando la biblioteca de Machine Learning de Python scikit-learn.

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

Dada una muestra de observaciones independientes, idénticamente distribuidas (iid) ((x_1, x_2, ldots, x_n) ) de una variable aleatoria de una distribución de fuente desconocida, la estimación de la densidad del núcleo 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 del kernel y (h ) es el parámetro de suavizado, también llamado ancho de banda. Más adelante en este artículo se analizan varios núcleos, pero solo para comprender las matemáticas, echemos un vistazo a un ejemplo simple.

Ejemplo de cálculo

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

x j =

[

−
2

Te puede interesar:Python: compruebe si el archivo o directorio está vacío

−
1

0

1

2

]

| 0 – x j | =

[

2

1

0

Te puede interesar:Uso de la biblioteca Plotly para la visualización interactiva de datos en Python

1

2

]

| 0 – x j h | =

[

0.2

0.1

0

0.1

0.2

Te puede interesar:Biblioteca Seaborn para la visualización de datos en Python: Parte 1

]

K (| 0 – x j h |) =

[

0.8

0.9

1

0.9

0.8

]

Reemplaza 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
$$

Te puede interesar:Introducción al estilo de codificación Python

Estimación de la densidad del kernel con Python

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

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 kernel, se generan datos sintéticos a partir de dos tipos diferentes de distribuciones. Uno es una distribución log-normal asimétrica y el otro es una distribución gaussiana. La siguiente función devuelve 2000 puntos de datos:

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 diagrama de dispersión de estos puntos a lo largo del eje y o podemos generar un histograma de estos puntos.

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()

Uso de KernelDensity 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 núcleo en cada punto. Los puntos de prueba vienen dados por:

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

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

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

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

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

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 densidad.

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()

Podemos ver claramente que aumentar el ancho de banda da como resultado una estimación más suave. Los valores de ancho de banda muy pequeños dan como resultado curvas puntiagudas y temblorosas, 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

los scikit-learn biblioteca permite la sintonización de la bandwidth parámetro mediante validación cruzada y devuelve el valor del parámetro que maximiza la probabilidad logarítmica de los datos. La función que podemos usar para lograr esto es GridSearchCV(), que requiere diferentes valores de la bandwidth parámetro.

Te puede interesar:Tutorial de Python para principiantes absolutos

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 utilizando el best_estimator_ campo de la GridSearchCV objeto.

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

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

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

Diferentes núcleos para la estimación de densidad

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

• kernel="cosine": (K (a; h) propto cos ( frac { pi a} {2h}) text {if} | a |
• kernel="epanechnikov": (K (a; h) propto 1 – frac {a ^ 2} {h ^ 2} )
• kernel="exponential": (K (a; h) propto exp (- frac {| a |} {h}) )
• kernel="gaussian": (K (a; h) propto exp (- frac {a ^ 2} {2h ^ 2}) )
• kernel="linear": (K (a; h) propto 1 – frac {| a |} {h} text {si} | a |
• kernel="tophat": (K (a; h) propto 1 text {si} | a |Una forma sencilla de comprender la forma en que funcionan estos núcleos es trazarlos. 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 grafique la densidad a lo largo del eje y. El siguiente código muestra todo el proceso:

  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()
  

  Experimentar 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.

  Nosotros podemos usar GridSearchCV(), como antes, para encontrar el óptimo bandwidth valor. Sin embargo, para cosine, lineary tophat granos GridSearchCV() puede dar una advertencia de tiempo de ejecución debido a algunas puntuaciones que resultan en -inf valores. Una forma posible de abordar este problema es escribir una función de puntuación personalizada para GridSearchCV().

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

  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()
  

  El modelo optimizado final

  El ejemplo anterior muestra cómo diferentes núcleos estiman la densidad de diferentes maneras. Un último paso 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 de la trama:

  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()
  

  Conclusión

  Estimación de la densidad de kernel usando scikit-learnbiblioteca de sklearn.neighbors se ha discutido en este artículo. Los ejemplos se dan para datos univariados, sin embargo, también se pueden aplicar a datos con múltiples dimensiones.

  Te puede interesar:Algoritmo de bosque aleatorio con Python y Scikit-Learn

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