Estimación de Densidad de Kernel

En estadística, la estimación de la densidad de Kernel (KDE) es la aplicación del suavizado de kernel para la estimación de la densidad de probabilidad, es decir, un método no paramétrico para estimar la función de densidad de probabilidad de una variable aleatoria basada en kernels como pesos. KDE responde a un problema fundamental de suavizado de datos en el que se hacen inferencias sobre la población, basándose en una muestra finita de datos. En algunos campos, como el procesamiento de señales y la econometría, también se denomina método de la ventana de Parzen-Rosenblatt, en honor a Emanuel Parzen y Murray Rosenblatt, a quienes se suele atribuir su creación independiente en su forma actual.​​Una de las aplicaciones más conocidas de la estimación de la densidad del núcleo es la estimación de las densidades marginales condicionales de clase de los datos cuando se utiliza un clasificador Bayes ingenuo,​​ que puede mejorar su precisión de predicción.​

Estimación de la densidad del núcleo de 100 números aleatorios distribuidos normalmente utilizando diferentes anchos de banda de suavizado.

Definición

Sean (x1, x2, ..., xn) muestras independientes e idénticamente distribuidas extraídas de alguna distribución univariante con una densidad desconocida ƒ en cualquier punto dado x. Nos interesa estimar la forma de esta función ƒ. Su estimador de densidad kernel es

${\widehat {f}}_{h}(x)={\frac {1}{n}}\sum _{i=1}^{n}K_{h}(x-x_{i})={\frac {1}{nh}}\sum _{i=1}^{n}K{\Big (}{\frac {x-x_{i}}{h}}{\Big )},$

Donde K es el núcleo (una función no negativa) y h > 0 es un parámetro de suavizado denominado ancho de banda. Un kernel con el subíndice h se denomina kernel escalado y se define como Kh(x) = 1/h K(x/h). Intuitivamente, uno quiere elegir h tan pequeño como permitan los datos; sin embargo, siempre hay un compromiso entre el sesgo del estimador y su varianza. La elección del ancho de banda se analiza con más detalle a continuación.

Se suele utilizar una serie de funciones de kernel: uniforme, triangular, biponderada, triponderada, Epanechnikov, normal y otras. El kernel de Epanechnikov es óptimo en el sentido del error cuadrático medio,​ aunque la pérdida de eficiencia es pequeña para los kernels enumerados anteriormente.​ Debido a sus convenientes propiedades matemáticas, a menudo se utiliza el kernel normal, que significa K(x) = ϕ(x), donde ϕ es la función de densidad normal estándar.

La construcción de una estimación de densidad de núcleo encuentra interpretaciones en campos ajenos a la estimación de densidad.​ Por ejemplo, en termodinámica, equivale a la cantidad de calor generada cuando se colocan kernel de calor (la solución fundamental de la ecuación del calor) en cada ubicación de punto de datos xi. Se utilizan métodos similares para construir operadores discretos de Laplace en nubes de puntos para el aprendizaje de variedades (por ejemplo, el mapa de difusión).

Ejemplo

Las estimaciones de densidad de kernel están estrechamente relacionadas con los histogramas, pero pueden dotarse de propiedades como la suavidad o la continuidad mediante el uso de un kernel adecuado. El siguiente diagrama basado en estos 6 puntos de datos ilustra esta relación:

Muestra	1	2	3	4	5	6
Valor	-2.1	-1.3	-0.4	1.9	5.1	6.2

Para el histograma, primero se divide el eje horizontal en subintervalos o "bins" que cubren el intervalo de los datos: En este caso, seis intervalos de 2. Cada vez que un dato cae dentro de este intervalo, se coloca en él una caja de 1/12 de altura. Si hay más de un dato dentro del mismo intervalo, las cajas se apilan unas encima de otras.

Para la estimación de la densidad del núcleo, se colocan núcleos normales con una desviación estándar de 1,5 (indicados por las líneas discontinuas rojas) en cada uno de los puntos de datos xi. Los núcleos se suman para obtener la estimación de la densidad del núcleo (curva azul continua). La suavidad de la estimación de la densidad del kernel (en comparación con la discreción del histograma) ilustra cómo las estimaciones de la densidad del kernel convergen más rápidamente a la verdadera densidad subyacente para variables aleatorias continuas.​

Comparación del histograma (izquierda) y la estimación de la densidad del kernel (derecha) construidos con los mismos datos. Los seis núcleos individuales son las curvas rojas discontinuas, y la estimación de la densidad del núcleo, las curvas azules. Los puntos de datos son el rug plot (diagrama de alfombra) en el eje horizontal.

Selección de ancho de banda

Estimación de la densidad del núcleo (KDE) con diferentes anchos de banda de una muestra aleatoria de 100 puntos de una distribución normal estándar. Gris: densidad real (normal estándar). Rojo: KDE con h=0,05. Negro: KDE con h=0,337. Verde: KDE con h=2.

El ancho de banda del kernel es un parámetro libre que influye mucho en la estimación resultante. Para ilustrar su efecto, tomamos una muestra aleatoria simulada de la distribución normal estándar (representada en los picos azules del rug plot en el eje horizontal). La curva gris es la densidad real (una densidad normal con media 0 y varianza 1). En comparación, la curva roja está poco suavizada, ya que contiene demasiados datos espurios debido a la utilización de un ancho de banda h = 0,05, que es demasiado pequeño. La curva verde está demasiado suavizada, ya que el ancho de banda h = 2 oculta gran parte de la estructura subyacente. La curva negra con un ancho de banda de h = 0,337 se considera óptimamente suavizada, ya que su estimación de la densidad se aproxima a la densidad real. Una situación extrema se da en el límite $h\to 0$ (sin suavizado), donde la estimación es una suma de n funciones delta centradas en las coordenadas de las muestras analizadas. En el otro límite extremo $h\to \infty$ la estimación conserva la forma del núcleo utilizado, centrada en la media de las muestras (completamente lisa).

El criterio de optimalidad más utilizado para seleccionar este parámetro es la función de riesgo L2 esperado, también denominada error cuadrático medio integrado:

$\operatorname {MISE} (h)=\operatorname {E} \!\left[\,\int ({\hat {f}}_{h}(x)-f(x))^{2}\,dx\right]$

Bajo supuestos débiles sobre ƒ y K, (ƒ es la función de densidad real, generalmente desconocida),​​

$\operatorname {MISE} (h)=\operatorname {AMISE} (h)+{\mathcal {o}}((nh)^{-1}+h^{4})$

donde o es la notación o pequeña, y n el tamaño de la muestra (como arriba). El AMISE es el MISE asintótico, es decir, los dos términos principales,

$\operatorname {AMISE} (h)={\frac {R(K)}{nh}}+{\frac {1}{4}}m_{2}(K)^{2}h^{4}R(f'')$

Donde $R(g)=\int g(x)^{2}\,dx$ para una función g, $m_{2}(K)=\int x^{2}K(x)\,dx$ y $f''$ es la segunda derivada de $f$ y $K$ es el kernel. El mínimo de este AMISE es la solución de esta ecuación diferencial.

${\frac {\partial }{\partial h}}\operatorname {AMISE} (h)=-{\frac {R(K)}{nh^{2}}}+m_{2}(K)^{2}h^{3}R(f'')=0$

o

$h_{\operatorname {AMISE} }={\frac {R(K)^{1/5}}{m_{2}(K)^{2/5}R(f'')^{1/5}}}n^{-1/5}=Cn^{-1/5}$

Ni las fórmulas AMISE ni las hAMISE pueden utilizarse directamente, ya que implican la función de densidad desconocida $f$ o su segunda derivada $f''$. Para superar esa dificultad, se han desarrollado diversos métodos automáticos basados en datos para seleccionar el ancho de banda. Se han realizado varios estudios de revisión para comparar su eficacia,​​​​​​​ con el consenso general de que los selectores de complemento​​​ y los selectores de validación cruzada​​​ son los más útiles en una amplia gama de conjuntos de datos.

Sustituyendo cualquier ancho de banda h que tenga el mismo orden asintótico n^-1/5 que h_AMISE en el AMISE se obtiene que AMISE(h) = O(n^-4/5), donde O es la notación big o. Puede demostrarse que, bajo supuestos débiles, no puede existir un estimador no paramétrico que converja a un ritmo más rápido que el estimador kernel.​ Obsérvese que el ritmo n^-4/5 es más lento que el típico ritmo de convergencia n^-1 de los métodos paramétricos.

Si el ancho de banda no se mantiene fijo, sino que varía en función de la ubicación de la estimación (estimador de globo) o de las muestras (estimador puntual), se obtiene un método especialmente potente denominado estimación de densidad kernel adaptativa o de ancho de banda variable.

La selección del ancho de banda para la estimación de la densidad del kernel de distribuciones con colas gruesas es relativamente difícil.​

Un estimador de ancho de banda basado en reglas empíricas

Si se utilizan funciones de base gaussiana para aproximar datos univariantes, y la densidad subyacente que se estima es gaussiana, la elección óptima para h (es decir, el ancho de banda que minimiza el error cuadrático medio integrado) es:​

$h=\left({\frac {4{\hat {\sigma }}^{5}}{3n}}\right)^{\frac {1}{5}}\approx 1.06\,{\hat {\sigma }}\,n^{-1/5},$

Un valor $h$ se considera más robusto cuando mejora el ajuste para distribuciones de cola larga y sesgadas o para distribuciones de mezcla bimodales. Esto suele hacerse empíricamente sustituyendo la desviación típica ${\hat {\sigma }}$ por el parámetro $A$ a continuación:

$A=\min \left({\hat {\sigma }},{\frac {IQR}{1.34}}\right)$ donde IQR es el rango intercuartílico.

Comparación entre la regla empírica y el ancho de banda de resolución de ecuaciones.

Otra modificación que mejorará el modelo es reducir el factor de 1,06 a 0,9. Entonces la fórmula final sería:

$h=0.9\,\min \left({\hat {\sigma }},{\frac {IQR}{1.34}}\right)\,n^{-{\frac {1}{5}}}$

Esta aproximación se denomina aproximación de distribución normal, aproximación gaussiana o regla empírica de Silverman.​ Aunque esta regla empírica es fácil de calcular, debe utilizarse con precaución, ya que puede dar lugar a estimaciones muy imprecisas cuando la densidad no se aproxima a la normalidad. Por ejemplo, al estimar el modelo de mezcla gaussiana bimodal:

${\frac {1}{2{\sqrt {2\pi }}}}e^{-{\frac {1}{2}}(x-10)^{2}}+{\frac {1}{2{\sqrt {2\pi }}}}e^{-{\frac {1}{2}}(x+10)^{2}}$

a partir de una muestra de 200 puntos, la figura de la derecha muestra la densidad real y dos estimaciones de densidad kernel: una que utiliza el ancho de banda de la regla empírica y otra que utiliza el ancho de banda de la solución de la ecuación.​​ La estimación basada en el ancho de banda de la regla empírica está muy suavizada.

Relación con el estimador de densidad de la función característica

Dada la muestra (x1, x2, ..., xn), es natural estimar la función característica φ(t) = E[e^itX] como:

${\widehat {\varphi }}(t)={\frac {1}{n}}\sum _{j=1}^{n}e^{itx_{j}}$

Conociendo la función característica, es posible hallar la correspondiente función de densidad de probabilidad mediante la fórmula de la transformada de Fourier. Una dificultad de la aplicación de esta fórmula de inversión es que conduce a una integral divergente, ya que la estimación $\scriptstyle {\widehat {\varphi }}(t)$ no es fiable para t's grandes. Para evitar este problema, el estimador $\scriptstyle {\widehat {\varphi }}(t)$ se multiplica por una función de amortiguación ψ_h(t) = ψ(ht) , que es igual a 1 en el origen y luego cae a 0 en el infinito. El "parámetro de ancho de banda" h controla la rapidez con la que intentamos amortiguar la función $\scriptstyle {\widehat {\varphi }}(t)$. En particular, cuando h es pequeño, entonces ψ_h(t) será aproximadamente uno para un amplio rango de t's, lo que significa que $\scriptstyle {\widehat {\varphi }}(t)$ permanece prácticamente inalterado en la región más importante de t's.

La elección más común para la función ψ es la función uniforme ψ(t) = 1{−1 ≤ t ≤ 1}, lo que significa truncar el intervalo de integración en la fórmula de inversión a [−1/h, 1/h] o la función gaussiana ψ (t) = e^πt2. Una vez elegida la función ψ, se puede aplicar la fórmula de inversión, y el estimador de densidad será:

${\begin{aligned}{\widehat {f}}(x)&={\frac {1}{2\pi }}\int _{-\infty }^{+\infty }{\widehat {\varphi }}(t)\psi _{h}(t)e^{-itx}\,dt={\frac {1}{2\pi }}\int _{-\infty }^{+\infty }{\frac {1}{n}}\sum _{j=1}^{n}e^{it(x_{j}-x)}\psi (ht)\,dt\\[5pt]&={\frac {1}{nh}}\sum _{j=1}^{n}{\frac {1}{2\pi }}\int _{-\infty }^{+\infty }e^{-i(ht){\frac {x-x_{j}}{h}}}\psi (ht)\,d(ht)={\frac {1}{nh}}\sum _{j=1}^{n}K{\Big (}{\frac {x-x_{j}}{h}}{\Big )},\end{aligned}}$

donde K es la transformada de Fourier de la función de amortiguación ψ. Así pues, el estimador de densidad kernel coincide con el estimador de densidad de función característica.

Características geométricas y topológicas

Podemos ampliar la definición del modo (global) a un sentido local y definir los modos locales:

$M=\{x:g(x)=0,\lambda _{1}(x)<0\}$

Concretamente, $M$ es la colección de puntos para los que la función de densidad se maximiza localmente. Un estimador natural de $M$ es un plug-in de KDE,​​ donde $g(x)$ y $\lambda _{1}(x)$ son la versión KDE de $g(x)$ y $\lambda _{1}(x)$. Bajo suposiciones suaves, $M_{c}$ es un estimador coherente de $M$. Obsérvese que se puede utilizar el algoritmo de desplazamiento de la media​ ​​para calcular el estimador $M_{c}$ numéricamente.

Aplicación estadística

Una lista no exhaustiva de implementaciones de software de estimadores de densidad kernel incluye:

• En la versión 4.4 de Analytica, la opción de suavizado para resultados PDF utiliza KDE, y a partir de las expresiones está disponible a través de la función Pdf incorporada.
• En C/C++, FIGTree es una biblioteca que puede utilizarse para calcular estimaciones de densidad de núcleos utilizando núcleos normales. Interfaz MATLAB disponible.
• En C++, libagf es una biblioteca para la estimación de la densidad del kernel variable.
• En C++, mlpack es una biblioteca que puede calcular KDE utilizando muchos núcleos diferentes. Permite establecer una tolerancia de error para un cálculo más rápido. Dispone de interfaces para Python y R.
• en C# y F#, Math.NET Numerics es una biblioteca de código abierto para cálculo numérico que incluye la estimación de la densidad del kernel.
• En CrimeStat, la estimación de la densidad del núcleo se implementa utilizando cinco funciones de núcleo diferentes: normal, uniforme, cuártica, exponencial negativa y triangular. Se dispone de rutinas de estimación de densidad de núcleo simple y doble. La estimación de la densidad del núcleo también se utiliza en la interpolación de una rutina Head Bang, en la estimación de una función de densidad bidimensional Journey-to-crime y en la estimación de una estimación Bayesiana Journey-to-crime tridimensional.
• En ELKI, las funciones de densidad del núcleo se encuentran en el paquete de.lmu.ifi.dbs.elki.math.statistics.kernelfunctions
• En los productos ESRI, la cartografía de densidad de kernel se gestiona fuera de la caja de herramientas Spatial Analyst y utiliza el kernel Quartic(biweight).
• En Excel, la Royal Society of Chemistry ha creado un complemento para la estimación de la densidad del núcleo basado en su Informe Técnico 4 del Comité de Métodos Analíticos.
• En gnuplot, la estimación de la densidad kernel se implementa mediante la opción smooth kdensity, el fichero de datos puede contener un peso y un ancho de banda para cada punto, o el ancho de banda puede establecerse automáticamente​ según la "regla empírica de Silverman" (véase más arriba).
• En Haskell, la densidad del núcleo se implementa en el paquete statistics.
• En IGOR Pro, la estimación de la densidad del kernel se implementa mediante la operación StatsKDE (añadida en Igor Pro 7.00). El ancho de banda puede ser especificado por el usuario o estimado mediante Silverman, Scott o Bowmann y Azzalini. Los tipos de núcleo son: Epanechnikov, Bi-peso, Tri-peso, Triangular, Gaussiano y Rectangular.
• En Java, el paquete de aprendizaje automático Weka proporciona weka.estimators.KernelEstimator, entre otros.
• En JavaScript, el paquete de visualización D3.js ofrece un paquete KDE en su paquete science.stats.
• En el JMP, la plataforma Graph Builder utiliza la estimación de densidad kernel para proporcionar gráficos de contorno y regiones de alta densidad (HDR) para densidades bivariantes, y gráficos de violín y HDR para densidades univariantes. Los controles deslizantes permiten al usuario variar el ancho de banda. Las plataformas Ajustar Y por X y Distribución también proporcionan estimaciones de densidad de kernel bivariantes y univariantes, respectivamente.
• En Julia, la estimación de la densidad del kernel se implementa en el paquete KernelDensity.jl.
• En KNIME, las distribuciones 1D y 2D Kernel Density pueden ser generadas y trazadas usando nodos de la contribución de la comunidad Vernalis, por ejemplo 1D Kernel Density Plot, entre otros. La implementación subyacente está escrita en Java.
• En MATLAB, la estimación de la densidad del kernel se implementa a través de la función ksdensity (Statistics Toolbox). A partir de la versión 2018a de MATLAB, se pueden especificar tanto el ancho de banda como el suavizador del kernel, incluyendo otras opciones como especificar el rango de la densidad del kernel.​ Alternativamente, en el MATLAB Central File Exchange está disponible un paquete de software gratuito de MATLAB que implementa un método de selección automática del ancho de banda​ para:
  • Datos unidimensionales
  • Datos bidimensionales
  • datos de n dimensiones

En estas páginas está disponible una caja de herramientas MATLAB gratuita con la implementación de la regresión kernel, la estimación de la densidad kernel, la estimación kernel de la función de riesgo y muchas otras (esta caja de herramientas forma parte del libro)​.

• En Mathematica, la estimación numérica de la densidad del kernel se implementa mediante la función SmoothKernelDistribution​ y la estimación simbólica se implementa mediante la función KernelMixtureDistribution​, las cuales proporcionan anchos de banda en función de los datos.
• En Minitab, la Royal Society of Chemistry ha creado una macro para ejecutar la estimación de densidad de kernel basada en su Informe técnico 4 del Comité de Métodos Analíticos​.
• En la Biblioteca NAG, la estimación de la densidad del kernel se implementa mediante la rutina g10ba (disponible tanto en la versión Fortran​ como en la versión C​ de la Biblioteca).
• En Nuklei, los métodos de densidad de núcleo C++ se centran en los datos del grupo euclidiano especial $SE(3)$
• En Octave, la estimación de la densidad kernel se implementa mediante la opción kernel_density (paquete econometrics).
• En Origin, el gráfico de densidad de kernel 2D puede realizarse desde su interfaz de usuario, y dos funciones, Ksdensity para 1D y Ks2density para 2D pueden utilizarse desde su código LabTalk, Python o C.
• En Perl, puede encontrarse una implementación en el módulo Statistics-KernelEstimation.
• En PHP, se puede encontrar una implementación en la biblioteca MathPHP
• En Python, existen muchas implementaciones: pyqt_fit.kde Module en el paquete PyQt-Fit, SciPy (scipy.stats.gaussian_kde), Statsmodels (KDEUnivariate y KDEMultivariate) y scikit-learn (KernelDensity) (véase la comparación)​. KDEpy soporta datos ponderados y su implementación FFT es órdenes de magnitud más rápida que las otras implementaciones. La biblioteca pandas de uso común​ ofrece soporte para el trazado kde a través del método plot (df.plot(kind='kde'))​). El paquete getdist para muestras MCMC ponderadas y correlacionadas admite ancho de banda optimizado, corrección de límites y métodos de orden superior para distribuciones 1D y 2D. Un paquete utilizado recientemente para la estimación de la densidad del núcleo es seaborn (import seaborn as sns , sns.kdeplot())​ También existe una implementación de KDE en la GPU.​
• En R, se implementa a través de la density en la distribución base, y la función bw.nrd0 se utiliza en el paquete stats, esta función utiliza la fórmula optimizada en el libro de Silverman. bkde en la biblioteca KernSmooth, ParetoDensityEstimation en la biblioteca DataVisualizations (para la estimación de la densidad de la distribución pareto), kde en la biblioteca ks, kde y dbckden en la biblioteca evmix (esta última para la estimación de la densidad del kernel con corrección de límites para soporte acotado), npudens en la biblioteca np (datos numéricos y categóricos), sm.density en la biblioteca sm. Para una implementación de la función kde.R, que no requiere la instalación de ningún paquete o biblioteca, véase kde.R. La biblioteca btb, dedicada al análisis urbano, implementa la estimación de densidad kernel mediante kernel_smoothing.
• En SAS, proc kde puede utilizarse para estimar densidades kernel univariantes y bivariantes.
• En Apache Spark, la clase KernelDensity()()​
• En Stata, se implementa a través de kdensity;​ por ejemplo histogram x, kdensity. Alternativamente, existe un módulo gratuito de Stata, KDENS​, que permite al usuario estimar funciones de densidad 1D o 2D.
• En Swift, se implementa a través de SwiftStats.KernelDensityEstimation en la biblioteca estadística de código abierto.

Véase también

• Kernel (estadísticas)
• Desplazamiento cuadrático medio
• Escala (espacial) : Las tripletas {(x, h, KDE con ancho de banda h evaluado en x: todo x, h > 0} forman una representación de espacio de escala de los datos.
• Cola (informática)
• Mean shift

Lectura adicional

• Härdle, Müller, Sperlich, Werwatz, Nonparametric and Semiparametric Methods, Springer-Verlag Berlin Heidelberg 2004, pp. 39–83

Enlaces externos

• Introduction to kernel density estimation A short tutorial which motivates kernel density estimators as an improvement over histograms.
• Kernel Bandwidth Optimization A free online tool that generates an optimized kernel density estimate.
• Free Online Software (Calculator) computes the Kernel Density Estimation for a data series according to the following Kernels: Gaussian, Epanechnikov, Rectangular, Triangular, Biweight, Cosine, and Optcosine.
• Kernel Density Estimation Applet An online interactive example of kernel density estimation. Requires .NET 3.0 or later.