Estimation de la densité du noyau - Kernel density estimation

Estimation de la densité du noyau de 100 nombres aléatoires normalement distribués en utilisant différentes bandes passantes de lissage.

En statistique , l' estimation de la densité par noyau ( KDE ) est un moyen non paramétrique d' estimer la fonction de densité de probabilité d' une variable aléatoire . L'estimation de la densité du noyau est un problème fondamental de lissage des données où des inférences sur la population sont faites, sur la base d'un échantillon de données fini . Dans certains domaines tels que le traitement du signal et l' économétrie, elle est également appelée méthode de la fenêtre de Parzen-Rosenblatt , d'après Emanuel Parzen et Murray Rosenblatt , qui sont généralement crédités de sa création indépendante sous sa forme actuelle. L'une des applications célèbres de l'estimation de la densité par noyau consiste à estimer les densités marginales conditionnelles de classe des données lors de l'utilisation d'un classificateur naïf de Bayes , ce qui peut améliorer la précision de sa prédiction.

Définition

Soit ( x ₁ , x ₂ , ..., x _n ) être indépendantes et identiquement distribuées des échantillons tirés de la distribution univariée avec une inconnue une densité ƒ à un moment donné x . Nous nous intéressons à l'estimation de la forme de cette fonction ƒ . Son estimateur de densité à noyau est

{\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 ) },

où K est le noyau — une fonction non négative — et h > 0 est un paramètre de lissage appelé bande passante. Un noyau avec l'indice h est appelé le noyau mis à l' échelle et défini comme K _h ( x ) = 1/ h K ( x / h ) . Intuitivement, on veut choisir h aussi petit que les données le permettent ; cependant, il y a toujours un compromis entre le biais de l'estimateur et sa variance. Le choix de la bande passante est discuté plus en détail ci-dessous.

Une gamme de fonctions du noyau est couramment utilisée : uniforme, triangulaire, bipoids, tripoids, Epanechnikov, normal et autres. Le noyau d'Epanechnikov est optimal au sens de l'erreur quadratique moyenne, bien que la perte d'efficacité soit faible pour les noyaux énumérés précédemment. En raison de ses propriétés mathématiques pratiques, le noyau normal est souvent utilisé, ce qui signifie que K ( x ) = φ ( x ) , où φ est la norme normale fonction de densité.

La construction d'une estimation de densité de noyau trouve des interprétations dans des domaines en dehors de l'estimation de densité. Par exemple, en thermodynamique , cela équivaut à la quantité de chaleur générée lorsque des noyaux de chaleur (la solution fondamentale de l' équation de la chaleur ) sont placés à chaque emplacement de point de données x _i . Des méthodes similaires sont utilisées pour construire des opérateurs de Laplace discrets sur des nuages de points pour l' apprentissage multiple (par exemple, carte de diffusion ).

Exemple

Les estimations de densité de noyau sont étroitement liées aux histogrammes , mais peuvent être dotées de propriétés telles que la régularité ou la continuité en utilisant un noyau approprié. Le diagramme ci-dessous basé sur ces 6 points de données illustre cette relation :

Échantillon	1	2	3	4	5	6
Valeur	-2.1	-1,3	-0,4	1.9	5.1	6.2

Pour l'histogramme, tout d'abord, l'axe horizontal est divisé en sous-intervalles ou cases qui couvrent la plage des données : dans ce cas, six cases chacune de largeur 2. Chaque fois qu'un point de données tombe à l'intérieur de cet intervalle, une case de hauteur 1 /12 y est placé. Si plusieurs points de données tombent dans le même bac, les boîtes sont empilées les unes sur les autres.

Pour l'estimation de la densité par noyau, des noyaux normaux avec un écart type de 2,25 (indiqués par les lignes pointillées rouges) sont placés sur chacun des points de données x _i . Les noyaux sont additionnés pour faire l'estimation de la densité du noyau (courbe bleue continue). La régularité de l'estimation de la densité du noyau (par rapport au caractère discret de l'histogramme) illustre comment les estimations de la densité du noyau convergent plus rapidement vers la véritable densité sous-jacente pour les variables aléatoires continues.

Comparaison de l'histogramme (à gauche) et de l'estimation de la densité du noyau (à droite) construits à partir des mêmes données. Les six noyaux individuels sont les courbes en pointillés rouges, la densité du noyau estime les courbes bleues. Les points de données sont le tracé du tapis sur l'axe horizontal.

Sélection de la bande passante

Estimation de la densité du noyau (KDE) avec différentes bandes passantes d'un échantillon aléatoire de 100 points à partir d'une distribution normale standard. Gris : densité réelle (standard normal). Rouge : KDE avec h=0,05. Noir : KDE avec h=0.337. Vert : KDE avec h=2.

La bande passante du noyau est un paramètre libre qui a une forte influence sur l'estimation résultante. Pour illustrer son effet, nous prenons un échantillon aléatoire simulé de la distribution normale standard (tracé au niveau des pointes bleues dans le tracé du tapis sur l'axe horizontal). La courbe grise est la vraie densité (une densité normale avec une moyenne 0 et une variance 1). En comparaison, la courbe rouge est sous- lissée car elle contient trop d'artefacts de données parasites résultant de l'utilisation d'une bande passante h = 0,05, ce qui est trop petit. La courbe verte est trop lissée car l'utilisation de la bande passante h = 2 masque une grande partie de la structure sous-jacente. La courbe noire avec une largeur de bande de h = 0,337 est considérée comme lissée de manière optimale puisque son estimation de densité est proche de la densité réelle. Une situation extrême est rencontrée à la limite (pas de lissage), où l'estimation est une somme de n fonctions delta centrées aux coordonnées des échantillons analysés. A l'autre extrême limite l'estimation conserve la forme du noyau utilisé, centrée sur la moyenne des échantillons (complètement lisse). ${\style d'affichage h\à 0}$ $h\to \infty$

Le critère d'optimalité le plus couramment utilisé pour sélectionner ce paramètre est la fonction de risque attendue L ₂ , également appelée erreur quadratique moyenne intégrée :

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

Sous des hypothèses faibles sur ƒ et K , ( ƒ est la fonction de densité réelle, généralement inconnue), MISE ( h ) = AMISE( h ) + o(1/(nh) + h ⁴ ) où o est la petite notation o , et n la taille de l'échantillon (comme ci-dessus). L'AMISE est l'Asymptotique MISE qui se compose des deux termes principaux

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

où pour une fonction g , et est la dérivée seconde de . Le minimum de cette AMISE est la solution de cette équation différentielle $R(g)=\int g(x)^{2}\,dx$ $m_{2}(K)=\int x^{2}K(x)\,dx$ ${\style d'affichage f''}$ ${\style d'affichage f}$

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

ou

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

Ni les formules AMISE ni h _AMISE ne peuvent être utilisées directement car elles impliquent la fonction de densité inconnue ou sa dérivée seconde , de sorte qu'une variété de méthodes automatiques basées sur les données ont été développées pour sélectionner la bande passante. De nombreuses études de revue ont été menées pour comparer leurs efficacités, avec le consensus général que les sélecteurs de plug-in et les sélecteurs de validation croisée sont les plus utiles sur un large éventail d'ensembles de données. ${\style d'affichage f}$ ${\style d'affichage f''}$

La substitution de toute bande passante h qui a le même ordre asymptotique n ^-1/5 que h _AMISE dans l'AMISE donne que AMISE( h ) = O ( n ^-4/5 ), où O est la grande notation o . On peut montrer que, sous des hypothèses faibles, il ne peut pas exister d'estimateur non paramétrique qui converge à un rythme plus rapide que l'estimateur à noyau. Notez que le taux n ^-4/5 est plus lent que le taux de convergence typique n ^-1 des méthodes paramétriques.

Si la bande passante n'est pas maintenue fixe, mais varie en fonction de l'emplacement de l'estimation (estimateur par ballonnet) ou des échantillons (estimateur ponctuel), cela produit une méthode particulièrement puissante appelée estimation de la densité du noyau à bande passante adaptative ou variable .

La sélection de la bande passante pour l'estimation de la densité du noyau des distributions à queue lourde est relativement difficile.

Un estimateur de bande passante empirique

Si des fonctions de base gaussiennes sont utilisées pour approximer des données univariées et que la densité sous-jacente estimée est gaussienne, le choix optimal pour h (c'est-à-dire la bande passante qui minimise l' erreur quadratique moyenne intégrée ) est :

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

Afin de rendre la valeur h plus robuste afin de bien adapter la distribution à la fois à longue traîne et asymétrique et à la distribution de mélange bimodal, il est préférable de substituer la valeur de par un autre paramètre A, qui est donné par : ${\hat {\sigma }}$

A = min (écart type, intervalle interquartile / 1,34).

Une autre modification qui améliorera le modèle est de réduire le facteur de 1,06 à 0,9. Alors la formule finale serait :

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

où est l' écart type des échantillons, n est la taille de l'échantillon. IQR est l'intervalle interquartile. ${\hat {\sigma }}$

Cette approximation est appelée approximation de la distribution normale , approximation gaussienne ou règle empirique de Silverman . Bien que cette règle empirique soit facile à calculer, elle doit être utilisée avec prudence car elle peut donner des estimations largement inexactes lorsque la densité n'est pas proche d'être normale. Par exemple, lors de l'estimation du modèle de mélange gaussien bimodal

Comparaison entre la règle empirique et la bande passante de résolution de l'équation.

\textstyle {\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}}

à partir d'un échantillon de 200 points. La figure de droite montre la densité réelle et deux estimations de densité de noyau, l'une utilisant la bande passante empirique et l'autre utilisant une bande passante de résolution de l'équation. L'estimation basée sur la bande passante de la règle empirique est considérablement surlissée.

Relation avec l'estimateur de densité de fonction caractéristique

Étant donné l'échantillon ( x ₁ , x ₂ , …, x _n ), il est naturel d'estimer la fonction caractéristique φ ( t ) = E[ e ^itX ] comme

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

Connaissant la fonction caractéristique, il est possible de trouver la fonction de densité de probabilité correspondante grâce à la formule de la transformée de Fourier . Une difficulté avec l'application de cette formule d'inversion est qu'elle conduit à une intégrale divergente, car l'estimation n'est pas fiable pour les grands t . Pour contourner ce problème, l'estimateur est multiplié par une fonction d'amortissement ψ _h ( t ) = ψ ( ht ) , qui est égale à 1 à l'origine puis tombe à 0 à l'infini. Le "paramètre de bande passante" h contrôle la vitesse à laquelle nous essayons d'atténuer la fonction . En particulier, lorsque h est petit, alors ψ _h ( t ) sera approximativement égal à un pour une large plage de t , ce qui signifie qu'il reste pratiquement inchangé dans la région la plus importante de t . $\scriptstyle {\widehat {\varphi }}(t)$ $\scriptstyle {\widehat {\varphi }}(t)$ $\scriptstyle {\widehat {\varphi }}(t)$ $\scriptstyle {\widehat {\varphi }}(t)$

Le choix le plus courant pour la fonction ψ est soit la fonction uniforme ψ ( t ) = 1 {−1 ≤ t ≤ 1 }, ce qui signifie effectivement tronquer l'intervalle d'intégration dans la formule d'inversion à [−1/ h , 1/ h ] , soit la fonction gaussienne ψ ( t ) = e ^{− $π$ t ²} . Une fois la fonction ψ choisie, la formule d'inversion peut être appliquée, et l'estimateur de densité sera

{\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}}

où K est la transformée de Fourier de la fonction d' amortissement ψ . Ainsi, l'estimateur de densité à noyau coïncide avec l'estimateur de densité de fonction caractéristique.

Caractéristiques géométriques et topologiques

On peut étendre la définition du mode (global) à un sens local et définir les modes locaux :

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

À savoir, est la collection de points pour lesquels la fonction de densité est localement maximisée. Un estimateur naturel de est un plug-in de KDE, où et sont la version KDE de et . Sous des hypothèses modérées, est un estimateur cohérent de . Notez que l'on peut utiliser l'algorithme de décalage moyen pour calculer l'estimateur numériquement. ${\style d'affichage M}$ ${\style d'affichage M}$ ${\style d'affichage g(x)}$ $\lambda _{1}(x)$ ${\style d'affichage g(x)}$ $\lambda _{1}(x)$ ${\style d'affichage M_{c}}$ ${\style d'affichage M}$ ${\style d'affichage M_{c}}$

Mise en œuvre statistique

Une liste non exhaustive d'implémentations logicielles d'estimateurs de densité à noyau comprend :

Dans Analytica version 4.4, l' option Lissage pour les résultats PDF utilise KDE, et à partir des expressions, elle est disponible via la Pdffonction intégrée.
En C / C++ , FIGTree est une bibliothèque qui peut être utilisée pour calculer des estimations de densité de noyau en utilisant des noyaux normaux. Interface MATLAB disponible.
En C++ , libagf est une bibliothèque d' estimation de densité de noyau variable .
En C++ , mlpack est une bibliothèque qui peut calculer KDE en utilisant de nombreux noyaux différents. Il permet de définir une tolérance d'erreur pour un calcul plus rapide. Des interfaces Python et R sont disponibles.
en C# et F# , Math.NET Numerics est une bibliothèque open source pour le calcul numérique qui inclut l' estimation de la densité du noyau
Dans CrimeStat , l'estimation de la densité par noyau est mise en œuvre à l'aide de cinq fonctions de noyau différentes : normale, uniforme, quartique, exponentielle négative et triangulaire. Des routines d'estimation de densité à noyau unique et à noyau double sont disponibles. L'estimation de la densité du noyau est également utilisée pour interpoler une routine Head Bang, pour estimer une fonction de densité bidimensionnelle du trajet jusqu'au crime et pour estimer une estimation bayésienne tridimensionnelle du trajet jusqu'au crime.
Dans ELKI , les fonctions de densité du noyau peuvent être trouvées dans le packagede.lmu.ifi.dbs.elki.math.statistics.kernelfunctions
Dans les produits ESRI , la cartographie de la densité du noyau est gérée à partir de la boîte à outils Spatial Analyst et utilise le noyau Quartic (biweight).
Dans Excel , la Royal Society of Chemistry a créé un complément pour exécuter l'estimation de la densité du noyau sur la base de son dossier technique du Comité des méthodes analytiques 4 .
Dans gnuplot , l'estimation de la densité du noyau est implémentée par l' smooth kdensityoption, le fichier de données peut contenir un poids et une bande passante pour chaque point, ou la bande passante peut être définie automatiquement selon la "règle empirique de Silverman" (voir ci-dessus).
Dans Haskell , la densité du noyau est implémentée dans le package de statistiques .
Dans IGOR Pro , l'estimation de la densité du noyau est implémentée par l' StatsKDEopération (ajoutée dans Igor Pro 7.00). La bande passante peut être spécifiée ou estimée par l'utilisateur au moyen de Silverman, Scott ou Bowmann et Azzalini. Les types de noyau sont : Epanechnikov, Bi-poids, Tri-poids, Triangulaire, Gaussien et Rectangulaire.
En Java , le package Weka (machine learning) fournit entre autres weka.estimators.KernelEstimator .
En JavaScript , le package de visualisation D3.js propose un package KDE dans son package science.stats.
Dans JMP , la plate-forme Graph Builder utilise l'estimation de la densité du noyau pour fournir des tracés de contour et des régions à haute densité (HDR) pour les densités bivariées, ainsi que des tracés de violon et des HDR pour les densités univariées. Les curseurs permettent à l'utilisateur de faire varier la bande passante. Des estimations de densité de noyau bivariées et univariées sont également fournies par les plates-formes Fit Y by X et Distribution, respectivement.
Dans Julia , l'estimation de la densité du noyau est implémentée dans le package KernelDensity.jl .
Dans MATLAB , l'estimation de la densité du noyau est implémentée via la ksdensityfonction (Statistics Toolbox). À partir de la version 2018a de MATLAB, la bande passante et le lissage du noyau peuvent être spécifiés, y compris d'autres options telles que la spécification de la plage de densité du noyau. Alternativement, un progiciel MATLAB gratuit qui implémente une méthode de sélection automatique de la bande passante est disponible à partir de MATLAB Central File Exchange pour
- Données unidimensionnelles
- Données bidimensionnelles
- Données n-dimensionnelles
  Une boîte à outils MATLAB gratuite avec implémentation de la régression du noyau, de l'estimation de la densité du noyau, de l'estimation du noyau de la fonction de risque et bien d'autres est disponible sur ces pages (cette boîte à outils fait partie du livre).
Dans Mathematica , l'estimation de la densité du noyau numérique est implémentée par la fonction SmoothKernelDistributionet l'estimation symbolique est implémentée à l'aide de la fonction, KernelMixtureDistributionqui fournissent toutes deux des bandes passantes pilotées par les données.
Dans Minitab , la Royal Society of Chemistry a créé une macro pour exécuter l'estimation de la densité du noyau sur la base de la note technique 4 du Comité des méthodes analytiques.
Dans la bibliothèque NAG , l'estimation de la densité du noyau est implémentée via la g10baroutine (disponible à la fois dans les versions Fortran et C de la bibliothèque).
Dans Nuklei , les méthodes de densité du noyau C++ se concentrent sur les données du groupe euclidien spécial . $SE(3)$
Dans Octave , l'estimation de la densité du noyau est implémentée par l' kernel_densityoption (package économétrie).
Dans origine , le tracé de la densité de 2D peut être fabriqué à partir de son interface utilisateur, et deux fonctions, pour Ksdensity 1D et 2D pour Ks2density peut être utilisé à partir de son LabTalk , Python ou C code.
En Perl , une implémentation peut être trouvée dans le module Statistics-KernelEstimation
En PHP , une implémentation peut être trouvée dans la bibliothèque MathPHP
En Python , de nombreuses implémentations existent : pyqt_fit.kde Module dans le package PyQt-Fit , SciPy ( scipy.stats.gaussian_kde), Statsmodels ( KDEUnivariateet KDEMultivariate) et Scikit-learn ( KernelDensity) (voir comparaison). KDEpy prend en charge les données pondérées et son implémentation FFT est de plusieurs ordres de grandeur plus rapide que les autres implémentations. La bibliothèque pandas couramment utilisée [1] offre un support pour le traçage de kde via la méthode plot ( df.plot(kind='kde')[2] ). Le package getdist pour les échantillons MCMC pondérés et corrélés prend en charge la bande passante optimisée, la correction des limites et les méthodes d'ordre supérieur pour les distributions 1D et 2D. Un package nouvellement utilisé pour l'estimation de la densité du noyau est seaborn ( import seaborn as sns, sns.kdeplot()). Une implémentation GPU de KDE existe également.
Dans R , il est implémenté densitydans la distribution de base et la bw.nrd0fonction est utilisée dans le package de statistiques, cette fonction utilise la formule optimisée du livre de Silverman. bkdedans la bibliothèque KernSmooth , ParetoDensityEstimationdans la bibliothèque de DataVisualizations (pour l' estimation de la densité de distribution de Pareto), kdedans la bibliothèque de ks , dkdenet dbckdendans la bibliothèque de evmix (ci pour limite l' estimation de la densité du noyau corrigé pour support borné), npudensdans la bibliothèque de np (données numérique et catégorique) , sm.densitydans la bibliothèque sm . Pour une implémentation de la kde.Rfonction, qui ne nécessite pas l'installation de packages ou de bibliothèques, consultez kde.R . La bibliothèque btb , dédiée à l'analyse urbaine, implémente l'estimation de la densité du noyau via kernel_smoothing.
Dans SAS , proc kdepeut être utilisé pour estimer les densités de noyau univariées et bivariées.
Dans Apache Spark , la KernelDensity()classe
Dans Stata , il est mis en œuvre par kdensity; par exemple histogram x, kdensity. Alternativement, un module Stata gratuit KDENS est disponible à partir d' ici permettant à un utilisateur d'estimer les fonctions de densité 1D ou 2D.
Dans Swift , il est implémenté SwiftStats.KernelDensityEstimationdans la bibliothèque de statistiques open source SwiftStats .

Voir également

Noyau (statistiques)
Lissage du noyau
Régression du noyau
Estimation de la densité (avec présentation d'autres exemples)
Décalage moyen
Espace d'échelle : Les triplets {( x , h , KDE avec une bande passante h évaluée à x : tous x , h > 0} forment une représentation d' espace d'échelle des données.
Estimation de la densité de noyau multivariée
Estimation de la densité du noyau variable
Casse tête/queue

Les références

Liens externes

Introduction à l'estimation de densité à noyau Un court tutoriel qui motive les estimateurs de densité à noyau comme une amélioration par rapport aux histogrammes.
Optimisation de la bande passante du noyau Un outil en ligne gratuit qui génère une estimation optimisée de la densité du noyau.
Un logiciel en ligne gratuit (Calculatrice) calcule l'estimation de la densité du noyau pour une série de données selon les noyaux suivants : gaussien, épanechnikov, rectangulaire, triangulaire, bipoids, cosinus et optcosinus.
Applet d'estimation de la densité du noyau Un exemple interactif en ligne d'estimation de la densité du noyau. Nécessite .NET 3.0 ou version ultérieure.

Languages

In other projects