Multivariate kernel density estimation: Difference between revisions

Revision as of 09:07, 14 September 2023

Latest revision as of 12:02, 17 June 2025

(9 intermediate revisions by 7 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 generalizationgeneralisation 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=978-0-387-94716-79780387947167}}</ref><ref name="chacon2018">{{Cite book| author=Chacón, J.E. and Duong, T. | title=Multivariate Kernel Smoothing and Its applications | publisher=Chapman & Hall/CRC | year=2018 | isbn=9781498763011}}</ref>

==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,&nbsp;−1.5): for the one on the right, we shift the anchor point by 0.125 in both directions to (−1.625,&nbsp;−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=978-0-412-24620-39780412246203 | pages=[https://archive.org/details/densityestimatio00silv_0/page/7 7–11] | url-access=registration | url=https://archive.org/details/densityestimatio00silv_0/page/7 }}</ref>

[[File:Synthetic data 2D histograms.png|thumb|center|500px|alt=Left. Histogram with anchor point at (−1.5,&nbsp;-1.5). Right. Histogram with anchor point at (−1.625,&nbsp;−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,&nbsp;-1.5). Right. Histogram with anchor point at (−1.625,&nbsp;−1.625). Both histograms have a bindwidthbinwidth of 0.5, so differences in appearances of the two histograms are due to the placement of the anchor point.]]

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">{{Cite book| author1=Wand, M.P | author2=Jones, M.C. | title=Kernel Smoothing | publisher=Chapman & Hall/CRC | ___location=London | year=1995 | isbn = 978-0-412-55270-0}}</ref>{{rp|36–39}} That the bandwidth matrix also induces an orientation is a basic difference between multivariate kernel density estimation from its univariate analogue since orientation is not defined for 1D kernels. This leads to the choice of the parameterizationparametrisation of this bandwidth matrix. The three main parameterizationparametrisation classes (in increasing order of complexity) are ''S'', the class of positive scalars times the identity matrix; ''D'', diagonal matrices with positive entries on the main diagonal; and ''F'', symmetric positive definite matrices. The ''S'' class kernels have the same amount of smoothing applied in all coordinate directions, ''D'' kernels allow different amounts of smoothing in each of the coordinates, and ''F'' kernels allow arbitrary amounts and orientation of the smoothing. Historically ''S'' and ''D'' kernels are the most widespread due to computational reasons, but research indicates that important gains in accuracy can be obtained using the more general ''F'' class kernels.<ref>{{cite journal | author1=Wand, M.P. | author2=Jones, M.C. | title=Comparison of smoothing parameterizations in bivariate kernel density estimation | journal=Journal of the American Statistical Association | year=1993 | volume=88 | issue=422 | pages=520–528 | doi=10.1080/01621459.1993.10476303 | jstor=2290332}}</ref><ref name="DH2003">{{Cite journal| doi=10.1080/10485250306039 | author1=Duong, T. | author2=Hazelton, M.L. | title=Plug-in bandwidth matrices for bivariate kernel density estimation | journal=Journal of Nonparametric Statistics | year=2003 | volume=15 | pages=17–30}}</ref>

[[File:Kernel parameterizationparametrisation class.png|thumb|center|500px|alt=Comparison of the three main bandwidth matrix parameterizationparametrisation 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 parameterizationparametrisation 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.]]

==Optimal bandwidth matrix selection==

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''^−4/(d+4)), hence the covergenceconvergence rate of the density estimator to ''f'' is ''O_p''(n^{−2/(''d''+4)}) where ''O_p'' denotes [[Big O in probability notation|order in probability]]. This establishes pointwise convergence. The functional covergenceconvergence is established similarly by considering the behaviour of the MISE, and noting that under sufficient regularity, integration does not affect the convergence rates.

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_p''(''n''^−''α''), ''α'' > 0 if

Line 103 ⟶ 104:

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_p''(''n''^{−2/(''d''+6)}) <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.

fix syntaxhighlight error==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.]]

Line 117 ⟶ 118:

H <- Hpi(x=faithful)

fhat <- kde(x=faithful, H=H)

plot(fhat, display="filled.contour2contour", drawpoints=TRUE, cex=0.5, pch=16, col.pt=1)

points(faithful, cex=0.5, pch=16)

</syntaxhighlight>

Line 154:

<syntaxhighlight lang="matlab" style="overflow:auto;">

clear all

% generate synthetic data

data=[randn(500, 2);

randn(500, 1) + 3.5, randn(500, 1);];

% call the routine, which has been saved in the current directory

[bandwidth, density, X, Y] = kde2d(data);

% plot the data and the density estimate

contour3(X, Y, density, 50), hold on

plot(data(:,1), data(:,2), 'r.', 'MarkerSize', 5)

</syntaxhighlight>

Line 188:

== 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- colored, monotone areas away from the first contiguous hypervolume correspond to additional contiguous hypervolumes (areas) with <math>|\hat{\varphi}|^2 \ge 4(N-1)N^{-2}</math>. The colors of these areas are arbitrary and only serve to visually differentiate nearby contiguous areas from one another.]]

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]]):

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://bitbucketgithub.orgcom/lblLBL-cascadeEESA/fastkde fastKDE]'' in the literature.<ref name=":22"/>

[[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/tduong/researchmvksa mvstat.net]<em>Multivariate AKernel collectionSmoothing ofand peer-reviewedIts articlesApplications</em>] ofis thea mathematicalcomprehensive detailsbook ofon many topics multivariatein kernel smoothing, including density estimation. andIncludes their[https://cran.r-project.org/web/packages/ks/index.html bandwidthks selectorspackage] oncode ansnippets {{mono|mvstat.net}}in web[[R_programming pagelanguage|R]].

* [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]].