Kernel density estimation

Example

Kernel density estimates are closely related to histograms, but can be endowed with properties such as smoothness or continuity by using a suitable kernel. The diagram below based on these 6 data points illustrates this relationship:

For the histogram, first, the horizontal axis is divided into sub-intervals or bins which cover the range of the data: In this case, six bins each of width 2. Whenever a data point falls inside this interval, a box of height 1/12 is placed there. If more than one data point falls inside the same bin, the boxes are stacked on top of each other.

For the kernel density estimate, normal kernels with a standard deviation of 1.5 (indicated by the red dashed lines) are placed on each of the data points x_i. The kernels are summed to make the kernel density estimate (solid blue curve). The smoothness of the kernel density estimate (compared to the discreteness of the histogram) illustrates how kernel density estimates converge faster to the true underlying density for continuous random variables.

Bandwidth selection

The bandwidth of the kernel is a free parameter which exhibits a strong influence on the resulting estimate. To illustrate its effect, we take a simulated random sample from the standard normal distribution (plotted at the blue spikes in the rug plot on the horizontal axis). The grey curve is the true density (a normal density with mean 0 and variance 1). In comparison, the red curve is undersmoothed since it contains too many spurious data artifacts arising from using a bandwidth h = 0.05, which is too small. The green curve is oversmoothed since using the bandwidth h = 2 obscures much of the underlying structure. The black curve with a bandwidth of h = 0.337 is considered to be optimally smoothed since its density estimate is close to the true density. An extreme situation is encountered in the limit ${\displaystyle h\to 0}$ (no smoothing), where the estimate is a sum of n delta functions centered at the coordinates of analyzed samples. In the other extreme limit ${\displaystyle h\to \mathrm{\infty }}$ the estimate retains the shape of the used kernel, centered on the mean of the samples (completely smooth).

The most common optimality criterion used to select this parameter is the expected L₂ risk function, also termed the mean integrated squared error:

$${\displaystyle \mathrm{MISE}(h)=\mathrm{E}\left[\int {\left(\hat{f}{}_h(x)-f(x)\right)}^{2}dx\right]}$$

Under weak assumptions on f and K, (f is the, generally unknown, real density function),

$${\displaystyle \mathrm{MISE}(h)=\mathrm{AMISE}(h)+\mathcal{o}\left((nh{)}^{-1}+h^4\right)}$$

where o is the little o notation, and n the sample size (as above). The AMISE is the asymptotic MISE, i. e. the two leading terms,

$${\displaystyle \mathrm{AMISE}(h)=\frac{R(K)}{nh}+\frac{1}{4}{m}_2(K{)}^2{h}^{4}R(f^{″})}$$

where $R(g)=\int g(x{)}^{2}dx$ for a function g, $m_2(K)=\int x^{2}K(x)dx$ and ${\displaystyle f^{″}}$ is the second derivative of ${\displaystyle f}$ and ${\displaystyle K}$ is the kernel. The minimum of this AMISE is the solution to this differential equation

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }h}\mathrm{AMISE}(h)=-\frac{R(K)}{n{h}^2}+m_2(K{)}^2{h}^{3}R(f^{″})=0}$$

or

$${\displaystyle h_{\mathrm{AMISE}}=\frac{R(K{)}^{1/5}}{m_2(K{)}^{2/5}R(f^{″}{)}^{1/5}}{n}^{-1/5}=C{n}^{-1/5}}$$

Neither the AMISE nor the h_AMISE formulas can be used directly since they involve the unknown density function ${\displaystyle f}$ or its second derivative ${\displaystyle f^{″}}$. To overcome that difficulty, a variety of automatic, data-based methods have been developed to select the bandwidth. Several review studies have been undertaken to compare their efficacies, with the general consensus that the plug-in selectors and cross validation selectors are the most useful over a wide range of data sets.

Substituting any bandwidth h which has the same asymptotic order n^−1/5 as h_AMISE into the AMISE gives that AMISE(h) = O(n^−4/5), where O is the big O notation. It can be shown that, under weak assumptions, there cannot exist a non-parametric estimator that converges at a faster rate than the kernel estimator. Note that the n^−4/5 rate is slower than the typical n⁻¹ convergence rate of parametric methods.

If the bandwidth is not held fixed, but is varied depending upon the location of either the estimate (balloon estimator) or the samples (pointwise estimator), this produces a particularly powerful method termed adaptive or variable bandwidth kernel density estimation.

Bandwidth selection for kernel density estimation of heavy-tailed distributions is relatively difficult.

A rule-of-thumb bandwidth estimator

If Gaussian basis functions are used to approximate univariate data, and the underlying density being estimated is Gaussian, the optimal choice for h (that is, the bandwidth that minimises the mean integrated squared error) is:

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

An ${\displaystyle h}$ value is considered more robust when it improves the fit for long-tailed and skewed distributions or for bimodal mixture distributions. This is often done empirically by replacing the standard deviation ${\displaystyle \hat{\sigma }}$ by the parameter ${\displaystyle A}$ below:

$${\displaystyle A=min\left(\hat{\sigma },\frac{\mathrm{I}\mathrm{Q}\mathrm{R}}{1.34}\right)}$$ where IQR is the interquartile range.

Another modification that will improve the model is to reduce the factor from 1.06 to 0.9. Then the final formula would be:

$${\displaystyle h=0.9min\left(\hat{\sigma },\frac{\mathrm{I}\mathrm{Q}\mathrm{R}}{1.34}\right)n^{-1/5}}$$ where ${\displaystyle n}$ is the sample size.

This approximation is termed the normal distribution approximation, Gaussian approximation, or Silverman's rule of thumb. While this rule of thumb is easy to compute, it should be used with caution as it can yield widely inaccurate estimates when the density is not close to being normal. For example, when estimating the bimodal Gaussian mixture model $${\displaystyle \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}}$$ from a sample of 200 points, the figure on the right shows the true density and two kernel density estimates — one using the rule-of-thumb bandwidth, and the other using a solve-the-equation bandwidth. The estimate based on the rule-of-thumb bandwidth is significantly oversmoothed.

Relation to the characteristic function density estimator

Given the sample (x₁, x₂, ..., x_n), it is natural to estimate the characteristic function φ(t) = E[e^itX] as $${\displaystyle \hat{\phi }(t)=\frac{1}{n}\sum _{j=1}^n{e}^{it{x}_j}}$$ Knowing the characteristic function, it is possible to find the corresponding probability density function through the Fourier transform formula. One difficulty with applying this inversion formula is that it leads to a diverging integral, since the estimate ${\displaystyle \hat{\phi }(t)}$ is unreliable for large t's. To circumvent this problem, the estimator ${\displaystyle \hat{\phi }(t)}$ is multiplied by a damping function ψ_h(t) = ψ(ht), which is equal to 1 at the origin and then falls to 0 at infinity. The "bandwidth parameter" h controls how fast we try to dampen the function ${\displaystyle \hat{\phi }(t)}$. In particular when h is small, then ψ_h(t) will be approximately one for a large range of t's, which means that ${\displaystyle \hat{\phi }(t)}$ remains practically unaltered in the most important region of t's.

The most common choice for function ψ is either the uniform function ψ(t) = 1{−1 ≤ t ≤ 1}, which effectively means truncating the interval of integration in the inversion formula to [−1/h, 1/h], or the Gaussian function ψ(t) = e^−πt2. Once the function ψ has been chosen, the inversion formula may be applied, and the density estimator will be $${\displaystyle \begin{array}{cc}\hat{f}(x)& =\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\hat{\phi }(t){\psi }_h(t)e^{-itx}dt\\ & =\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{1}{n}\sum _{j=1}^n{e}^{it(x_j-x)}\psi (ht)dt\\ & =\frac{1}{nh}\sum _{j=1}^n\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{e}^{-i(ht)\frac{x-x_j}{h}}\psi (ht)d(ht)\\ & =\frac{1}{nh}\sum _{j=1}^{n}K\left(\frac{x-x_j}{h}\right),\end{array}}$$

where K is the Fourier transform of the damping function ψ. Thus the kernel density estimator coincides with the characteristic function density estimator.

Kernel quantum state estimation

Kernel density estimation has been applied to continuous-variable quantum state tomography, in which repeated homodyne measurements produce samples of rotated field quadratures. For real measurement parameters ${\displaystyle (\mu ,\nu )}$, the quadrature observable is

${\displaystyle X_{\mu ,\nu }=\mu q+\nu p.}$

Its measurement outcomes follow an ordinary probability density ${\displaystyle \mathcal{W}(x\mid \mu ,\nu )}$, known as a symplectic tomogram. The collection of tomograms over different measurement settings is informationally complete and therefore provides an alternative representation of the quantum state.

For independent quadrature measurements ${\displaystyle X_{1,\mu ,\nu },\dots ,X_{n,\mu ,\nu }}$ obtained at a fixed setting, the tomogram can be estimated by

${\displaystyle {\hat{\mathcal{W}}}_{n,h}(x\mid \mu ,\nu )=\frac{1}{nh}\sum _{j=1}^n{K}_{\mu ,\nu }\left(\frac{x-X_{j,\mu ,\nu }}{h}\right),}$

where ${\displaystyle K_{\mu ,\nu }}$ is a kernel and ${\displaystyle h}$ is the bandwidth. This treatment does not require the quadrature density to belong to a predetermined parametric family and can therefore accommodate multimodal and non-Gaussian measurement distributions.

The characteristic function of a tomogram is its Fourier transform,

${\displaystyle \phi (t;\mu ,\nu )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathcal{W}(x\mid \mu ,\nu )e^{itx}dx.}$

For every fixed ${\displaystyle (\mu ,\nu )}$, this is the characteristic function of the measured random variable ${\displaystyle X_{\mu ,\nu }}$. It always exists and satisfies the standard conditions ${\displaystyle \phi (0;\mu ,\nu )=1}$, ${\displaystyle |\phi (t;\mu ,\nu )|\le 1}$, and ${\displaystyle \phi (-t;\mu ,\nu )=\phi (t;\mu ,\nu {)}^{\ast }}$. In the tomographic representation, the values ${\displaystyle \phi (1;\mu ,\nu )}$ over phase space also coincide, up to the adopted conventions, with the characteristic function associated with the Wigner representation of the state.

Taking the Fourier transform of the kernel density estimator gives a kernel estimator of the characteristic function,

${\displaystyle {\hat{\phi }}_{n,h}(t;\mu ,\nu )={\hat{\phi }}_{X,n}(t;\mu ,\nu ){\phi }_K(th;\mu ,\nu ),}$

where

${\displaystyle {\hat{\phi }}_{X,n}(t;\mu ,\nu )=\frac{1}{n}\sum _{j=1}^n{e}^{it{X}_{j,\mu ,\nu }}}$

is the empirical characteristic function and

${\displaystyle {\phi }_K(th;\mu ,\nu )={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{K}_{\mu ,\nu }(z)e^{ithz}dz}$

is the known characteristic function of the rescaled kernel. Thus, the estimator separates into a data-dependent empirical factor and a deterministic smoothing factor controlled by the kernel and bandwidth. The kernel factor suppresses high-frequency contributions, providing a regularized version of the empirical characteristic function.

The characteristic-function representation can be used directly for quantum-state reconstruction. For example, the density matrix in the coordinate basis can be recovered from

${\displaystyle \rho (y,y^{\prime })=\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\phi (1;\mu ,y-y^{\prime })\exp \left[-\frac{i\mu (y+y^{\prime })}{2}\right]d\mu .}$

It can also be used to calculate trace quantities without first reconstructing a phase-space quasiprobability distribution. In particular, two density operators satisfy the Parseval-type relation

${\displaystyle \mathrm{tr}({\rho }_1{\rho }_2)=\frac{1}{2\pi }{\iint }_{{\mathbb{R}}^2}{\phi }_1(1;\mu ,\nu ){\phi }_2(1;\mu ,\nu {)}^{\ast }d\mu d\nu .}$

When the two states are identical, this expression gives the purity ${\displaystyle \mathrm{tr}({\rho }^2)}$; for different states it gives their Hilbert–Schmidt overlap. Related integral expressions can be constructed for higher-order traces, state distances, and phase-space quasiprobability distributions.

The characteristic-function formulation also permits correction for some measurement-noise models. If the observed quadrature is modeled as ${\displaystyle Z=\kappa X+(1-\kappa )Y}$, where ${\displaystyle Y}$ is independent noise with a known nonvanishing characteristic function, then the signal characteristic function can be obtained through a deconvolution,

${\displaystyle {\phi }_X(t)=\frac{{\phi }_Z(t/\kappa )}{{\phi }_Y((1-\kappa )t/\kappa )}.}$

A kernel factor can subsequently be included to regularize the resulting estimate. This correction requires the noise characteristic function not to vanish at the frequencies used in the reconstruction.

The bias and variance of the KCFE are

${\displaystyle \mathrm{Bias}\left[{\hat{\phi }}_{n,h}(t;\mu ,\nu )\right]=\left[{\phi }_K(th;\mu ,\nu )-1\right]\phi (t;\mu ,\nu )}$

and

${\displaystyle \mathrm{Var}\left[{\hat{\phi }}_{n,h}(t;\mu ,\nu )\right]=\frac{|{\phi }_K(th;\mu ,\nu ){|}^2}{n}\left[1-|\phi (t;\mu ,\nu ){|}^2\right].}$

Thus its pointwise mean squared error is

${\displaystyle \begin{array}{cc}\mathrm{MSE}\left[{\hat{\phi }}_{n,h}(t;\mu ,\nu )\right]=& |{\phi }_K(th;\mu ,\nu )-1{|}^2|\phi (t;\mu ,\nu ){|}^2\\ & +\frac{|{\phi }_K(th;\mu ,\nu ){|}^2}{n}\left[1-|\phi (t;\mu ,\nu ){|}^2\right].\end{array}}$

For the Gaussian kernel

${\displaystyle {\phi }_K(th;\mu ,\nu )=\exp \left[-\frac{t^2{h}^2({\mu }^2+{\nu }^2)}{4}\right].}$

The resulting error is

${\displaystyle \begin{array}{cc}\mathrm{MSE}\left[{\hat{\phi }}_{n,h}^G(t;\mu ,\nu )\right]=& {\left(\exp \left[-\frac{t^2{h}^2({\mu }^2+{\nu }^2)}{4}\right]-1\right)}^2|\phi (t;\mu ,\nu ){|}^2\\ & +\frac{1}{n}\exp \left[-\frac{t^2{h}^2({\mu }^2+{\nu }^2)}{4}\right]\left[1-|\phi (t;\mu ,\nu ){|}^2\right].\end{array}}$

For ${\displaystyle 0<|\phi (t;\mu ,\nu )|<1}$, minimizing this expression gives a characteristic-function bandwidth of order

${\displaystyle h_{\phi }=O(n^{-1/2}),}$

and the minimized error satisfies

${\displaystyle \mathrm{MSE}\left[{\hat{\phi }}_{n,h_{\phi }}^G(t;\mu ,\nu )\right]=\frac{1-|\phi (t;\mu ,\nu ){|}^2}{n}-\frac{{\left[1-|\phi (t;\mu ,\nu ){|}^2\right]}^2}{4n^2|\phi (t;\mu ,\nu ){|}^2}=O(n^{-1}).}$

The characteristic-function estimator therefore converges pointwise at the parametric ${\displaystyle n^{-1}}$ mean squared error rate, although the underlying tomogram is estimated nonparametrically.

Collections of characteristic-function estimates at several measurement settings form the input to the kernel quantum state estimation procedure. If ${\displaystyle n}$ measurements are taken at each setting and ${\displaystyle N_{\mu }=N_{\nu }=(\log n{)}^2}$, the numbers of state preparations used for density-matrix and two-dimensional trace estimation are

${\displaystyle T_{\mu }=n{N}_{\mu },T_{\mu ,\nu }=n{N}_{\mu }{N}_{\nu }.}$

For the coordinate-space density matrix, the reported uniform pointwise risk is

${\displaystyle \underset{y,y^{\prime }}{sup}\mathbb{E}\left[{\left|\hat{\rho }(y,y^{\prime })-\rho (y,y^{\prime })\right|}^2\right]=\tilde{O}\left(T_{\mu }^{-1}\right).}$

For the overlap of two quantum states, the reported mean squared error is

${\displaystyle \mathbb{E}\left[{\left|\widehat{\mathrm{tr}({\rho }_1{\rho }_2)}-\mathrm{tr}({\rho }_1{\rho }_2)\right|}^2\right]=\tilde{O}\left(T_{\mu ,\nu }^{-1}\right),}$

where ${\displaystyle \tilde{O}}$ suppresses logarithmic factors arising from the number of measurement settings.

Geometric and topological features

We can extend the definition of the (global) mode to a local sense and define the local modes:

$${\displaystyle M=\{x:g(x)=0\mid {\lambda }_1(x)<0\}}$$

Namely, ${\displaystyle M}$ is the collection of points for which the density function is locally maximized. A natural estimator of ${\displaystyle M}$ is a plug-in from KDE, where ${\displaystyle g(x)}$ and ${\displaystyle {\lambda }_1(x)}$ are KDE version of ${\displaystyle g(x)}$ and ${\displaystyle {\lambda }_1(x)}$. Under mild assumptions, ${\displaystyle M_c}$ is a consistent estimator of ${\displaystyle M}$. Note that one can use the mean shift algorithm to compute the estimator ${\displaystyle M_c}$ numerically.

Statistical implementation

A non-exhaustive list of software implementations of kernel density estimators includes:

• In Analytica release 4.4, the Smoothing option for PDF results uses KDE, and from expressions it is available via the built-in Pdf function.
• In C/C++, FIGTree is a library that can be used to compute kernel density estimates using normal kernels. MATLAB interface available.
• In C++, libagf is a library for variable kernel density estimation.
• In C++, mlpack is a library that can compute KDE using many different kernels. It allows to set an error tolerance for faster computation. Python and R interfaces are available.
• in C# and F#, Math.NET Numerics is an open source library for numerical computation which includes kernel density estimation
• In CrimeStat, kernel density estimation is implemented using five different kernel functions – normal, uniform, quartic, negative exponential, and triangular. Both single- and dual-kernel density estimate routines are available. Kernel density estimation is also used in interpolating a Head Bang routine, in estimating a two-dimensional Journey-to-crime density function, and in estimating a three-dimensional Bayesian Journey-to-crime estimate.
• In ELKI, kernel density functions can be found in the package de.lmu.ifi.dbs.elki.math.statistics.kernelfunctions
• In ESRI products, kernel density mapping is managed out of the Spatial Analyst toolbox and uses the Quartic(biweight) kernel.
• In Excel, the Royal Society of Chemistry has created an add-in to run kernel density estimation based on their Analytical Methods Committee Technical Brief 4.
• In gnuplot, kernel density estimation is implemented by the smooth kdensity option, the datafile can contain a weight and bandwidth for each point, or the bandwidth can be set automatically according to "Silverman's rule of thumb" (see above).
• In Haskell, kernel density is implemented in the statistics package.
• In IGOR Pro, kernel density estimation is implemented by the StatsKDE operation (added in Igor Pro 7.00). Bandwidth can be user specified or estimated by means of Silverman, Scott or Bowmann and Azzalini. Kernel types are: Epanechnikov, Bi-weight, Tri-weight, Triangular, Gaussian and Rectangular.
• In Java, the Weka machine learning package provides weka.estimators.KernelEstimator, among others.
• In JavaScript, the visualization package D3.js offers a KDE package in its science.stats package.
• In JMP, the Graph Builder platform utilizes kernel density estimation to provide contour plots and high density regions (HDRs) for bivariate densities, and violin plots and HDRs for univariate densities. Sliders allow the user to vary the bandwidth. Bivariate and univariate kernel density estimates are also provided by the Fit Y by X and Distribution platforms, respectively.
• In Julia, kernel density estimation is implemented in the KernelDensity.jl package.
• In KNIME, 1D and 2D Kernel Density distributions can be generated and plotted using nodes from the Vernalis community contribution, e.g. 1D Kernel Density Plot, among others. The underlying implementation is written in Java.
• In MATLAB, kernel density estimation is implemented through the ksdensity function (Statistics Toolbox). As of the 2018a release of MATLAB, both the bandwidth and kernel smoother can be specified, including other options such as specifying the range of the kernel density. Alternatively, a free MATLAB software package which implements an automatic bandwidth selection method is available from the MATLAB Central File Exchange for
  • 1-dimensional data
  • 2-dimensional data
  • n-dimensional data Archived 2017-02-09 at the Wayback Machine
    A free MATLAB toolbox with implementation of kernel regression, kernel density estimation, kernel estimation of hazard function and many others is available on these pages (this toolbox is a part of the book ).
• In Mathematica, numeric kernel density estimation is implemented by the function SmoothKernelDistribution and symbolic estimation is implemented using the function KernelMixtureDistribution both of which provide data-driven bandwidths.
• In Minitab, the Royal Society of Chemistry has created a macro to run kernel density estimation based on their Analytical Methods Committee Technical Brief 4.
• In the NAG Library, kernel density estimation is implemented via the g10ba routine (available in both the Fortran and the C versions of the Library).
• In Nuklei, C++ kernel density methods focus on data from the Special Euclidean group ${\displaystyle SE(3)}$.
• In Octave, kernel density estimation is implemented by the kernel_density option (econometrics package).
• In Origin, 2D kernel density plot can be made from its user interface, and two functions, Ksdensity for 1D and Ks2density for 2D can be used from its LabTalk Archived 2014-10-19 at the Wayback Machine, Python, or C code.
• In Perl, an implementation can be found in the Statistics-KernelEstimation module
• In PHP, an implementation can be found in the MathPHP library
• In Python, many implementations exist: pyqt_fit.kde Module in the PyQt-Fit package, SciPy (scipy.stats.gaussian_kde), Statsmodels (KDEUnivariate and KDEMultivariate), and scikit-learn (KernelDensity) (see comparison). KDEpy supports weighted data and its FFT implementation is orders of magnitude faster than the other implementations. The commonly used pandas library offers support for kde plotting through the plot method (df.plot(kind='kde')). The getdist package for weighted and correlated MCMC samples supports optimized bandwidth, boundary correction and higher-order methods for 1D and 2D distributions. One newly used package for kernel density estimation is seaborn ( import seaborn as sns , sns.kdeplot() ). A GPU implementation of KDE also exists.
• In R, it is implemented through density in the base distribution, and bw.nrd0 function is used in stats package, this function uses the optimized formula in Silverman's book. bkde in the KernSmooth library, ParetoDensityEstimation in the DataVisualizations library (for pareto distribution density estimation), kde in the ks library, dkden and dbckden in the evmix library (latter for boundary corrected kernel density estimation for bounded support), npudens in the np library (numeric and categorical data), sm.density in the sm library. For an implementation of the kde.R function, which does not require installing any packages or libraries, see kde.R. The btb library, dedicated to urban analysis, implements kernel density estimation through kernel_smoothing.
• In SAS, proc kde can be used to estimate univariate and bivariate kernel densities.
• In Apache Spark, the KernelDensity() class
• In Stata, it is implemented through kdensity; for example histogram x, kdensity. Alternatively a free Stata module KDENS is available allowing a user to estimate 1D or 2D density functions.
• In Swift, it is implemented through SwiftStats.KernelDensityEstimation in the open-source statistics library SwiftStats.