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


    1

    0

    1

    2

    ]

    | 0 – x j | =

    [

    2

    1

    0

    1

    2

    ]

    | 0 – x j h | =

    [

    0.2

    0.1

    0

    0.1

    0.2

    ]

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

    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.

    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.

      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.

       

    Etiquetas:

    Deja una respuesta

    Tu direcci贸n de correo electr贸nico no ser谩 publicada. Los campos obligatorios est谩n marcados con *