Content deleted Content added
| (21 intermediate revisions by 14 users not shown) | |||
Line 1:
{{Short description|Concept in statistics mathematics}}
[[Kernel density estimation]] is a [[nonparametric]] technique for [[density estimation]] i.e., estimation of [[probability density function]]s, which is one of the fundamental questions in [[statistics]]. It can be viewed as a generalisation of [[histogram]] density estimation with improved statistical properties. Apart from histograms, other types of density estimators include [[parametric statistics|parametric]], [[spline interpolation|spline]], [[wavelet]] and [[Fourier series]]. Kernel density estimators were first introduced in the scientific literature for [[univariate]] data in the 1950s and 1960s<ref>{{Cite journal| doi=10.1214/aoms/1177728190 | last=Rosenblatt | first=M.| title=Remarks on some nonparametric estimates of a density function | journal=Annals of Mathematical Statistics | year=1956 | volume=27 | issue=3 | pages=832–837| doi-access=free }}</ref><ref>{{Cite journal| doi=10.1214/aoms/1177704472| last=Parzen | first=E.| title=On estimation of a probability density function and mode | journal=Annals of Mathematical Statistics| year=1962 | volume=33 | issue=3 | pages=1065–1076| doi-access=free }}</ref> and subsequently have been widely adopted. It was soon recognised that analogous estimators for multivariate data would be an important addition to [[multivariate statistics]]. Based on research carried out in the 1990s and 2000s, '''multivariate kernel density estimation''' has reached a level of maturity comparable to its univariate counterparts.<ref name="WJ1995">{{Cite book| author1=Wand, M.P | author2=Jones, M.C. | title=Kernel Smoothing | publisher=Chapman & Hall/CRC | ___location=London | year=1995 | isbn = 9780412552700}}</ref><ref name="simonoff1996">{{Cite book| author=Simonoff, J.S. | title=Smoothing Methods in Statistics | publisher=Springer | year=1996 | isbn=
==Motivation==
We take an illustrative [[Synthetic data|synthetic]] [[bivariate data|bivariate]] data set of 50 points to illustrate the construction of histograms. This requires the choice of an anchor point (the lower left corner of the histogram grid). For the histogram on the left, we choose (−1.5, −1.5): for the one on the right, we shift the anchor point by 0.125 in both directions to (−1.625, −1.625). Both histograms have a binwidth of 0.5, so any differences are due to the change in the anchor point only. The colour-coding indicates the number of data points which fall into a bin: 0=white, 1=pale yellow, 2=bright yellow, 3=orange, 4=red. The left histogram appears to indicate that the upper half has a higher density than the lower half, whereas the reverse is the case for the right-hand histogram, confirming that histograms are highly sensitive to the placement of the anchor point.<ref>{{Cite book | author=Silverman, B.W. | title=Density Estimation for Statistics and Data Analysis | publisher=Chapman & Hall/CRC | year=1986 | isbn=
[[File:Synthetic data 2D histograms.png|thumb|center|500px|alt=Left. Histogram with anchor point at (−1.5, -1.5). Right. Histogram with anchor point at (−1.625, −1.625). Both histograms have a bin width of 0.5, so differences in appearances of the two histograms are due to the placement of the anchor point.|Comparison of 2D histograms. Left. Histogram with anchor point at (−1.5, -1.5). Right. Histogram with anchor point at (−1.625, −1.625). Both histograms have a
One possible solution to this anchor point placement problem is to remove the histogram binning grid completely. In the left figure below, a kernel (represented by the grey lines) is centred at each of the 50 data points above. The result of summing these kernels is given on the right figure, which is a kernel density estimate. The most striking difference between kernel density estimates and histograms is that the former are easier to interpret since they do not contain artifices induced by a binning grid.
Line 24 ⟶ 25:
* <math>K_\mathbf{H}(\mathbf{x})=|\mathbf{H}|^{-1/2}K(\mathbf{H}^{-1/2}\mathbf{x} )</math>.
The choice of the kernel function ''K'' is not crucial to the accuracy of kernel density estimators, so we use the standard [[multivariate normal distribution|multivariate normal]] kernel throughout: <math display="inline">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} }</math>, where H plays the role of the [[covariance matrix]]. On the other hand, the choice of the bandwidth matrix <strong>H</strong> is the single most important factor affecting its accuracy since it controls the amount and orientation of smoothing induced.<ref name="WJ1995"
[[File:Kernel parametrisation class.png|thumb|center|500px|alt=Comparison of the three main bandwidth matrix parametrisation classes. Left. S positive scalar times the identity matrix. Centre. D diagonal matrix with positive entries on the main diagonal. Right. F symmetric positive definite matrix.|Comparison of the three main bandwidth matrix parametrisation classes. Left. ''S'' positive scalar times the identity matrix. Centre. ''D'' diagonal matrix with positive entries on the main diagonal. Right. ''F'' symmetric positive definite matrix.]]
Line 74 ⟶ 75:
+ K_{2\mathbf{G}}) (\mathbf{X}_i - \mathbf{X}_j)</math>
Thus <math>\hat{\mathbf{H}}_{\operatorname{SCV}} = \operatorname{argmin}_{\mathbf{H} \in F} \, \operatorname{SCV} (\mathbf{H})</math> is the SCV selector.<ref name="DH2005">{{Cite journal| doi=10.1111/j.1467-9469.2005.00445.x | author1=Duong, T. | author2=Hazelton, M.L. | title=Cross validation bandwidth matrices for multivariate kernel density estimation | journal=Scandinavian Journal of Statistics | year=2005 | volume=32 | issue=3 | pages=485–506}}</ref><ref>{{Cite journal| doi=10.1007/BF01205233 | author1=Hall, P. | author2=Marron, J. | author3=Park, B. | title=Smoothed cross-validation | journal=Probability Theory and Related Fields | year=1992 | volume=92 | pages=1–20| doi-access=free }}</ref>
These references also contain algorithms on optimal estimation of the pilot bandwidth matrix <strong>G</strong> and establish that <math>\hat{\mathbf{H}}_{\operatorname{SCV}}</math> converges in probability to '''H'''<sub>AMISE</sub>.
=== Rule of thumb ===
Silverman's rule of thumb suggests using <math>\sqrt{\mathbf{H}_{ii}} = \left(\frac{4}{d+2}\right)^{\frac{1}{d+4}} n^{\frac{-1}{d+4}} \sigma_i</math>, where <math>\sigma_i</math> is the standard deviation of the ith variable and <math>d</math> is the number of dimensions, and <math>\mathbf{H}_{ij} = 0, i\neq j</math>. Scott's rule is <math>\sqrt{\mathbf{H}_{ii}} = n^{\frac{-1}{d+4}} \sigma_i</math>.
==Asymptotic analysis==
Line 89 ⟶ 90:
where * is the [[convolution]] operator between two functions, and
:<math>\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}).</math>
For these two expressions to be well-defined, we require that all elements of '''H''' tend to 0 and that ''n''<sup>−1</sup> |'''H'''|<sup>−1/2</sup> tends to 0 as ''n'' tends to infinity. Assuming these two conditions, we see that the expected value tends to the true density ''f'' i.e. the kernel density estimator is asymptotically [[Bias of an estimator|unbiased]]; and that the variance tends to zero. Using the standard mean squared value decomposition
Line 95 ⟶ 96:
:<math>\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</math>
we have that the MSE tends to 0, implying that the kernel density estimator is (mean square) consistent and hence converges in probability to the true density ''f''. The rate of convergence of the MSE to 0 is the necessarily the same as the MISE rate noted previously ''O''(''n''<sup>−4/(d+4)</sup>), hence the
For the data-based bandwidth selectors considered, the target is the AMISE bandwidth matrix. We say that a data-based selector converges to the AMISE selector at relative rate ''O<sub>p</sub>''(''n''<sup>−''α''</sup>), ''α'' > 0 if
Line 101 ⟶ 102:
:<math>\operatorname{vec} (\hat{\mathbf{H}} - \mathbf{H}_{\operatorname{AMISE}}) = O(n^{-2\alpha}) \operatorname{vec} \mathbf{H}_{\operatorname{AMISE}}.</math>
It has been established that the plug-in and smoothed cross validation selectors (given a single pilot bandwidth '''G''') both converge at a relative rate of ''O<sub>p</sub>''(''n''<sup>−2/(''d''+6)</sup>) <ref name="DH2005" /><ref>{{Cite journal| doi=10.1016/j.jmva.2004.04.004 | author1=Duong, T. | author2=Hazelton, M.L. | title=Convergence rates for unconstrained bandwidth matrix selectors in multivariate kernel density estimation | journal=Journal of Multivariate Analysis | year=2005 | volume=93 | issue=2 | pages=417–433| doi-access=free }}</ref> i.e., both these data-based selectors are consistent estimators.
==Density estimation with a full bandwidth matrix==
[[File:Old Faithful Geyser KDE with plugin bandwidth.png|thumb|250px|alt=Old Faithful Geyser data kernel density estimate with plug-in bandwidth matrix.|Old Faithful Geyser data kernel density estimate with plug-in bandwidth matrix.]]
The [https://cran.r-project.org/web/packages/ks/index.html ks package]<ref>{{Cite journal| author1=Duong, T. | title=ks: Kernel density estimation and kernel discriminant analysis in R | journal=Journal of Statistical Software | year=2007 | volume=21 | issue = 7 | doi=10.18637/jss.v021.i07 | doi-access=free }}</ref> in [[R programming language|R]] implements the plug-in and smoothed cross validation selectors (amongst others). This dataset (included in the base distribution of R) contains
272 records with two measurements each: the duration time of an eruption (minutes) and the
waiting time until the next eruption (minutes) of the [[Old Faithful Geyser]] in Yellowstone National Park, USA.
Line 112 ⟶ 113:
The code fragment computes the kernel density estimate with the plug-in bandwidth matrix <math>\hat{\mathbf{H}}_{\operatorname{PI}} = \begin{bmatrix}0.052 & 0.510 \\ 0.510 & 8.882\end{bmatrix}.</math> Again, the coloured contours correspond to the smallest region which contains the respective probability mass: red = 25%, orange + red = 50%, yellow + orange + red = 75%. To compute the SCV selector, <code>Hpi</code> is replaced with <code>Hscv</code>. This is not displayed here since it is mostly similar to the plug-in estimate for this example.
<
library(ks)
data(faithful)
H <- Hpi(x=faithful)
fhat <- kde(x=faithful, H=H)
plot(fhat, display="filled.
</syntaxhighlight>
==Density estimation with a diagonal bandwidth matrix==
Line 153:
in the current directory.
<
</syntaxhighlight>
==Alternative optimality criteria==
Line 171:
: <math>\operatorname{MIAE} (\mathbf{H}) = \operatorname{E}\, \int |\hat{f}_\mathbf{H} (\mathbf{x}) - f(\mathbf{x})| \, d\mathbf{x}.</math>
Its mathematical analysis is considerably more difficult than the MISE ones. In practice, the gain appears not to be significant.<ref>{{cite journal | author1=Hall, P. | author2=Wand, M.P. | title=Minimizing L<sub>1</sub> distance in nonparametric density estimation | journal = Journal of Multivariate Analysis | year=1988 | volume=26 | pages=59–88 | doi=10.1016/0047-259X(88)90073-5| doi-access=free }}</ref> The ''L<sub>∞</sub>'' norm is the Mean Uniform Absolute Error
: <math>\operatorname{MUAE} (\mathbf{H}) = \operatorname{E}\, \operatorname{sup}_{\mathbf{x}} |\hat{f}_\mathbf{H} (\mathbf{x}) - f(\mathbf{x})|.</math>
Line 183:
: <math>\operatorname{MH} (\mathbf{H}) = \operatorname{E} \int (\hat{f}_\mathbf{H} (\mathbf{x})^{1/2} - f(\mathbf{x})^{1/2})^2 \, d\mathbf{x} .</math>
The KL can be estimated using a cross-validation method, although KL cross-validation selectors can be sub-optimal even if it remains [[Consistent estimator|consistent]] for bounded density functions.<ref>{{cite journal | author=Hall, P. | title=On Kullback-Leibler loss and density estimation | journal=Annals of Statistics | volume=15 | issue=4 | year=1989 | pages=589–605 | doi=10.1214/aos/1176350606| doi-access=free }}</ref> MH selectors have been briefly examined in the literature.<ref>{{cite journal | author1=Ahmad, I.A. | author2=Mugdadi, A.R. | title=Weighted Hellinger distance as an error criterion for bandwidth selection in kernel estimation | journal=Journal of Nonparametric Statistics | volume=18 | issue=2 | year=2006 | pages=215–226 | doi=10.1080/10485250600712008}}</ref>
All these optimality criteria are distance based measures, and do not always correspond to more intuitive notions of closeness, so more visual criteria have been developed in response to this concern.<ref>{{cite journal | author1=Marron, J.S. | author2=Tsybakov, A. | title=Visual error criteria for qualitative smoothing | journal = Journal of the American Statistical Association | year=1996 | volume=90 | issue=430 | pages=499–507 | doi=10.2307/2291060 | jstor=2291060}}</ref>
== Objective and data-driven kernel selection ==
[[File:Empirical Characteristic Function.jpg|alt=An x-shaped region of empirical characteristic function in Fourier space.|thumb|Demonstration of the filter function <math>I_{\vec{A}}(\vec{t})</math>. The square of the empirical distribution function <math>|\hat{\varphi}|^2</math> from ''N''=10,000 samples of the ‘transition distribution’ discussed in Section 3.2 (and shown in Fig. 4), for <math>|\hat{\varphi}|^2 \ge 4(N-1)N^{-2}</math>. There are two color schemes present in this figure. The predominantly dark, multicolored colored ‘X-shaped’ region in the center corresponds to values of <math>|\hat{\varphi}|^2</math> for the lowest contiguous hypervolume (the area containing the origin); the colorbar at right applies to colors in this region. The lightly
Recent research has shown that the kernel and its bandwidth can both be optimally and objectively chosen from the input data itself without making any assumptions about the form of the distribution.<ref name=":0">{{Cite journal|last = Bernacchia|first = Alberto|last2 = Pigolotti|first2 = Simone|date = 2011-06-01|title = Self-consistent method for density estimation|journal = Journal of the Royal Statistical Society, Series B|language = en|volume = 73|issue = 3|pages = 407–422|doi = 10.1111/j.1467-9868.2011.00772.x|issn = 1467-9868|arxiv = 0908.3856}}</ref> The resulting kernel density estimate converges rapidly to the true probability distribution as samples are added: at a rate close to the <math>n^{-1}</math> expected for parametric estimators.<ref name=":0" /><ref name=":1">{{Cite journal|last = O’Brien|first = Travis A.|last2 = Collins|first2 = William D.|last3 = Rauscher|first3 = Sara A.|last4 = Ringler|first4 = Todd D.|date = 2014-11-01|title = Reducing the computational cost of the ECF using a nuFFT: A fast and objective probability density estimation method|journal = Computational Statistics & Data Analysis|volume = 79|pages = 222–234|doi = 10.1016/j.csda.2014.06.002|doi-access = free}}</ref><ref name=":22">{{Cite journal|last = O’Brien|first = Travis A.|last2 = Kashinath|first2 = Karthik|last3 = Cavanaugh|first3 = Nicholas R.|last4 = Collins|first4 = William D.|last5 = O’Brien|first5 = John P.|title = A fast and objective multidimensional kernel density estimation method: fastKDE|journal = Computational Statistics & Data Analysis|volume = 101|pages = 148–160|doi = 10.1016/j.csda.2016.02.014|year = 2016|url = https://escholarship.org/content/qt9g56181p/qt9g56181p.pdf?t=p7qvyp|doi-access = free}}</ref> This kernel estimator works for univariate and multivariate samples alike. The optimal kernel is defined in Fourier space—as the optimal damping function <math>\hat{\psi_h}(\vec{t})</math> (the Fourier transform of the kernel <math>\hat{K}(\vec{x})</math> )-- in terms of the Fourier transform of the data <math>\hat{\varphi}(\vec{t})</math>, the ''[[Characteristic function (probability theory)|empirical characteristic function]]'' (see [[Kernel density estimation]]):
<math>\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]</math> <ref name=":22"/>
Line 197:
where, ''N'' is the number of data points, ''d'' is the number of dimensions (variables), and <math>I_{\vec{A}}(\vec{t})</math> is a filter that is equal to 1 for 'accepted frequencies' and 0 otherwise. There are various ways to define this filter function, and a simple one that works for univariate or multivariate samples is called the 'lowest contiguous hypervolume filter'; <math>I_{\vec{A}}(\vec{t})</math> is chosen such that the only accepted frequencies are a contiguous subset of frequencies surrounding the origin for which <math>|\hat{\varphi}(\vec{t})|^2 \ge 4(N-1)N^{-2}</math> (see <ref name=":22"/> for a discussion of this and other filter functions).
Note that direct calculation of the ''empirical characteristic function'' (ECF) is slow, since it essentially involves a direct Fourier transform of the data samples. However, it has been found that the ECF can be approximated accurately using a [[Non-uniform discrete Fourier transform|non-uniform fast Fourier transform]] (nuFFT) method,<ref name=":1" /><ref name=":22"/> which increases the calculation speed by several orders of magnitude (depending on the dimensionality of the problem). The combination of this objective KDE method and the nuFFT-based ECF approximation has been referred to as ''[https://
[[File:FastKDE_example.jpg|alt=A demonstration of fastKDE relative to a sample PDF. (a) True PDF, (b) a good representation with fastKDE, and (c) a slightly blurry representation.|none|thumb|664x664px|A non-trivial mixture of normal distributions: (a) the underlying PDF, (b) a fastKDE estimate on 1,000,000 samples, and (c) a fastKDE estimate on 10,000 samples.]]
Line 208:
==External links==
* [http://www.mvstat.net/
* [http://www.mathworks.com/matlabcentral/fileexchange/17204-kernel-density-estimation kde2d.m] A [[Matlab]] function for bivariate kernel density estimation.
* [http://libagf.sf.net libagf] A [[C++]] library for multivariate, [[variable bandwidth kernel density estimation]].
| |||