In statistics, kernel density estimation (KDE) is a non-parametric way to estimate the probability density function of a random variable. Kernel density estimation is a fundamental data smoothing problem where inferences about the population are made, based on a finite data sample. In some fields such as signal processing and econometrics it is also termed the Parzen-Rosenblatt window method, after Emanuel Parzen and Murray Rosenblatt, who are usually credited with independently creating it in its current form. One of the famous applications of kernel density estimation is in estimating the class-conditional marginal densities of data when using a naive Bayes classifier, which can improve its prediction accuracy dramatically.
Let (x1, x2, ..., xn) be a univariate independent and identically distributed sample drawn from some distribution with an unknown density ?. We are interested in estimating the shape of this function ?. Its kernel density estimator is
where K is the kernel -- a non-negative function -- and is a smoothing parameter called the bandwidth. A kernel with subscript h is called the scaled kernel and defined as . Intuitively one wants to choose h as small as the data will allow; however, there is always a trade-off between the bias of the estimator and its variance. The choice of bandwidth is discussed in more detail below.
A range of kernel functions are commonly used: uniform, triangular, biweight, triweight, Epanechnikov, normal, and others. The Epanechnikov kernel is optimal in a mean square error sense, though the loss of efficiency is small for the kernels listed previously. Due to its convenient mathematical properties, the normal kernel is often used, which means , where ? is the standard normal density function.
The construction of a kernel density estimate finds interpretations in fields outside of density estimation. For example, in thermodynamics, this is equivalent to the amount of heat generated when heat kernels (the fundamental solution to the heat equation) are placed at each data point locations xi. Similar methods are used to construct discrete Laplace operators on point clouds for manifold learning (e.g. diffusion map).
Kernel density estimates are closely related to histograms, but can be endowed with properties such as smoothness or continuity by using a suitable kernel. To see this, we compare the construction of histogram and kernel density estimators, using these 6 data points:
For the histogram, first the horizontal axis is divided into sub-intervals or bins which cover the range of the data. In this case, we have six bins each of width 2. Whenever a data point falls inside this interval, we place a box of height 1/12. If more than one data point falls inside the same bin, we stack the boxes on top of each other.
For the kernel density estimate, we place a normal kernel with standard deviation 2.25 (indicated by the red dashed lines) on each of the data points xi. The kernels are summed to make the kernel density estimate (solid blue curve). The smoothness of the kernel density estimate is evident compared to the discreteness of the histogram, as kernel density estimates converge faster to the true underlying density for continuous random variables.
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 (no smoothing), where the estimate is a sum of n delta functions centered at the coordinates of analyzed samples. In the other extreme limit the estimate retains the shape of the used kernel, centered on the mean of the samples (completely smooth).
Under weak assumptions on ? and K, (? is the, generally unknown, real density function), MISE (h) = AMISE(h) + o(1/(nh) + h4) where o is the little o notation. The AMISE is the Asymptotic MISE which consists of the two leading terms
where for a function g, and ?'' is the second derivative of ?. The minimum of this AMISE is the solution to this differential equation
Neither the AMISE nor the hAMISE formulas are able to be used directly since they involve the unknown density function ? or its second derivative ?'', so a variety of automatic, data-based methods have been developed for selecting the bandwidth. Many review studies have been carried out 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 hAMISE 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-1 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.
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:
In order to make the h value more robust to make the fitness well for both long-tailed and skew distribution and bimodal mixture distribution, it is better to substitute the value of with another parameter A, which is given by:
Another modification that will improve the model is to reduce the factor from 1.06 to 0.9. Then the final formula would be:
where is the standard deviation of the samples, n is the sample size. IQR is the interquartile range. This approximation is termed the normal distribution approximation, Gaussian approximation, or Silverman's (1986) 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, consider estimating the bimodal Gaussian mixture:
from a sample of 200 points. The figure on the right below 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. The Matlab script for this example uses kde.m and is given below.
1 % Data 2 randn('seed',1) % Used for reproducibility 3 data = [randn(100,1)-10; randn(100,1)+10]; % Two Normals mixed 4 % True 5 phi = @(x) exp(-.5*x.^2)/sqrt(2*pi); % Normal Density 6 tpdf = @(x) phi(x+10)/2+phi(x-10)/2; % True Density 7 % Kernel 8 h = std(data)*(4/3/numel(data))^(1/5); % Bandwidth estimated by Silverman's Rule of Thumb 9 kernel = @(x) mean(phi((x-data)/h)/h); % Kernel Density 10 kpdf = @(x) arrayfun(kernel,x); % Elementwise application 11 % Plot 12 figure(2), clf, hold on 13 x = linspace(-25,+25,1000); % Linear Space 14 plot(x,tpdf(x)) % Plot True Density 15 plot(x,kpdf(x)) % Plot Kernel Density with Silverman's Rule of Thumb 16 kde(data) % Plot Kernel Density with Solve-the-Equation Bandwidth
1 #The same code with R language 2 #` Data 3 set.seed(1) #Used for reproducibility 4 data = c(rnorm(100,-10,1),rnorm(100,10,1)) #Two Normals mixed 5 #` True 6 phi = function(x) exp(-.5*x^2)/sqrt(2*pi) #Normal Density 7 tpdf = function(x) phi(x+10)/2+phi(x-10)/2 #True Density 8 #` Kernel 9 h = sd(data)*(4/3/length(data))^(1/5) #Bandwidth estimated by Silverman's Rule of Thumb 10 Kernel2 = function(x) mean(phi((x-data)/h)/h) #Kernel Density 11 kpdf = function(x) sapply(x,Kernel2) #Elementwise application 12 #` Plot 13 x=seq(-25,25,length=1000) #Linear Space 14 plot(x,tpdf(x),type="l",ylim=c(0,0.23),col="red") #Plot True Density 15 par(new=T) 16 plot(x,kpdf(x),type="l",ylim=c(0,0.23),xlab="",ylab="",axes=F) #Plot Kernel Density with Silverman's Rule of Thumb
Given the sample (x1, x2, ..., xn), it is natural to estimate the characteristic function as
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 is unreliable for large t's. To circumvent this problem, the estimator is multiplied by a damping function , 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 . In particular when h is small, then ?h(t) will be approximately one for a large range of t's, which means that remains practically unaltered in the most important region of t's.
The most common choice for function ? is either the uniform function
where K is the Fourier transform of the damping function ?. Thus the kernel density estimator coincides with the characteristic function density estimator.
We can extend the definition of the (global) mode to a local sense and define the local modes:
Namely, is the collection of points for which the density function is locally maximized. A natural estimator of is a plug-in from KDE, where and are KDE version of and . Under mild assumptions, is a consistent estimator of . Note that one can use the mean shift algorithm to compute the estimator numerically.
A non-exhaustive list of software implementations of kernel density estimators includes:
smooth kdensityoption, 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).
StatsKDEoperation (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.
ksdensityfunction (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
SmoothKernelDistributionhere and symbolic estimation is implemented using the function
KernelMixtureDistributionhere both of which provide data-driven bandwidths.
g10baroutine (available in both the Fortran and the C versions of the Library).
kernel_densityoption (econometrics package).
scipy.signal.parzen), Statsmodels (
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,
densityin the base distribution, and
bw.nrd0function is used in stats package, this function uses the optimized formula in Silverman's book.
bkdein the KernSmooth library,
ParetoDensityEstimationin the AdaptGauss library (for pareto distribution density estimation),
kdein the ks library,
dbckdenin the evmix library (latter for boundary corrected kernel density estimation for bounded support),
npudensin the np library (numeric and categorical data),
sm.densityin the sm library. For an implementation of the
kde.Rfunction, which does not require installing any packages or libraries, see kde.R. The btb library, dedicated to urban analysis, implements kernel density estimation through
proc kdecan be used to estimate univariate and bivariate kernel densities.
KernelDensityclass (see official documentation for more details )
kdensity; for example
histogram x, kdensity. Alternatively a free Stata module KDENS is available from here allowing a user to estimate 1D or 2D density functions.
SwiftStats.KernelDensityEstimationin the open-source statistics library SwiftStats.