Estimación de densidad de núcleo multivariante

La estimación de densidad kernel es una técnica no paramétrica para la estimación de densidad , es decir, la estimación de funciones de densidad de probabilidad , que es una de las cuestiones fundamentales en estadística . Puede considerarse una generalización de la estimación de densidad de histogramas con propiedades estadísticas mejoradas. Además de los histogramas, otros tipos de estimadores de densidad incluyen los paramétricos , splines , wavelets y series de Fourier . Los estimadores de densidad kernel se introdujeron por primera vez en la literatura científica para datos univariados en las décadas de 1950 y 1960 ^{[ 1 ] [ 2 ]} y posteriormente se han adoptado ampliamente. Pronto se reconoció que los estimadores análogos para datos multivariados serían una importante adición a la estadística multivariada . Basándose en la investigación llevada a cabo en las décadas de 1990 y 2000, la estimación de densidad kernel multivariada ha alcanzado un nivel de madurez comparable al de sus contrapartes univariadas. ^{[ 3 ] [ 4 ] [ 5 ]}

Utilizamos un conjunto de datos sintéticos bivariados de 50 puntos para ilustrar la construcción de histogramas. Esto requiere la elección de un punto de anclaje (la esquina inferior izquierda de la cuadrícula del histograma). Para el histograma de la izquierda, elegimos (−1.5,  −1.5); para el de la derecha, desplazamos el punto de anclaje 0.125 en ambas direcciones a (−1.625,  −1.625). Ambos histogramas tienen un ancho de intervalo de 0.5, por lo que cualquier diferencia se debe únicamente al cambio en el punto de anclaje. La codificación por colores indica la cantidad de puntos de datos que caen en un intervalo: 0=blanco, 1=amarillo pálido, 2=amarillo brillante, 3=naranja, 4=rojo. El histograma de la izquierda parece indicar que la mitad superior tiene una densidad mayor que la mitad inferior, mientras que ocurre lo contrario con el histograma de la derecha, lo que confirma que los histogramas son muy sensibles a la ubicación del punto de anclaje. ^{[ 6 ]}

Una posible solución a este problema de ubicación de puntos de anclaje es eliminar por completo la cuadrícula de agrupamiento del histograma. En la figura de la izquierda, un núcleo (representado por las líneas grises) está centrado en cada uno de los 50 puntos de datos anteriores. El resultado de la suma de estos núcleos se muestra en la figura de la derecha, que es una estimación de la densidad del núcleo. La diferencia más notable entre las estimaciones de densidad del núcleo y los histogramas es que las primeras son más fáciles de interpretar, ya que no contienen artificios inducidos por una cuadrícula de agrupamiento. Los contornos de color corresponden a la región más pequeña que contiene la masa de probabilidad respectiva: rojo = 25%, naranja + rojo = 50%, amarillo + naranja + rojo = 75%, lo que indica que una única región central contiene la mayor densidad.

El objetivo de la estimación de densidad es tomar una muestra finita de datos y realizar inferencias sobre la función de densidad de probabilidad subyacente en todo el espacio, incluso donde no se observan datos. En la estimación de densidad de kernel, la contribución de cada punto de datos se suaviza desde un único punto hacia una región del espacio que lo rodea. La agregación de las contribuciones suavizadas individualmente proporciona una visión general de la estructura de los datos y su función de densidad. En los detalles que se presentan a continuación, demostramos que este enfoque conduce a una estimación razonable de la función de densidad subyacente.

La figura anterior es una representación gráfica de la estimación de densidad del núcleo, que ahora definimos de manera exacta. Sean x ₁ , x ₂ , ..., x _n una muestra de vectores aleatorios de d variables extraídos de una distribución común descrita por la función de densidad ƒ . La estimación de densidad del núcleo se define como

${\displaystyle {\hat {f}}_{\mathbf {H} }(\mathbf {x} )={\frac {1}{n}}\sum _{i=1}^{n}K_{\mathbf {H} }(\mathbf {x} -\mathbf {x} _{i})}$

dónde

• x = ( x ₁ , x ₂ , …, x _d ) ^T , x _i = ( x _{i 1} , x _{i 2} , …, x _id ) ^T , i = 1, 2, …, n son d- vectores;
• H es la matriz de ancho de banda (o suavizado) d×d que es simétrica y definida positiva ;
• K es la función núcleo , que es una densidad multivariada simétrica;
• ${\displaystyle K_{\mathbf {H} }(\mathbf {x} )=|\mathbf {H} |^{-1/2}K(\mathbf {H} ^{-1/2}\mathbf {x} )}$.

La elección de la función núcleo K no es crucial para la precisión de los estimadores de densidad de núcleo, por lo que utilizamos el núcleo normal multivariado estándar en todo momento:${\textstyle K_{\mathbf {H} }(\mathbf {x} )={(2\pi )^{-d/2}}\mathbf {|H|} ^{-1/2}e^{-{\frac {1}{2}}\mathbf {x^{T}} \mathbf {H^{-1}} \mathbf {x} }}$donde H desempeña el papel de la matriz de covarianza . Por otro lado, la elección de la matriz de ancho de banda H es el factor más importante que afecta su precisión, ya que controla la cantidad y la orientación del suavizado inducido. ^{[ 3 ] : 36–39} Que la matriz de ancho de banda también induzca una orientación es una diferencia básica entre la estimación de densidad de kernel multivariante y su análogo univariante, ya que la orientación no está definida para kernels 1D. Esto lleva a la elección de la parametrización de esta matriz de ancho de banda. Las tres clases principales de parametrización (en orden creciente de complejidad) son S , la clase de escalares positivos multiplicados por la matriz identidad; D , matrices diagonales con entradas positivas en la diagonal principal ; y F , matrices simétricas definidas positivas. Los kernels de la clase S tienen la misma cantidad de suavizado aplicada en todas las direcciones de coordenadas, los kernels D permiten diferentes cantidades de suavizado en cada una de las coordenadas, y los kernels F permiten cantidades y orientación arbitrarias del suavizado. Históricamente, los núcleos S y D son los más extendidos debido a razones computacionales, pero las investigaciones indican que se pueden obtener mejoras importantes en la precisión utilizando los núcleos de clase F más generales . ^{[ 7 ] [ 8 ]}

El criterio de optimalidad más utilizado para seleccionar una matriz de ancho de banda es el MISE o error cuadrático medio integrado.

${\displaystyle \operatorname {MISE} (\mathbf {H} )=\operatorname {E} \!\left[\,\int ({\hat {f}}_{\mathbf {H} }(\mathbf {x} )-f(\mathbf {x} ))^{2}\,d\mathbf {x} \;\right].}$

En general, esta función no posee una expresión analítica , por lo que se suele utilizar su aproximación asintótica (AMISE) como indicador.

${\displaystyle \operatorname {AMISE} (\mathbf {H} )=n^{-1}|\mathbf {H} |^{-1/2}R(K)+{\tfrac {1}{4}}m_{2}(K)^{2}(\operatorname {vec} ^{T}\mathbf {H} )\mathbf {\Psi } _{4}(\operatorname {vec} \mathbf {H} )}$

dónde

• ${\displaystyle R(K)=\int K(\mathbf {x} )^{2}\,d\mathbf {x} }$, con R ( K ) = (4 π ) ^{−d /2} cuando K es un núcleo normal
• ${\displaystyle \int \mathbf {x} \mathbf {x} ^{T}K(\mathbf {x} )\,d\mathbf {x} =m_{2}(K)\mathbf {I} _{d}}$,

donde I _d es la matriz identidad d × d , con m ₂ = 1 para el núcleo normal.

• D ² ƒ es la matriz hessiana d × d de las derivadas parciales de segundo orden de ƒ
• ${\displaystyle \mathbf {\Psi } _{4}=\int (\operatorname {vec} \,\operatorname {D} ^{2}f(\mathbf {x} ))(\operatorname {vec} ^{T}\operatorname {D} ^{2}f(\mathbf {x} ))\,d\mathbf {x} }$es una matriz d ² × d ² de derivadas parciales de cuarto orden integradas de ƒ
• vec es el operador vectorial que apila las columnas de una matriz en un único vector, por ejemplo${\displaystyle \operatorname {vec} {\begin{bmatrix}a&c\\b&d\end{bmatrix}}={\begin{bmatrix}a&b&c&d\end{bmatrix}}^{T}.}$

La calidad de la aproximación AMISE a MISE ^{[ 3 ] : 97} viene dada por

${\displaystyle \operatorname {MISE} (\mathbf {H} )=\operatorname {AMISE} (\mathbf {H} )+o(n^{-1}|\mathbf {H} |^{-1/2}+\operatorname {tr} \,\mathbf {H} ^{2})}$

donde o indica la notación habitual de o pequeña . Heurísticamente, esta afirmación implica que el AMISE es una "buena" aproximación del MISE cuando el tamaño de la muestra n → ∞.

Se puede demostrar que cualquier selector de ancho de banda razonable H tiene H = O ( n ^{−2/( d +4)} ) donde la notación O grande se aplica elemento a elemento. Sustituyendo esto en la fórmula MISE se obtiene que el MISE óptimo es O ( n ^{−4/( d +4)} ). ^{[ 3 ] : 99–100} Por lo tanto, cuando n → ∞, el MISE → 0, es decir, la estimación de densidad del núcleo converge en media cuadrática y, por lo tanto, también en probabilidad a la densidad verdadera f . Estos modos de convergencia confirman la afirmación en la sección de motivación de que los métodos de núcleo conducen a estimadores de densidad razonables. Un selector de ancho de banda óptimo ideal es

${\displaystyle \mathbf {H} _{\operatorname {AMISE} }=\operatorname {argmin} _{\mathbf {H} \in F}\,\operatorname {AMISE} (\mathbf {H} ).}$

Dado que este selector ideal contiene la función de densidad desconocida ƒ , no puede utilizarse directamente. Las diversas variantes de selectores de ancho de banda basados ​​en datos surgen de los diferentes estimadores del AMISE. Nos centramos en dos clases de selectores que han demostrado ser los más aplicables en la práctica: la validación cruzada suavizada y los selectores de sustitución.

En la sección de selección de ancho de banda óptimo, presentamos el MISE. Su construcción se basa en el valor esperado y la varianza del estimador de densidad ^{[ 3 ] : 97}

${\displaystyle \operatorname {E} {\hat {f}}(\mathbf {x}$ ;\mathbf {H} )=K_{\mathbf {H} }*f(\mathbf {x} )=f(\mathbf {x} )+{\frac {1}{2}}m_{2}(K)\operatorname {tr} (\mathbf {H} \operatorname {D} ^{2}f(\mathbf {x} ))+o(\operatorname {tr} \,\mathbf {H} )}

donde * es el operador de convolución entre dos funciones, y

${\displaystyle \operatorname {Var} {\hat {f}}(\mathbf {x}$ ;\mathbf {H} )=n^{-1}|\mathbf {H} |^{-1/2}R(K)f(\mathbf {x} )+o(n^{-1}|\mathbf {H} |^{-1/2}).}

Para que estas dos expresiones estén bien definidas, requerimos que todos los elementos de H tiendan a 0 y que n ⁻¹ | H | ^−1/2 tienda a 0 cuando n tienda a infinito. Suponiendo estas dos condiciones, vemos que el valor esperado tiende a la densidad verdadera f , es decir, el estimador de densidad de núcleo es asintóticamente insesgado ; y que la varianza tiende a cero. Usando la descomposición estándar del valor cuadrático medio

${\displaystyle \operatorname {MSE} \,{\hat {f}}(\mathbf {x}$ ;\mathbf {H} )=\operatorname {Var} {\hat {f}}(\mathbf {x}  ;\mathbf {H} )+[\operatorname {E} {\hat {f}}(\mathbf {x}  ;\mathbf {H} )-f(\mathbf {x} )]^{2}}

Tenemos que el MSE tiende a 0, lo que implica que el estimador de densidad del núcleo es consistente (en media cuadrática) y, por lo tanto, converge en probabilidad a la densidad verdadera f . La tasa de convergencia del MSE a 0 es necesariamente la misma que la tasa del MISE mencionada anteriormente O ( n ^−4/(d+4) ), por lo tanto, la tasa de convergencia del estimador de densidad a f es O _p (n ^{−2/( d +4)} ) donde O _p denota el orden en probabilidad . Esto establece la convergencia puntual. La convergencia funcional se establece de manera similar al considerar el comportamiento del MISE y observar que, bajo suficiente regularidad, la integración no afecta las tasas de convergencia.

Para los selectores de ancho de banda basados ​​en datos considerados, el objetivo es la matriz de ancho de banda AMISE. Decimos que un selector basado en datos converge al selector AMISE a una tasa relativa O _p ( n ^{− α} ), α > 0 si

${\displaystyle \operatorname {vec} ({\hat {\mathbf {H} }}-\mathbf {H} _{\operatorname {AMISE} })=O(n^{-2\alpha })\operatorname {vec} \mathbf {H} _{\operatorname {AMISE} }.}$

Se ha establecido que los selectores de validación cruzada de inserción y suavizada (dado un único ancho de banda piloto G ) convergen a una tasa relativa de O _p ( n ^{−2/( d +6)} ) ^{[ 10 ] [ 12 ]} es decir, ambos selectores basados ​​en datos son estimadores consistentes.

El paquete ks ^{[ 13 ]} en R implementa los selectores de validación cruzada suavizada y de complemento (entre otros). Este conjunto de datos (incluido en la distribución base de R) contiene 272 registros con dos mediciones cada uno: la duración de una erupción (en minutos) y el tiempo de espera hasta la siguiente erupción (en minutos) del géiser Old Faithful en el Parque Nacional de Yellowstone, EE. UU.

El fragmento de código calcula la estimación de la densidad del kernel con la matriz de ancho de banda de entrada.${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {PI} }={\begin{bmatrix}0.052&0.510\\0.510&8.882\end{bmatrix}}.}$Nuevamente, los contornos de color corresponden a la región más pequeña que contiene la masa de probabilidad respectiva: rojo = 25%, naranja + rojo = 50%, amarillo + naranja + rojo = 75%. Para calcular el selector SCV, Hpise reemplaza con Hscv. Esto no se muestra aquí ya que es muy similar a la estimación de la entrada para este ejemplo.

biblioteca ( ks ) datos ( faithful ) H <- Hpi ( x = faithful ) fhat <- kde ( x = faithful , H = H ) plot ( fhat , display = "filled.contour" , drawpoints = TRUE , cex = 0.5 , pch = 16 , col.pt = 1 )

Consideramos estimar la densidad de la mezcla gaussiana (4 π ) ⁻¹ exp(− 1 ⁄ 2 ( x ₁ ² + x ^{2 2} ₎ ) + (4 π ) ⁻¹ exp(− 1 ⁄ 2 (( x ₁ - 3.5) ² + x ₂ ² )) , a partir de 500 puntos generados aleatoriamente. Empleamos la rutina de Matlab para datos bidimensionales . La rutina es un método de selección automática de ancho de banda diseñado específicamente para un núcleo gaussiano de segundo orden. ^{[ 14 ]} La figura muestra la estimación de densidad conjunta que resulta de usar el ancho de banda seleccionado automáticamente.

Script de Matlab para el ejemplo

Escriba los siguientes comandos en Matlab después de descargar y guardar la función kde2d.m en el directorio actual.

clear all % generar datos sintéticos data =[ randn ( 500 , 2 ); randn ( 500 , 1 ) + 3.5 , randn ( 500 , 1 );]; % llamar a la rutina, que se ha guardado en el directorio actual [ ancho de banda , densidad , X , Y ] = kde2d ( datos ); % graficar los datos y la estimación de densidad contour3 ( X , Y , densidad , 50 ), hold on plot ( datos (:, 1 ), datos (:, 2 ), 'r.' , 'MarkerSize' , 5 )

El MISE es la distancia L2 integrada esperada entre la estimación de densidad y la función de densidad verdadera f . Es la más utilizada, principalmente debido a su manejabilidad y a que la mayoría del software implementa selectores de ancho de banda basados ​​en MISE. Existen criterios de optimalidad alternativos que intentan cubrir casos en los que MISE no es una medida apropiada. _[ ^{4 ] : ​​34–37 , 78} La medida L1 equivalente , _el Error Absoluto Integrado Medio, es

${\displaystyle \operatorname {MIAE} (\mathbf {H} )=\operatorname {E} \,\int |{\hat {f}}_{\mathbf {H} }(\mathbf {x} )-f(\mathbf {x} )|\,d\mathbf {x} .}$

Su análisis matemático es considerablemente más difícil que el de los MISE. En la práctica, la ganancia no parece ser significativa. ^{[ 15 ]} La norma L _∞ es el error absoluto uniforme medio

${\displaystyle \operatorname {MUAE} (\mathbf {H} )=\operatorname {E} \,\operatorname {sup} _{\mathbf {x} }|{\hat {f}}_{\mathbf {H} }(\mathbf {x} )-f(\mathbf {x} )|.}$

que se ha investigado solo brevemente. ^{[ 16 ]} Los criterios de error de verosimilitud incluyen aquellos basados ​​en la divergencia media de Kullback-Leibler

${\displaystyle \operatorname {MKL} (\mathbf {H} )=\int f(\mathbf {x} )\,\operatorname {log} [f(\mathbf {x} )]\,d\mathbf {x} -\operatorname {E} \int f(\mathbf {x} )\,\operatorname {log} [{\hat {f}}(\mathbf {x}$ ;\mathbf {H} )]\,d\mathbf {x} }

y la distancia media de Hellinger

${\displaystyle \operatorname {MH} (\mathbf {H} )=\operatorname {E} \int ({\hat {f}}_{\mathbf {H} }(\mathbf {x} )^{1/2}-f(\mathbf {x} )^{1/2})^{2}\,d\mathbf {x} .}$

El KL se puede estimar utilizando un método de validación cruzada, aunque los selectores de validación cruzada de KL pueden ser subóptimos incluso si se mantienen consistentes para funciones de densidad acotadas. ^{[ 17 ]} Los selectores MH se han examinado brevemente en la literatura. ^{[ 18 ]}

Todos estos criterios de optimización son medidas basadas en la distancia y no siempre se corresponden con nociones más intuitivas de cercanía, por lo que se han desarrollado criterios más visuales en respuesta a esta preocupación. ^{[ 19 ]}

Investigaciones recientes han demostrado que tanto el núcleo como su ancho de banda pueden elegirse de manera óptima y objetiva a partir de los propios datos de entrada sin hacer suposiciones sobre la forma de la distribución. ^{[ 20 ]} La estimación de densidad del núcleo resultante converge rápidamente a la verdadera distribución de probabilidad a medida que se agregan muestras: a una tasa cercana a la${\displaystyle n^{-1}}$esperado para estimadores paramétricos. ^{[ 20 ] [ 21 ] [ 22 ]} Este estimador de núcleo funciona tanto para muestras univariadas como multivariadas. El núcleo óptimo se define en el espacio de Fourier, como la función de amortiguación óptima.${\displaystyle {\hat {\psi _{h}}}({\vec {t}})}$(la transformada de Fourier del núcleo${\displaystyle {\hat {K}}({\vec {x}})}$)-- en términos de la transformada de Fourier de los datos${\displaystyle {\hat {\varphi }}({\vec {t}})}$, la función característica empírica (véase Estimación de densidad de núcleo ):

${\displaystyle {\hat {\psi _{h}}}({\vec {t}})\equiv {\frac {N}{2(N-1)}}\left[1+{\sqrt {1-{\frac {4(N-1)}{N^{2}|{\hat {\varphi }}({\vec {t}})|^{2}}}}}I_{\vec {A}}({\vec {t}})\right]}$^{[ 22 ]}

${\displaystyle {\hat {f}}(x)={\frac {1}{(2\pi )^{d}}}\int {\hat {\varphi }}({\vec {t}})\psi _{h}({\vec {t}})e^{-i{\vec {t}}\cdot {\vec {x}}}d{\vec {t}}}$

donde N es el número de puntos de datos, d es el número de dimensiones (variables) y${\displaystyle I_{\vec {A}}({\vec {t}})}$es un filtro que es igual a 1 para las "frecuencias aceptadas" y 0 en caso contrario. Hay varias formas de definir esta función de filtro, y una simple que funciona para muestras univariadas o multivariadas se llama "filtro de hipervolumen contiguo más bajo";${\displaystyle I_{\vec {A}}({\vec {t}})}$se elige de tal manera que las únicas frecuencias aceptadas sean un subconjunto contiguo de frecuencias que rodean el origen para el cual${\displaystyle |{\hat {\varphi }}({\vec {t}})|^{2}\geq 4(N-1)N^{-2}}$(véase ^{[ 22 ]} para un análisis de esta y otras funciones de filtro).

Cabe señalar que el cálculo directo de la función característica empírica (FCE) es lento, ya que implica esencialmente una transformada de Fourier directa de las muestras de datos. Sin embargo, se ha descubierto que la FCE puede aproximarse con precisión utilizando un método de transformada rápida de Fourier no uniforme (nuFFT), ^{[ 21 ] [ 22 ]} que aumenta la velocidad de cálculo en varios órdenes de magnitud (dependiendo de la dimensionalidad del problema). La combinación de este método KDE objetivo y la aproximación de la FCE basada en nuFFT se ha denominado fastKDE en la literatura. ^{[ 22 ]}

• Estimación de densidad de núcleo – estimación de densidad de núcleo univariada. 
• Estimación de densidad de núcleo variable : estimación de densidades multivariadas utilizando el núcleo con ancho de banda variable 

• Suavizado de núcleo multivariado y sus aplicaciones es un libro exhaustivo sobre diversos temas relacionados con el suavizado de núcleo, incluida la estimación de densidad. Incluyefragmentosde código del paquete ks en R.
• kde2d.m Una función de Matlab para la estimación de la densidad del núcleo bivariado.
• libagf Una biblioteca de C++ para la estimación de densidad de kernel multivariada de ancho de banda variable .
• akde.m Un archivo m de Matlab para la estimación de densidad de kernel multivariada de ancho de banda variable .
• helit y el módulo pyqt_fit.kde del paquete PyQt-Fit son bibliotecas de Python para la estimación de la densidad del núcleo multivariante.