Multivariate kernel density estimation

Kernel density estimation is one of the most popular techniques for density estimation i.e., estimation of probability density functions, which is one of the fundamental questions in statistics. 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^[1][2] 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.^[3]

Motivation

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.

 

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

The previous figure is a graphical representation of kernel density estimate. We now define it in an exact manner. Let ${\displaystyle {\mathbf {X}}_{1},{\mathbf {X}}_{2},\dots ,{\mathbf {X}}_{n}}$  be a d-variate random sample drawn from a common density function f. The kernel density estimate is defined to be

${\displaystyle {\hat {f}}_{\mathbf {H}}({\mathbf {x}})=n^{-1}|{\mathbf {H}}|^{-1/2}\sum _{i=1}^{n}K_{\mathbf {H}}({\mathbf {x}}-{\mathbf {X}}_{i})}$ 

where

• ${\displaystyle {\mathbf {x}}=(x_{1},x_{2},\dots ,x_{d})^{T}}$ , ${\displaystyle {\mathbf {X}}_{i}=(X_{i1},X_{i2},\dots ,X_{id})^{T},i=1,2,\dots ,n}$  are d-vectors
• K is the kernel function which is a symmetric density function, with ${\displaystyle K_{\mathbf {H}}({\mathbf {x}})=|{\mathbf {H}}|^{-1/2}K({\mathbf {H}}^{-1/2}{\mathbf {x}})}$ 
• H is the bandwidth (or smoothing) matrix which is a symmetric, positive definite d x d matrix.

The choice of the kernel function K is not crucial to the accuracy of kernel density estimators, so we use the standard multivariate normal or Gaussian density function as our kernel K throughout. Whereas the choice of the bandwidth matrix H is the single most important factor affecting its accuracy ^[4](pp. 36-39). (The bandwith for a kernel density estimate is analogous to the binwidth of histograms).

${\displaystyle K({\mathbf {x}})=(2\pi )^{-d/2}\exp(-{\tfrac {1}{2}}\,{\mathbf {x}}^{T}{\mathbf {x}}).}$ 

Optimal bandwidth matrix selection

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

${\displaystyle \operatorname {MISE} ({\mathbf {H}})=E\int [{\hat {f}}_{\mathbf {H}}({\mathbf {x}})-f({\mathbf {x}})]^{2}\,d{\mathbf {x}}.}$ 

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

${\displaystyle \operatorname {AMISE} ({\mathbf {H}})=n^{-1}|{\mathbf {H}}|^{-1/2}R(K)+{\tfrac {1}{4}}m_{2}(K)^{2}(\operatorname {vec} ^{T}{\mathbf {H}}){\mathbf {\Psi }}_{4}(\operatorname {vec} ^{T}{\mathbf {H}})}$ 

where

• ${\displaystyle R(K)=\int K({\mathbf {x} })^{2}\,d{\mathbf {x} }}$ , with ${\displaystyle R(K)=(4\pi )^{-d/2}}$  when ${\displaystyle K}$  is a normal kernel
• ${\displaystyle \int {\mathbf {x}}{\mathbf {x}}^{T}K({\mathbf {x}})^{2}\,d{\mathbf {x}}=m_{2}(K){\mathbf {I}}_{d}}$ , with I_d being the d x d identity matrix, with m₂ = 1 for the normal kernel
• ${\displaystyle \operatorname {D} ^{2}f}$  is the d x d Hessian matrix of second order partial derivatives of ${\displaystyle f}$ 
• ${\displaystyle {\mathbf {\Psi } }_{4}=\int (\operatorname {vec} \,\operatorname {D} ^{2}f({\mathbf {x} }))(\operatorname {vec} ^{T}\operatorname {D} ^{2}f({\mathbf {x} }))\,d{\mathbf {x} }}$  is a d² x d² matrix of integrated fourth order partial derivatives of f
• vec is the vector operator which stacks the columns of a matrix into a single vector e.g. ${\displaystyle \operatorname {vec} {\begin{bmatrix}a&c\\b&d\end{bmatrix}}={\begin{bmatrix}a&b&c&d\end{bmatrix}}^{T}.}$ 

This formula of the AMISE is due to ^[4](p. 97). The quality of the AMISE approximation to the MISE is given by

${\displaystyle \operatorname {MISE} ({\mathbf {H}})=\operatorname {AMISE} ({\mathbf {H}})+o(n^{-1}|{\mathbf {H}}|^{-1/2})+O(\operatorname {tr} \,{\mathbf {H}}^{2})}$ 

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 n → ∞.

The ideal optimal bandwidth selector is ${\displaystyle {\mathbf {H}}_{\operatorname {AMISE} }=\operatorname {argmin} _{{\mathbf {H}}\in F}\,\operatorname {AMISE} ({\mathbf {H}})}$  where F is the space of all symmetric, positive definite matrices. Since this ideal selector contains the unknown density function f, 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.

Plug-in

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

${\displaystyle \operatorname {PI} ({\mathbf {H} })=n^{-1}|{\mathbf {H} }|^{-1/2}R(K)+{\tfrac {1}{4}}m_{2}(K)^{2}(\operatorname {vec} ^{T}{\mathbf {H} }){\hat {\mathbf {\Psi } }}_{4}({\mathbf {G} })(\operatorname {vec} \,{\mathbf {H} })}$ 

where ${\displaystyle {\hat {\mathbf {\Psi }}}_{4}({\mathbf {G}})=n^{-1}\sum _{i=1}^{n}\sum _{j=1}^{n}[(\operatorname {vec} \,\operatorname {D} ^{2})(\operatorname {vec} ^{T}\operatorname {D} ^{2})]K_{\mathbf {G}}({\mathbf {X}}_{i}-{\mathbf {X}}_{j})}$ . Thus ${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {PI} }=\operatorname {argmin} _{{\mathbf {H} }\in F}\,\operatorname {PI} ({\mathbf {H} })}$  is the plug-in selector^[5][6]. These references also contain algorithms on optimal estimation of the pilot bandwidth matrix G.

Smoothed cross validation

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

${\displaystyle \operatorname {SCV} ({\mathbf {H} })=n^{-1}|{\mathbf {H} }|^{-1/2}R(K)+n^{-2}\sum _{i=1}^{n}\sum _{j=1}^{n}(K_{2{\mathbf {H} }+2{\mathbf {G} }}-2K_{{\mathbf {H} }+2{\mathbf {G} }}+K_{2{\mathbf {G} }})({\mathbf {X} }_{i}-{\mathbf {X} }_{j})}$ 

Thus ${\displaystyle {\hat {\mathbf {H} }}_{\operatorname {SCV} }=\operatorname {argmin} _{{\mathbf {H} }\in F}\,\operatorname {SCV} ({\mathbf {H} })}$  is the SCV selector^[7][8]. These references also contain algorithms on optimal estimation of the pilot bandwidth matrix G.

Computer implementation

The ks package^[9] in 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 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, Hpi is replaced with Hscv. This is not displayed here since it is mostly similar to the plug-in estimate for this data set.

library(ks)
data(faithful)
H <- Hpi(x=faithful)
fhat <- kde(x=faithful, H=H)
plot(fhat, display="filled.contour2")
points(faithful)

 

References

1. Rosenblatt, M. (1956). "Remarks on some nonparametric estimates of a density function". Annals of Mathematical Statistics. 27: 832–837. doi:10.1214/aoms/1177728190.
2. Parzen, E. (1962). "On estimation of a probability density function and mode". Annals of Mathematical Statistics. 33: 1065–1076. doi:10.1214/aoms/1177704472.
3. Simonoff, J.S. (1996). Smoothing Methods in Statistics. Springer. ISBN 0387947167.
4. Wand, M.P; Jones, M.C. (1995). Kernel Smoothing. London: Chapman & Hall. ISBN 0412552701. Cite error: The named reference "WJ1995" was defined multiple times with different content (see the help page).
5. Wand, M.P.; Jones, M.C. (1994). "Multivariate plug-in bandwidth selection". Computational Statistics. 9: 97–177.
6. Duong, T.; Hazelton, M.L. (2003). "Plug-in bandwidth matrices for bivariate kernel density estimation". Journal of Nonparametric Statistics. 15: 17–30. doi:10.1080/10485250306039.
7. Hall, P.; Marron, J.; Park, B. (1992). "Smoothed cross-validation". Probability Theory and Related Fields. 92: 1–20. doi:10.1007/BF01205233.
8. Duong, T.; Hazelton, M.L. (2005). "Cross validation bandwidth matrices for multivariate kernel density estimation". Scandinavian Journal of Statistics. 32: 485–506. doi:10.1111/j.1467-9469.2005.00445.x.
9. Duong, T. (2007). "ks: Kernel density estimation and kernel discriminant analysis in R". Journal of Statistical Software. 21(7).

External links

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

See also

• Kernel density estimation Univariate kernel density estimation.