Multivariate kernel density estimation: Difference between revisions

Revision as of 00:08, 16 September 2010

Latest revision as of 12:02, 17 June 2025

(155 intermediate revisions by 85 users not shown)

Line 1:

{{Short description|Concept in statistics mathematics}}

{{Userspace draft|source=ArticleWizard|date=September 2010}}

[[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=9780387947167}}</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==

[[Kernel density estimation]] is one of the most popular techniques for [[density estimation]] i.e., estimation of [[probability density function|probability density functions]], which is one of the fundamental questions in [[statistics]].

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=9780412246203 | pages=[https://archive.org/details/densityestimatio00silv_0/page/7 7–11] | url-access=registration | url=https://archive.org/details/densityestimatio00silv_0/page/7 }}</ref>

It can be viewed as a generalisation of [[histogram]] density estimation with improved statistical properties.

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 | pages=832-837}}</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 | pages=1065-1076}}</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 their univariate counterparts.<ref>{{cite book | author=Simonoff, J.S. | title=Smoothing Methods in Statistics | publisher=Springer | date=1996 | isbn=0387947167}}</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 binwidth 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.

== Motivation ==

The coloured contours correspond to the smallest region which contains the respective probability mass: red = 25%, orange + red = 50%, yellow + orange + red = 75%, thus indicating that a single central region contains the highest density.

To motivate the definition of multivariate kernel density estimators, we take as an illustrative synthetic [[bivariate]] data set of 50 points. First we create a histogram. 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 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 it is the reverse is the case. It is well-known that histograms are highly sensitive the placement of the anchor point.

[[File:Synthetic data 2D KDE.png|thumb|center|500px|alt=Left. Individual kernels. Right. Kernel density estimate.|Construction of 2D kernel density estimate. Left. Individual kernels. Right. Kernel density estimate.]]

[[Image:Synthetic data 2D histograms.png|center]]

The goal of density estimation is to take a finite sample of data and to make inferences about the underlying probability density function everywhere, including where no data are observed. In kernel density estimation, the contribution of each data point is smoothed out from a single point into a region of space surrounding it. Aggregating the individually smoothed contributions gives an overall picture of the structure of the data and its density function. In the details to follow, we show that this approach leads to a reasonable estimate of the underlying density function.

One possible solution to this anchor point placement problem to remove the histogram binning grid completely. In the left figure below, a kernel (represented by the dashed grey lines) is placed at the same 50 data points above. Then the result of summing these kernels is given on the right figure, which is a kernel density estimate.

==Definition==

[[Image:Synthetic data 2D KDE.png|center]]

The previous figure is a graphical representation of kernel density estimate, which we now define in an exact manner. Let '''x'''₁, '''x'''₂, ..., '''x'''_''n'' be a [[random sample|sample]] of ''d''-variate [[random vector]]s drawn from a common distribution described by the [[probability density function|density function]] ''ƒ''. The kernel density estimate is defined to be

: <math>

\hat{f}_\mathbf{H}(\mathbf{x})= \frac1n \sum_{i=1}^n K_\mathbf{H} (\mathbf{x} - \mathbf{x}_i)

</math>

where

* {{nowrap|'''x''' {{=}} (''x''₁, ''x''₂, …, ''x_d'')^''T''}}, {{nowrap|'''x'''_''i'' {{=}} (''x''_''i''1, ''x''_''i''2, …, ''x_id'')^''T'', ''i'' {{=}} 1, 2, …, ''n''}} are ''d''-vectors;

* '''H''' is the bandwidth (or smoothing) ''d×d'' matrix which is [[symmetric matrix|symmetric]] and [[positive definite matrix|positive definite]];

* ''K'' is the [[kernel (statistics)|kernel]] function which is a symmetric multivariate density;

* <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"/>{{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 parametrisation of this bandwidth matrix. The three main parametrisation 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>

== Definition ==

The previous figure is a graphical representation of kernel density estimate. We now define it in an exact manner.

Let <math>\bold{X}_1, \bold{X}_2, \dots, \bold{X}_n</math> be a <em>d</em>-variate random sample drawn from a common density function <em>f</em>. The kernel density estimate is defined to be

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

<math>\hat{f}_\bold{H}(\bold{x})= n^{-1} |\bold{H}|^{-1/2} \sum_{i=1}^n K_\bold{H} (\bold{x} - \bold{X}_i)</math>

==Optimal bandwidth matrix selection==

where

The most commonly used optimality criterion for selecting a bandwidth matrix is the MISE or [[mean integrated squared error]]

<ul>

<li><math>\bold{x} = (x_1, x_2, \dots, x_d)^T</math>, <math>\bold{X}_i = (X_{i1}, X_{i2}, \dots, X_{id})^T, i=1, 2, \dots, n</math> are <em>d</em>-vectors

<li><em>K</em> is the kernel function which is a symmetric density function, with <math>K_\bold{H}(\bold{x}) = |\bold{H}|^{-1/2} K(\bold{H}^{-1/2} \bold{x})</math>

<li><strong>H</strong> is the bandwidth (or smoothing) matrix which is a symmetric, [[positive definite matrix|positive definite]] <em>d x d</em> matrix.

</ul>

: <math>\operatorname{MISE} (\mathbf{H}) = \operatorname{E}\!\left[\, \int (\hat{f}_\mathbf{H} (\mathbf{x}) - f(\mathbf{x}))^2 \, d\mathbf{x} \;\right].</math>

The choice of the kernel function <em>K</em> is not crucial to the accuracy of kernel density estimators, so we use the standard [[multivariate normal distribution|multivariate normal]] or Gaussian density function as our kernel <em>K</em> throughout.

Whereas the choice of the bandwidth matrix <strong>H</strong> is the single most important factor affecting its accuracy <ref name="WJ1995">{{cite book | author1=Wand, M.P | author2=Jones, M.C. | title=Kernel Smoothing | publisher=Chapman & Hall | ___location=London | date=1995 | isbn = 0412552701}}</ref>(pp. 36-39). (The bandwith for a kernel density estimate is analogous to the binwidth of histograms).

This in general does not possess a [[closed-form expression]], so it is usual to use its asymptotic approximation (AMISE) as a proxy

<math>K (\bold{x}) = (2\pi)^{-d/2} \exp(-\tfrac{1}{2} \, \bold{x}^T \bold{x}).</math>

== Optimal bandwidth matrix selection ==

The most commonly used optimality criterion for selecting a bandwidth matrix is the MISE or Mean Integrated Squared Error

: <math>\operatorname{MISEAMISE} (\boldmathbf{H}) = E \int [\hatn^{f-1}_ |\boldmathbf{H}|^{-1/2} R(\bold{x}K) -+ f(\boldtfrac{1}{x4} m_2(K)]^2 \, d\bold{x} .</math>

(\operatorname{vec}^T \mathbf{H}) \mathbf{\Psi}_4 (\operatorname{vec} \mathbf{H})</math>

where

This is in general does not possess a closed form expression, so it is usual to use its asymptotic approximation (AMISE) as a proxy

* <math>R(K) = \int K(\mathbf{x})^2 \, d\mathbf{x}</math>, with {{nowrap|''R''(''K'') {{=}} (4''π'')^''−d''/2}} when ''K'' is a normal kernel

* <math>\int \mathbf{x} \mathbf{x}^T K(\mathbf{x}) \, d\mathbf{x} = m_2(K) \mathbf{I}_d</math>,

:with <strong>I</strong>_d being the ''d × d'' [[identity matrix]], with ''m''₂ = 1 for the normal kernel

* D²''ƒ'' is the ''d × d'' Hessian matrix of second order partial derivatives of ''ƒ''

* <math>\mathbf{\Psi}_4 = \int (\operatorname{vec} \, \operatorname{D}^2 f(\mathbf{x})) (\operatorname{vec}^T \operatorname{D}^2 f(\mathbf{x})) \, d\mathbf{x}</math> is a ''d''² × ''d''² matrix of integrated fourth order partial derivatives of ''ƒ''

* vec is the vector operator which stacks the columns of a matrix into a single vector e.g. <math>\operatorname{vec}\begin{bmatrix}a & c \\ b & d\end{bmatrix} = \begin{bmatrix}a & b & c & d\end{bmatrix}^T.</math>

The quality of the AMISE approximation to the MISE<ref name="WJ1995"/>{{rp|97}} is given by

<math>\operatorname{AMISE} (\bold{H}) = n^{-1} |\bold{H}|^{-1/2} R(K) + \tfrac{1}{4} m_2(K)^2

(\operatorname{vec}^T \bold{H}) \bold{\Psi}_4 (\operatorname{vec}^T \bold{H})</math>

: <math>\operatorname{MISE} (\mathbf{H}) = \operatorname{AMISE} (\mathbf{H}) + o(n^{-1} |\mathbf{H}|^{-1/2} + \operatorname{tr} \, \mathbf{H}^2)</math>

where

<ul>

<li><math>R(K) = \int K(\bold{x})^2 \, d\bold{x}</math>, with <math>R(K) = (4 \pi)^{-d/2}</math> when <math>K</math> is a normal kernel

<li><math>\int \bold{x} \bold{x}^T K(\bold{x})^2 \, d\bold{x} = m_2(K) \bold{I}_d</math>,

with <strong>I</strong>_d being the <em>d x d</em> [[identity matrix]], with <em>m</em>₂ = 1 for the normal kernel

<li><math>\operatorname{D}^2 f</math> is the <em>d x d</em> Hessian matrix of second order partial derivatives of <math>f</math>

<li><math>\bold{\Psi}_4 = \int (\operatorname{vec} \, \operatorname{D}^2 f(\bold{x})) (\operatorname{vec}^T \operatorname{D}^2 f(\bold{x})) \, d\bold{x}</math> is a <em>d</em>² x <em>d</em>² matrix of integrated fourth order

partial derivatives of <em>f</em>

<li>vec is the vector operator which stacks the columns of a matrix into a single vector e.g. <math>\operatorname{vec}\begin{bmatrix}a & c \\ b & d\end{bmatrix} = \begin{bmatrix}a & b & c & d\end{bmatrix}^T.</math>

</ul>

This formula of the AMISE is due to <ref name="WJ1995">/</ref>(p. 97). The quality of the AMISE approximation to the MISE is given by

where ''o'' indicates the usual [[big O notation|small o notation]]. Heuristically this statement implies that the AMISE is a 'good' approximation of the MISE as the sample size <var>n</var> → ∞.

<math>\operatorname{MISE} (\bold{H}) = \operatorname{AMISE} (\bold{H}) + o(n^{-1} |\bold{H}|^{-1/2}) + O(\operatorname{tr} \, \bold{H}^2)</math>

It can be shown that any reasonable bandwidth selector '''H''' has '''H''' = ''O''(''n''^{−2/(''d''+4)}) where the [[big O notation]] is applied elementwise. Substituting this into the MISE formula yields that the optimal MISE is ''O''(''n''^{−4/(''d''+4)}).<ref name="WJ1995"/>{{rp|99–100}} Thus as ''n'' → ∞, the MISE → 0, i.e. the kernel density estimate [[convergence in mean|converges in mean square]] and thus also in probability to the true density ''f''. These modes of convergence are confirmation of the statement in the motivation section that kernel methods lead to reasonable density estimators. An ideal optimal bandwidth selector is

where o, O indicate the usual small and [[big O notation]]. Heuristically this statement implies that the AMISE is a 'good' approximation of the MISE as the sample size <em>n → ∞<em>.

: <math>\mathbf{H}_{\operatorname{AMISE}} = \operatorname{argmin}_{\mathbf{H} \in F} \, \operatorname{AMISE} (\mathbf{H}).</math>

The ideal optimal bandwidth selector is

<math>\bold{H}_{\operatorname{AMISE}} = \operatorname{argmin}_{\bold{H} \in F} \, \operatorname{AMISE} (\bold{H})</math>

where <em>F</em> is the space of all symmetric, positive definite matrices.

Since this ideal selector contains the unknown density function <em>f</em>, it cannot be used directly. The many different varieties of data-based bandwidth selectors arise from the different estimators of the AMISE. We concentrate on two classes of selectors which have been shown to be the most widely applicable in practise: smoothed cross validation and plug-in selectors.

Since this ideal selector contains the unknown density function ''ƒ'', it cannot be used directly. The many different varieties of data-based bandwidth selectors arise from the different estimators of the AMISE. We concentrate on two classes of selectors which have been shown to be the most widely applicable in practice: smoothed cross validation and plug-in selectors.

=== Plug-in ===

The plug-in (PI) selector of the AMISE is formed by replacing the Hessian matrix by its estimator

===Plug-in===

<math>\operatorname{PI}(\bold{H}) = n^{-1} |\bold{H}|^{-1/2} R(K) + \tfrac{1}{4} m_2(K)^2

The plug-in (PI) estimate of the AMISE is formed by replacing '''Ψ'''₄ by its estimator <math>\hat{\mathbf{\Psi}}_4</math>

(\operatorname{vec}^T \bold{H}) \hat{\bold{\Psi}}_4 (\bold{G}) (\operatorname{vec} \, \bold{H})</math>

where: <math>\hatoperatorname{\bold{\PsiPI}}_4 (\boldmathbf{GH}) = n^{-1} |\sum_mathbf{i=H}|^{-1/2} R(K) + \tfrac{1}{4} m_2(K)^n2

(\operatorname{vec}^T \mathbf{H}) \hat{\mathbf{\Psi}}_4 (\mathbf{G}) (\operatorname{vec} \, \mathbf{H})</math>

\sum_{j=1}^n [(\operatorname{vec} \, \operatorname{D}^2) (\operatorname{vec}^T \operatorname{D}^2)] K_\bold{G} (\bold{X}_i - \bold{X}_j)</math>. Thus <math>\hat{\bold{H}}_{\operatorname{PI}} = \operatorname{argmin}_{\bold{H} \in F} \, \operatorname{PI} (\bold{H})</math> is the plug-in selector<ref>{{cite journal | author1=Wand, M.P. | author2=Jones, M.C. | title=Multivariate plug-in bandwidth selection | journal=Computational Statistics | year=1994 | volume=9 | pages=97-177}}</ref><ref>{{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>. These references also contain algorithms on optimal estimation of the pilot bandwidth matrix <strong>G</strong>.

where <math>\hat{\mathbf{\Psi}}_4 (\mathbf{G}) = n^{-2} \sum_{i=1}^n

=== Smoothed cross validation ===

\sum_{j=1}^n [(\operatorname{vec} \, \operatorname{D}^2) (\operatorname{vec}^T \operatorname{D}^2)] K_\mathbf{G} (\mathbf{X}_i - \mathbf{X}_j)</math>. Thus <math>\hat{\mathbf{H}}_{\operatorname{PI}} = \operatorname{argmin}_{\mathbf{H} \in F} \, \operatorname{PI} (\mathbf{H})</math> is the plug-in selector.<ref>{{Cite journal| author1=Wand, M.P. | author2=Jones, M.C. | title=Multivariate plug-in bandwidth selection | journal=Computational Statistics | year=1994 | volume=9 | pages=97–177}}</ref><ref name="DH2005" /> These references also contain algorithms on optimal estimation of the pilot bandwidth matrix <strong>G</strong> and establish that <math>\hat{\mathbf{H}}_{\operatorname{PI}}</math> [[convergence in probability|converges in probability]] to '''H'''_AMISE.

Smoothed cross validation (SCV) is a subset of a larger class of [[cross-validation_(statistics) | cross validation]] techniques. The SCV estimator differs from the plug-in estimator in the second term

<math>\operatorname{SCV}(\bold{H}) = n^{-1} |\bold{H}|^{-1/2} R(K) +

n^{-2} \sum_{i=1}^n \sum_{j=1}^n (K_{2\bold{H} +2\bold{G}} - 2K_{\bold{H} +2\bold{G}}

+ K_{2\bold{G}}) (\bold{X}_i - \bold{X}_j)</math>

===Smoothed cross validation===

Thus <math>\hat{\bold{H}}_{\operatorname{SCV}} = \operatorname{argmin}_{\bold{H} \in F} \, \operatorname{SCV} (\bold{H})</math> is the SCV selector<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}}</ref><ref>{{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 | pages=485-506}}</ref>.

Smoothed cross validation (SCV) is a subset of a larger class of [[cross-validation (statistics)|cross validation]] techniques. The SCV estimator differs from the plug-in estimator in the second term

These references also contain algorithms on optimal estimation of the pilot bandwidth matrix <strong>G</strong>.

: <math>\operatorname{SCV}(\mathbf{H}) = n^{-1} |\mathbf{H}|^{-1/2} R(K) +

== Computer implementation==

n^{-2} \sum_{i=1}^n \sum_{j=1}^n (K_{2\mathbf{H} +2\mathbf{G}} - 2K_{\mathbf{H} +2\mathbf{G}}

The [http://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(7) | url=http://www.jstatsoft.org/v21/i07}}</ref> in [[R programming language|R]] implements the plug-in and smoothed cross validation selectors. This example is based on the [[Old Faithful Geyser]] in Yellowstone National Park, USA. This dataset contains

+ K_{2\mathbf{G}}) (\mathbf{X}_i - \mathbf{X}_j)</math>

272 records with two measurements each: the eruption duration time (minutes) and the

waiting time until the next eruption (minutes), and is contained in the base distribution of R. This code fragment computes the kernel density estimate with the plug-in bandwidth matrix. The coloured contours correspond to the smallest region which contains that corresponding 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 data set.

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'''_AMISE.

<pre>

=== 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==

In the optimal bandwidth selection section, we introduced the MISE. Its construction relies on the [[expected value]] and the [[variance]] of the density estimator<ref name="WJ1995" />{{rp|97}}

:<math>\operatorname{E} \hat{f}(\mathbf{x};\mathbf{H}) = K_\mathbf{H} * f (\mathbf{x}) = f(\mathbf{x}) + \frac{1}{2} m_2(K)\operatorname{tr} (\mathbf{H} \operatorname{D}^2 f(\mathbf{x})) + o(\operatorname{tr} \, \mathbf{H})</math>

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''⁻¹ |'''H'''|^−1/2 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

:<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 convergence 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 convergence 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

:<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_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.

==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.

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.

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

library(ks)

data(faithful)

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)

</syntaxhighlight>

points(faithful)

</pre>

[[Image:Old==Density Faithful Geyser KDEestimation with plugina diagonal bandwidth.png|center]] matrix==

[[File:Bivariate example.png|thumb|250px|alt=Kernel density estimate with diagonal bandwidth for synthetic normal mixture data. |Kernel density estimate with diagonal bandwidth for synthetic normal mixture data.]]

We consider estimating the density of the Gaussian mixture

== References ==

{{math|(4''π'')⁻¹&thinsp;exp(−{{frac|2}} (''x''₁² + ''x''₂²))

{{Reflist}}

+ (4''π'')⁻¹&thinsp;exp(−{{frac|2}} ((''x''₁ - 3.5)² + ''x''₂²))}},

from 500 randomly generated points. We employ the Matlab routine for

[http://www.mathworks.com/matlabcentral/fileexchange/17204 2-dimensional data].

The routine is an automatic bandwidth selection method specifically designed

for a second order Gaussian kernel.<ref>{{Cite journal

| author1 = Botev, Z.I.

| author2 = Grotowski, J.F.

| author3 = Kroese, D.P.

| title = Kernel density estimation via diffusion

| journal = [[Annals of Statistics]]

| volume = 38

| issue = 5

| pages = 2916–2957

| year = 2010

| doi = 10.1214/10-AOS799

| arxiv = 1011.2602}}

</ref>

The figure shows the joint density estimate that results from using the automatically selected bandwidth.

'''Matlab script for the example'''

== External links ==

* [http://www.mvstat.net/tduong/research www.mvstat.net/tduong/research] A collection of peer-reviewed articles of the mathematical details of multivariate kernel density estimation and their bandwidth selectors.

Type the following commands in Matlab after

[http://www.mathworks.com/matlabcentral/fileexchange/17204 downloading]

and saving the function kde2d.m

in the current directory.

<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>

==Alternative optimality criteria==

The MISE is the expected integrated ''L₂'' distance between the density estimate and the true density function ''f''. It is the most widely used, mostly due to its tractability and most software implement MISE-based bandwidth selectors.

There are alternative optimality criteria, which attempt to cover cases where MISE is not an appropriate measure.<ref name="simonoff1996" />{{rp|34–37,78}} The equivalent ''L₁'' measure, Mean Integrated Absolute Error, is

: <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₁ 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_∞'' 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>

which has been investigated only briefly.<ref>{{cite journal | author1=Cao, R. | author2=Cuevas, A. | author3=Manteiga, W.G.| title=A comparative study of several smoothing methods in density estimation | journal = Computational Statistics and Data Analysis | year=1994 | volume=17 | issue=2 | pages=153–176 | doi=10.1016/0167-9473(92)00066-Z}}</ref> Likelihood error criteria include those based on the Mean [[Kullback–Leibler divergence]]

: <math>\operatorname{MKL} (\mathbf{H}) = \int f(\mathbf{x}) \, \operatorname{log} [f(\mathbf{x})] \, d\mathbf{x} - \operatorname{E} \int f(\mathbf{x}) \, \operatorname{log} [\hat{f}(\mathbf{x};\mathbf{H})] \, d\mathbf{x}</math>

and the Mean [[Hellinger distance]]

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

<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"/>

<math>\hat{f}(x) = \frac{1}{(2\pi)^d} \int \hat\varphi(\vec{t})\psi_h(\vec{t}) e^{-i\vec{t} \cdot \vec{x}}d\vec{t}</math>

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://github.com/LBL-EESA/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.]]

==See also==

* [[Kernel density estimation]]&nbsp;&ndash; Univariateunivariate kernel density estimation.

* [[Variable kernel density estimation]]&nbsp;&ndash; estimation of multivariate densities using the kernel with variable bandwidth

==References==

{{Reflist}}

==External links==

* [http://www.mvstat.net/mvksa <em>Multivariate Kernel Smoothing and Its Applications</em>] is a comprehensive book on many topics in kernel smoothing, including density estimation. Includes [https://cran.r-project.org/web/packages/ks/index.html ks package] code snippets in [[R_programming language|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]].

* [http://www.mathworks.com/matlabcentral/fileexchange/58312-kernel-density-estimator-for-high-dimensions akde.m] A [[Matlab]] m-file for multivariate, [[variable bandwidth kernel density estimation]].

* [https://github.com/thaines/helit/tree/master/ms helit] and [http://pythonhosted.org/PyQt-Fit/mod_kde.html pyqt_fit.kde Module] in the [https://pypi.python.org/packages/source/P/PyQt-Fit/PyQt-Fit-1.3.4.tar.gz PyQt-Fit package] are [[Python (programming language)|Python]] libraries for multivariate kernel density estimation.

[[Category:Estimation of densities]]

[[Category:Nonparametric statistics]]

[[Category:Computational statistics]]

[[Category:Multivariate statistics]]

[[Category:Articles with example MATLAB/Octave code]]