Multivariate locally adaptive density estimation.
Multivariate versions of variable bandwidth kernel density estimators can lead to improvement over kernel density estimators using global bandwidth choices. These estimators are more flexible and better able to model complex (multimodal) densities. In this work, two variable bandwidth estimators are discussed: the balloon estimator which varies the smoothing matrix with each estimation point and the sample point estimator which uses a different smoothing matrix for each data point. A binned version of the sample point estimator is developed that, for various situations in low to moderate dimensions, exhibits less error (MISE) than the fixed bandwidth estimator and the balloon estimator. A practical implementation of the sample point estimator is shown through simulation and example to do a better job at reconstructing features of the underlying density than fixed bandwidth estimators. Computational details, including parameterization of the smoothing matrix, are discussed throughout.
Introduction
The kernel density estimator has become a staple in the data analysis tool box largely
because of the flexibility and generality of the method. Much research has been done on the
theoretical properties of the kernel estimator and the superiority over such naive estimators
as histograms is well-established. With cross-validation (Rudemo, 1982; Bowman, 1984)
and plug-in (Sheather and Jones, 1991) rules for choosing smoothing parameters, data-
driven, automatic bandwidth selection has achieved a certain level of maturity.
Rules for choosing bandwidths are generally based on the relatively simple but important
idea of balancing bias and variance globally. This idea works well for most regular
densities, i.e. those densities that are unimodal and are not extremely skewed. Consider
the multivariate kernel density estimator given by
=n
is a d-variate random sample with density f . The kernel, K, is taken to be
a d-variate density function with
I d .
The contours of the kernel are restricted to be spherically symmetric and the smoothing
parameter, h, controls the size of the kernels.
Straightforward asymptotic approximations yield integrated squared bias (ISB) and
integrated variance (IV) equal to
nh d
and
Z
R K 2 (w)dw, tr indicates the trace of a matrix, and 5 2 f(x) is the Hessian
(matrix of second partial derivatives) of f . Combining these terms yields an estimate of
the mean integrated squared error,
dx. The optimal bandwidth
can then be easily derived and is equal to
oe 4
This choice of h yields an
For a univariate standard normal distribution, the globally optimal smoothing parameter
for a sample of size 0:421. This value reflects the best trade-off between
variance and bias for this density and sample size. Interestingly, the optimal smoothing
parameter for the region near the mode, (\Gamma1; 1), is slightly smaller than
the global bandwidth. For either of the tails (-1,-1) and (1,1), the optimal bandwidth
is larger, reflecting the sparseness of data and an attempt control the higher
variability in the tails. Ultimately, the global bandwidth is closer to the bandwidth optimal
for the mode as the largest contribution to the error comes from the area near the mode.
The effect of the tails is not substantial enough to greatly influence the global bandwidth.
As the dimensionality increases however, the so-called curse of dimensionality becomes
influential. The effect on multivariate density estimates is well noted (Friedman and Stue-
zle, 1981, and Scott and Wand, 1991). Due to the sparseness of data in higher dimensions,
multivariate neighborhoods are generally empty, particularly in the "tails" of the density
(Exercise 4.1 in Wand and Jones, 1995). Furthermore, there is less probability mass near
the mode. Not only does this have the effect of slowing the convergence rate of the MISE
as dimensionality increases, but the relative contributions of variance and bias change. In
the univariate setting, the ratio of bias to variance is 4:1. For general dimensions, this
ratio is 4:d. Hence, as the dimension increases, larger and larger bandwidths are necessary
to control the increased variability, particularly the contributions from the tails. However,
this also has adverse effect of averaging away features near the modes.
Sain and Scott (1996) give another example of where this trade-off between variance
and bias fails. Consider a bimodal normal mixture of the form
1=4OE 1=3 normal density with mean  and variance oe 2 . This
mixture is characteristic of a density that is difficult for the kernel estimator in that the
modes are of equal height but have differing scales. The globally optimal smoothing parameter
of a sample of size 0:248. Considering the sample sizes associated with
each mode, the optimal smoothing parameters are
Again, the global value reflects an attempt to find some sort of middle ground between
what is optimal for each mode. However, this value of h will undersmooth or oversmooth,
depending on the mode, and the analyst will be faced with the difficult decision of choosing
which features are real and which are noise.
Graphical ideas such those proposed by Minnotte and Scott (1993) and Chaudhuri and
Marron (1997) have proven to be effective aids for exploring structure in univariate multi-modal
densities. Both approaches use a family of kernel estimates, where the bandwidth is
allowed to range from small to large. These methods are attractive in that they allow the
user to change "resolution" by focusing on large bandwidths to gain insight into general
structure and then smaller bandwidths for finer detail. However, they serve to emphasize
the lesson that a single smoothing parameter can be ineffective for more complex densities.
There are proposals to vary the amount of smoothing in some fashion. Intuitively, the
task seems straightforward. More smoothing is necessary where data are scarce (i.e. tails
and valleys) and less smoothing is necessary near modes. However, it is not so easy as to
model the bandwidth function on the level of the underlying density alone. Some attempt
to account for the curvature must be made as well (Sain and Scott, 1996). Actually, it has
proven difficult to gain understanding of exactly how to vary the bandwidth and, further-
more, how to apply this understanding to varying the bandwidth in practice, particularly
in the multivariate case.
In this paper, variable bandwidth estimators will be studied in the multivariate setting,
building on the work of Terrell and Scott (1992) and Sain and Scott (1996). Following a brief
introduction to the issues involved with a more general fixed bandwidth multivariate kernel
estimator (Section 2), a comparison of two common approaches to variable bandwidth
estimators is presented (Sections 3 and 4). Finally, a brief discussion of some of the issues
and problems of variable bandwidth estimation in practice is presented in Section 5.
Fixed Bandwidth Estimators
The general multivariate kernel density estimator is given by
K(H \Gamma1=2
Unlike (1) the smoothing parameter, H, is actually a symmetric positive definite matrix
that is analogous to the covariance matrix of K. Hence, KH
Straightforward asymptotic analysis involving multivariate Taylor's series expansions
can yield an estimate of the MISE. Following Scott (1992) or Wand and Jones (1995), an
asymptotic approximation of the MISE is given by
Z
indicates the determinant of a matrix.
Choosing the form of H depends on the complexity of the underlying density and the
number of parameters that must be estimated. Wand and Jones (1993) give an excellent
discussion of the issues in the bivariate case. There are three primary classes for
parameterizing the smoothing matrix. Assuming that K is a multivariate normal kernel,
i.e. x=2), the first class was introduced in (1). This kernel has
contours that are restricted to be spherically symmetric and it has only one independent
smoothing parameter. In keeping with the notation of Wand and Jones (1993), this class
will be denoted as H 1 where H
The second class, H 2 , allows kernels that are ellipsoidal. However, the axes are restricted
to be parallel to the coordinate axis. Here, d independent smoothing parameters are allowed
diag indicates a diagonal matrix. This
class is commonly referred to as the product kernel estimator and allows different amounts
of smoothing in each dimension. The optimal bandwidths, h are still proportional
to n \Gamma1=(d+4) and the AMISE is of the same order as H 1 . However, some gain can be obtained
by using different smoothing parameters for each dimension, especially when the scales of
the variables differ.
The final class, in which a full smoothing matrix is employed, involves d(d
independent smoothing parameters and is denoted as H 3 . This class allows ellipsoidal
kernels of arbitrary orientation and is given by
Unfortunately, there is no closed form expression for the optimal smoothing matrix in this
case. Numerical methods are required to find H for a candidate density and a particular
sample size.
If the underlying density has spherically symmetric contours (such as the standard nor-
mal), then the optimal smoothing matrix for each class will be the same, i.e. proportional
to the identity matrix. However, it is an interesting exercise to compare the smoothing
parameterizations in other cases in order to evaluate the gain in terms of the AMISE, as
well as developing some understanding how the smoothing matrix adapts to the shape of
the density for each class. Consider the case of a bivariate normal density centered at the
origin with oe 0:6. The sample size is In this case, the
optimal bandwidth matrices for each class are given by
1:0 0:0
0:0 1:0
0:5 0:0
0:0 2:0
0:625 0:75
0:75 2:5
for smoothing matrices in H 1 , H 2 , and H 3 , respectively and a standard normal kernel. Note
that these smoothing matrices are written in the form
h controls the size of the kernel and A controls the orientation (see Scott, 1992). The most
restrictive case of a single smoothing parameter yields a smaller kernel as it tries to adjust
to the structure of the density, e.g. different scales and the correlation. Both H 2 and H 3
correctly adjust for scale as the squared smoothing parameters are four times larger in the
second dimension. The unrestricted kernel, H 3 , has the largest area suggesting a trade-off
between the flexibility to model the orientation and scale of the density and the size of the
kernel. The unrestricted and the product kernel estimators yield AMISE values that are
30% and 16% smaller than using just one smoothing parameter.
In practice, a common approach is to scale the data so that the sample variances are
the same in each dimension or to sphere the data in which a linear transformation is
applied that yields an identity sample covariance matrix. These approaches are essentially
dimension reduction ideas as they allow for the use of a single smoothing parameter
with the transformed data. However, Wand and Jones (1993) urge caution when using
these approaches as they are not guaranteed to give the correct transformation or rotation
to achieve the gains possible by using the full smoothing matrix. Many authors have
noted that for most densities, in particular unimodal ones, allowing different amounts
of smoothing for each dimension (the product kernel estimator) is adequate. With more
complex densities, especially multimodal ones, the situation is less clear, although rotations
can help if the structure of the distribution can be aligned with the coordinate axis (Wand
and Jones, 1993).
3 Locally Adaptive Density Estimation
It is conceivable to consider a bandwidth function that adapts to not only the point of
estimation, but also the observed data points and the shape of the underlying density. As
a matter of practice, however, two simplified versions have been studied. The first varies the
bandwidth at each estimation point and is referred to as the balloon estimator. The second
varies the bandwidth for each data point and is referred as the sample-point estimator. See
Jones (1990) for a detailed comparison of the two estimators in the univariate case.
3.1 Balloon Estimators
The general form of the balloon estimator is given by
where H(x) is the smoothing matrix for the estimation point x. Considered pointwise,
the construction of this estimator is the same as the fixed bandwidth estimator given in
(2). For each point at which the density is to be estimated, kernels of the same size and
orientation are centered at each data point and the density estimate is computed by taking
the average of the heights of the kernels at the estimation point.
This type estimator was introduced by Loftsgaarden and Quesenberry (1965) as the kth
nearest neighbor estimator. While not generally considered a kernel estimator directly, it
can be be written as in (4) by taking K to be a uniform density on the unit sphere and by
restricting the smoothing matrix to H 1 . Thus, the bandwidth function can be written as
which measures the distance from x to the kth nearest data point.
Much has been written about the kth nearest neighbor estimator and it seems clear
that it is not an effective density estimator in the univariate case. The bandwidth function
is discontinuous and these discontinuities manifest themselves directly in the density
estimate (Silverman, 1986). Furthermore, the estimator suffers from severe bias problems,
particularly in the tails (Mack and Rosenblatt, 1979; Hall, 1983; Terrell and Scott, 1992).
However, Terrell and Scott (1992) show that the kth nearest neighbor estimator improves
as dimensionality increases and will perform well in dimensions greater than 4.
Terrell and Scott (1992) also study error properties of the general balloon estimator and
found some remarkable results. Applied pointwise, the balloon estimator behaves just like
the fixed bandwidth estimator. The authors show that by choosing the orientation of the
smoothing matrix appropriately, the bias can be dramatically reduced. In fact, the bias
will be of higher-order. In the following, this important result of Terrell and Scott (1992)
will be outlined with some extensions not shown in that paper that will allow comparison
with the sample point estimator.
Consider using (4) pointwise with 1. For notational conve-
nience, the argument of H will be dropped in the following. Here, h controls the size of
the kernel and A the orientation. Then, an asymptotic approximation for the bias can be
written as
where S x is the matrix of second partial derivatives of the underlying density. Terrell and
Scott (1992) look at three cases based on the form of S x . To demonstrate, consider f to
be a bivariate standard normal density and K a bivariate standard normal kernel.
In the first, S x is assumed to be positive or negative definite. No choice of A will
make the asymptotic bias equal to zero and the best possible asymptotic mean square
error is O(n \Gamma4=(d+4) ), the same order as the fixed bandwidth estimator. In the example,
this case corresponds to points inside the unit ball,
A is chosen to minimize the bias and is a solution to AA
x . Given the
form of A, h is computing by minimizing the resulting pointwise AMSE, resulting in
The second case corresponds to the density being saddle shaped with some eigenvalues
of S x being positive and some negative. In this case, it can be shown that by properly
choosing A, the leading of the bias can be made to be zero. In the example, this
corresponds to points outside the unit ball, jjxjj ? 1. Choose A to be of the form
is the absolute value of the eigenvalues of S x , O is the matrix of
eigenvectors of S x , and a is chosen such that finding A in this fashion, h
can be chosen to minimize the AMSE. To determine the AMSE and the resulting value of
h, the fourth order terms in the Taylor's series expansion of the expectation in the bias
must be considered. The expectation is given by
K(h
Z
Z Z
is the ijth element of A, and f i represents the ith
partial derivative of f with respect to x i . By the assumptions on the kernel and the proper
choice of A, the expectation reduces to
f ijkl
Z Z
yielding a bias of O(h 4 ). Writing the bias as
where C(f; K) is a constant dependent on the mixed fourth-order partial derivatives of f
evaluated at x and the mixed fourth-order moments of K and combining the bias with the
asymptotic expression for the variance,
the pointwise MSE for points outside the unit ball is given by
Taking the partial derivative with respect to h, setting to zero, and solving for h gives
The third case corresponds to S x being positive or negative semidefinite, with some
eigenvalues equal to zero. Hence, the Hessian is of lower rank and the bias is of higher
order (with proper choice of A). This corresponds to points in the example on the unit
circle,
Unfortunately, there is no closed form expression for the MISE that can be computed
from these expressions for the optimal pointwise smoothing matrices. However, as Terrell
and Scott (1992) point out, the asymptotic MSE can be numerically integrated to study
the behavior of the balloon estimator.
3.2 Sample-Point Estimators
The multivariate sample-point estimator is given by
where H(x i ) is the smoothing matrix associate with x i . The sample-point estimator still
places a kernel at each data point, but these kernels each have their own size and orientation
regardless of where the density is to be estimated.
This type of estimator was introduced by Breiman, Meisel, and Purcell (1977) who
suggested using H(x i is the distance from x i to the kth nearest
data point. Asymptotically, this is equivalent to choosing h(x i d is the
dimension of the data.
Abramson (1982a,b) suggested the square-root law, i.e. setting
in practice, using a pilot estimate of the density to calibrate the bandwidth function. This
formulation of the bandwidth function has been popular in no small part due to the early
results that show that the bias associated with the square-root law was of higher order
(Silverman, 1986; Hall and Marron, 1988; Jones, 1990). However, recent work has shown
that this early result may not always hold due to bias contributions from the tails (Hall,
1992; McKay, 1993; Terrell and Scott, 1992; Hall, Hu, and Marron, 1994; Sain and Scott,
1996).
The square-root law suffers from bias properties that can diminish the gains of varying
the bandwidth asymptotically, but it also suffers from a certain inflexibility by restricting
the bandwidth function to be only a function of the height of the density. Sain and Scott
introduce a binned version of the sample-point estimator that uses a piecewise
constant bandwidth function in an attempt to provide a more general study of the properties
of the estimator. In that work, the authors showed that the estimator did not exhibit a
higher-order MISE but it did lead to substantial improvement over the fixed bandwidth
estimator in the univariate case.
A multivariate version of the binned sample-point estimator is given by
is the number of data points in the jth bin, t j is the center of the jth bin, and
is the smoothing matrix associated with the jth bin. In general, an equally spaced
mesh of points is laid down over the support of the density to define the bins, although
other binning rules such as the linear binning defined in Hall and Wand (1996) could be
considered.
Binning has been used in density estimation for a variety of reasons. Hall (1982) studied
rounded and truncated data. Silverman (1982), Hardle and Scott (1992), Wand (1994), and
Hall and Wand (1996) use binning as a device that can radically reduce computing time.
Scott and Sheather (1985) build on the result of Hall (1982) and show that binning results
in an inflated bias but that the MISE is insensitive to reasonable amounts of binning.
For adaptive estimation, binning becomes much more than a computational tool. Computing
the bias of an adaptive estimator is difficult, if not impossible, without specifying
the form of the bandwidth function as was done with the square root law. Through binning,
the expectation of the estimator in (6) is easy to compute by noting that the only random
quantities are the n These counts can be thought of as a realization of a
multinomial distribution with parameters
R
denotes that jth bin.
Assuming that K is a multivariate normal kernel and following Sain and Scott (1996),
the MISE of the multivariate binned-sample point estimator is given by
Z
Z
Z
\Gamman
Z
Z
Note that the normal integrals follow directly from formulae such as
those presented in the appendices to Wand and Jones (1995). By specifying f , usually as
some sort of normal mixture, the probabilities can be computed and then
the MISE function optimized over the collection of smoothing matrices H
4 A Comparison of Variable Bandwidth Methods
Comparing optimal bandwidth functions for the balloon and the sample-point estimator is a
difficult task, especially in the multivariate setting. However, developing an understanding
of how optimal smoothing parameters behave will yield insight, particularly when designing
practical algorithms. It would seem reasonable that there is some sort of fundamental
relationship between the two methods, at least asymptotically. Sain and Scott (1996)
noticed some similarities between the sample-point estimator and the so-called zero-bias
bandwidths discussed by Hazelton (1996), Sain and Scott (1997), and Devroye (1998).
However, for finite multivariate samples, some differences appear.
To illustrate, let f to be a bivariate standard normal and let K also be standard
normal. Figure 1 display optimal kernels for the fixed bandwidth estimator, the balloon
estimator and the binned sample-point estimator using an equally spaced mesh with 9 bins
per dimension. Two versions of the sample-point estimator are considered. One restricts
the contours of the kernels to be circular (H 1 ) while the second allows unrestricted size
and orientation (H 3 ). Optimal kernels were computed using the formulae and techniques
discussed in the previous sections.
Four cases are considered, corresponding to bins centered at the origin and along the
x-axis at 0.9, 1.8, and 2.7. These correspond to the mesh for the sample-point estimator
being laid down on the square defined by (\Gamma3; 3). Since the contours of the
density are radially symmetric, it follows that the contours of the optimal kernels centered
at different bins, but at the same distance from the origin, will behave similarly.

Figure

1: Optimal kernels for a fixed bandwidth estimator (solid line), balloon estimator
(dotted lines), sample-point estimator with the smoothing matrix restricted to H 1 (short
dashed lines), and sample-point estimator with an unrestricted smoothing matrix (long
dashed lines). Ellipses are 95% contours of the kernels.
The optimal kernel for the fixed bandwidth estimator has contours that are circular
regardless of where the density is being estimated, which data point, or which parameterization
of the smoothing matrix is used. The 95% contour of this kernel is plotted in Figure
1 as a solid line. It is centered at the origin as a reference and an indication of scale for
the variable bandwidth kernels.
When examining the optimal kernels for the balloon estimator and the sample-point
estimator, it is important to keep in mind the construction of the estimator. For the
balloon estimators, the plotted kernels are centered at the estimation point, but the density
is constructed by centering the kernels at each data point. For the binned-sample point
estimator, the density is estimated by placing a weighted kernel at the center of the bin
with each bin having a different size and orientation. In practice, the kernels are assumed
to be the same for each data point in the bin.
First, consider the kernels at the origin (jjxjj = 0:0) in Figure 1. All of the kernels are
circular which is no surprise considering the structure of this density at the origin. However,
all of the variable bandwidth kernels are smaller than the fixed bandwidth kernel. In the
next plot, 0:9, the center of the bin is still within the unit circle. The fixed bandwidth
kernel is still much larger than the sample-point kernels which are still roughly circular.
The balloon estimator kernel is already adjusting to the local curvature and the major axis
of the ellipse is parallel to the line between the origin to the center of the bin.
The third plot, shows the kernels for the bin just outside the unit circle. Now
the variable bandwidth kernels are much larger than the fixed bandwidth kernel. However,
a distinct change in orientation is noted between the balloon kernel and the unrestricted
sample-point kernel. Recall the balloon kernel has higher-order bias in this region and is
based on the fourth-order derivatives of the density. Hence, the orientation changes, with
the major axis of the ellipse being perpendicular to the line between the origin and the
center of the bin. The major axis of the unrestricted kernel retains the lower-order bias
behavior with the major axis of the kernel parallel to the origin.
Moving even further out into the tails of the density, 2:7, the variable bandwidth
kernels dwarf the fixed bandwidth kernel, with the restricted sample-point kernel being
exceptionally large. This phenomenon is explained by the sample-point estimator with the
restricted kernels attempt to minimize variability in the tails by trading a much larger size
for the correct orientation. Note that the kernel for the sample-point estimator is much
larger than the balloon estimator kernel.
The binned sample-point estimator allows head-to-head comparisons between a general
sample-point estimator and the fixed bandwidth estimator or the balloon estimator. In
actuality, the binned sample-point estimator is not entirely general in that the smoothing
matrix is assumed to be constant for all data points in a particular bin. However, through
the optimization, the smoothing matrices are allowed to adapt to level and curvature,
something that the square-root law cannot accomplish.
Estimator Bins per Dimension
Fixed

Table

1: MISE values for fixed and sample-point adaptive methods using a normal kernel
with H 1 and a N(0; I d ) density. Sample size is
Scott and Sheather (1985) and Hall and Wand (1996) showed that binning can inflate
the bias. Sain and Scott (1996) showed that the binned sample-point estimator needed
enough bins to counteract this bias and give sufficient flexibility to improve on the fixed
bandwidth estimator. The same is true in the multivariate case. In fact, Table 1 shows that
more bins per dimension may be needed to achieve that same gain. Table 1 is based on a
multivariate standard normal density and 100. For the fixed bandwidth method, the
MISE was computed using the expressions derived by Worton (1989) and Marron and Wand
(1992). The kernels for the binned sample-point estimator were restricted to be spherically
simplicity. An equally spaced mesh was laid down over the approximate
support of the density (assumed to be (-3,3) in each dimension for consistency).
What is clear from Table 1 is that considerable gain in terms of the MISE can be
achieved with a relatively small number of bins per dimension and restricting the shape
of the kernels. For using only 5 bins was actually worse than the fixed bandwidth
estimator. However, for using 5 bins per dimension (125 total bins) was able to
achieve some gain over the fixed bandwidth estimator.
Using bins achieved a 44% gain in MISE for
This suggests that it may be necessary to use more bins per dimension or use the more
general smoothing matrices for each bin. In Table 2 MISE values are also computed for
the bivariate standard normal. However, a sample-point estimator included
that allows unrestricted smoothing matrices. A considerable gain in terms of the MISE is
realized as the restriction on the kernels is removed. The restricted kernel sample-point
estimator with bins leads to a 29% improvement in the MISE while
the unrestricted smoothing matrices leads to a 51% reduction in the MISE.
Fixed
Balloon

Table

2: MISE values for fixed, sample-point, and balloon estimators using a normal kernel
and a N(0; I 2 ). For the sample-point estimator,
The binned sample-point estimator also allows a direct comparison with the balloon
estimator. In this case, the unrestricted sample-point estimator actually gives a slight
improvement over the balloon estimator. This implies that the sample-point estimator is
doing a better job even near the mode which is the only place the estimators actually
compete as the balloon estimator is of higher order in the tails. The results also suggest
that it is not a small sample curiosity as the improvement is also achieved for the more
moderate sample size of n = 1000. However, for larger d, the results will likely change
as the impact of the tails (where the balloon estimator is of higher order) becomes more
pronounced.
While varying the bandwidth can aid in combating the curse of dimensionality, the
real power of the sample-point estimator may lie in the ability to model more complex,
multimodal distributions. Consider the normal mixture density,
labeled "K" in Wand and Jones (1993). This is a bivariate, trimodal density comprised of
three normal components with parameters
p3
p3
Contours of this mixture are shown in the left frame of Figure 2. This density poses
problems for the fixed bandwidth kernel, even with unrestricted orientation because the
covariance structure of the entire density is not far from the identity. The modes, however,
have varying heights and different covariance structures.

Table

3 shows a comparison of the MISE for a fixed kernel estimator with a smoothing
Estimator MISE
Fixed

Table

3: MISE values for fixed and sample-point estimators using a normal kernel
and a trimodal normal mixture. Sample size is
matrix in H 1 and a unrestrictive smoothing matrix in H 3 . The sample point estimator
with bins per dimension is also considered for the two smoothing matrix struc-
tures. The results for the two parameterizations of the fixed bandwidth estimator are not
terribly different with a slight improvement going to the unrestricted smoothing matrix.
However, allowing the kernels to vary can lead to significant reduction in the MISE. Even
just allowing the size of the kernels to vary can give significant gains with the restricted
sample-point estimator giving a 35% reduction in the MISE. Contours of the optimal band-width
function are shown in the right frame of Figure 2. The roughness of the contours
is due to interpolating the piecewise constant bandwidth function. However, it is clear
from the contours that the size of the kernels mirror the shape of the density. What is
difficult to discern from the plot, however, is that the kernels are not adapting simply based
on heights alone. Curvature is also taken into account as the kernels near the left mode
(first component of the mixture) are much smaller than the those near the top right mode
(second component of the mixture) despite the similar heights of these modes. Finally, the
unrestricted sample-point estimator yields even more promising results and has a MISE
that is 62% smaller than the unrestricted fixed estimator.
5 Practical Algorithms
The previous section shows that varying the bandwidth in some fashion can lead to substantial
theoretical gains. Designing a practical, data-based algorithm to realize these gains
in practice is a difficult task. There are many more things to consider than in the univariate
case, not least of these is estimating the bandwidth function and parameterization of the
Figure

2: Contours of normal mixture distribution (left frame) and contours of the optimal
bandwidth function (right frame) for the binned sample-point estimator
kernels.
For example, the balloon estimator is attractive for a variety of reasons. It has tremendous
potential when good pointwise estimates of the density in the tails are required. How-
ever, the orientation of the balloon kernel leading to a higher-order bias requires knowledge
of the fourth order derivatives. A pilot estimate could be used to calibrate the kernels, as
suggested by Terrell and Scott (1992), but estimating fourth-order derivatives is generally
perceived as harder than estimating the density itself and the gains from adaptivity would
be overwhelmed from the error in estimating the derivatives.
For the sample-point estimator, there are currently no asymptotic approximations to
the MISE that would allow a plug-in style bandwidth selector. This is due in part to
the traditional asymptotic approximations and the inability to model the behavior of the
optimal bandwidth function (Sain and Scott, 1996). However, least-squares cross-validation
offers a method that is an unbiased estimate of the MISE and does not rely on an asymptotic
approximation. It was also shown in Sain and Scott (1996) to perform well in the univariate
setting at estimating the optimal bandwidth function. Hence, it will be used as a test case
to illustrate some of the difficulties that may be encountered.
Least-squares cross-validation, also referred to as unbiased cross-validation (UCV), was
introduced by Rudemo (1982) and Bowman (1984) and is an approximation of the integrated
squared error. On average, this criterion is equal to the MISE; hence the term
unbiased cross-validation. In this setting, the UCV function is given as
where
and
Note that this is not a fully binned version of UCV because the actual data are used in
the estimation of the cross-product. As with the expression for the MISE, the bin counts
are computed and then the criterion minimized over the space of the smoothing matrices,
m.
As an illustration, consider the lipid data of Scott, et al. (1978). These data consist of
measurements of cholesterol and triglycerides for 320 men diagnosed with coronary artery
disease. The original paper showed that the data were bimodal, indicating an increased
risk for heart disease associated with increased cholesterol level. Further study of the data
reveals the potential of a third mode. Figure 3 shows a scatterplot of the data with the
mesh used to bin the data for the sample point estimator shown by the dotted lines. An
11\Theta11 mesh was used; only the bins with data are indicated on the plot. Note that the
data were scaled to have zero mean and unit variance in each dimension.

Figure

4 shows perspective and contour plots for a fixed bandwidth estimator using a
smoothing matrix in H 1 . The smoothing parameter was selected by cross-validation. The
estimate shows evidence of two distinct modes, while the third mode is marginal. There is
also a considerable amount of noise exhibited by spurious modes in the tails.
Figure

3: Scatterplot of cholesterol and triglyceride levels for males with heart
disease. The grid indicates the mesh used for the sample-point estimator. The data have
been scaled to have zero mean and unit variance in each dimension.
-22
-22
Figure

4: Perspective and contour plots of the fixed bandwidth (H 1 ) estimate for the lipid
data.
The bandwidths for the sample-point estimate with smoothing matrices restricted to
are shown in Figure 5. The observed pattern is consistent with intuition, having larger
bandwidths in the tails and smaller near the modes. The estimate is shown in Figure 6 and
initially appears oversmoothed. However, the modes are actually more pronounced (higher)
and there are three clear modes in the adaptive estimate. There is also little evidence of
the spurious modes in the tails of the distribution. The resulting estimator achieves the
goals of the adaptive procedure, namely to correctly emphasize features (both heights and
positions of modes) while minimizing noise and spurious modes.

Figure

7 shows the sample-point estimate using unrestricted kernels (H 3 ) chosen by
UCV. Unfortunately, this estimate does not reflect the promise of the MISE results in the
previous section. The estimate is quite rough with many spurious modes. To determine
the cause, consider an examination of the estimated kernels. Those kernels near the modes
are shown in Figure 8. Here, the source of the difficulties can be found. Many of the
appear close to degenerate and most seem to have wildly varying orientations.
This may be an artifact of the well-established aggressive behavior of cross-validation to
yield highly variable bandwidths that are often much smaller than optimal. Unfortunately,
there are just too many parameters to estimate, in this case 48\Theta3=144, for such a moderate
sample size. Clearly, attempting to estimate the fully parameterized smoothing matrices
will require much more data or a more stable procedure.
One way to reduce the dimensionality of the problem is based on the notion that the
fixed bandwidth procedure is adequate for unimodal distributions. Hence a pilot estimate
could be used to calibrate the bandwidth function, using the same bandwidth for each data
point near a sample mode. An example is shown in Figure 9. Here the fixed bandwidth
estimate from Figure 4 is used to partition the data. There are 16 sample modes in the
estimate. Each data point is assigned to the nearest mode and each mode is assigned
a separate smoothing parameter. This has the effect of defining a piecewise constant
bandwidth function; only now the mesh is not equally spaced.

Figure

9 shows the estimated bandwidths using the restricted smoothing matrices (H 1 ).
The behavior is as expected with smaller bandwidths near dominant modes and larger
bandwidths out in the tails. The estimate using these bandwidths is plotted in Figure
2.59 1.86
2.18 1.73 0.7 0.8 1.55 2.09
2.09 2.02
Figure

5: Estimated bandwidths for the restricted sample-point estimate.
-22
Figure

Perspective and contour plots of the sample-point estimate (H 1 ) for the lipid
data.
Figure

7: Perspective and contour plots of the estimated density using a fully parameterized
bandwidth matrix (H 3 ).
-22
Figure

8: Estimated kernels for fully parameterized bandwidth matrix.
1.378 2.2210.3480.4332.9562.3431.901
Figure

9: Left plot: modes and partitions for calibrating the smoothing matrix. Right plot:
estimated bandwidths based on the partitions.
-22
Figure

10: Perspective and contour plots of the estimated density using the bandwidth
function calibrated by the pilot estimate.
and it shows the two modes clearly and the third mode is apparent. Most of the
variability in the tails is also diminished. While not as good as the sample-point estimate
shown in Figure 6, this estimate is still an improvement over the fixed bandwidth estimate,
particularly in the tails where the number of spurious modes is reduced. One interpretation
of the result is that some improvement can be gained by using a more restrictive bandwidth
function that puts the adaptivity where it may be needed the most, that is out in the tails.
6 Conclusions
There are two primary points to be made from this work. First, locally adaptive density
estimation that allows smoothing to vary in some fashion can lead to substantial gains over
fixed bandwidth density estimation and, to some extent, combat the curse of dimensionality.
The results are even more promising when the underlying density is complex, exhibiting
multiple modes with differences in scale and orientation. Second, it is possible to achieve
gains in practice, but it is not necessarily easy and overparameterizing the smoothing
matrices is a concern.
More work is certainly required in understanding the complex theory behind locally
adaptive methods and their connections to mixture models. More work is certainly required
in building practical algorithms. Least-squares cross-validation holds a lot of promise, but
it certainly has problems as well.
This research is blurring the lines between parametric and nonparametric density es-
timation. What has been proposed here is certainly inspired by kernel estimators and
motivated by a desire to learn more about variable bandwidth kernel estimators and applications
in higher dimensions. Kernel estimators already possess a great flexibility, and
locally adaptive kernel estimators add significantly to an already rich class. There is some
concern as the kernels that not necessarily local. However, there are not any additional
assumptions on the underlying density.
What has been proposed here also has the flavor of mixture models. Inspired by the
desire to parametrically model data that appears to be generated from several subpopula-
tions, mixture models have the surprising ability to model distributions whose components
are not necessarily in the parametric family of the components of the mixture (see Sain,
et al., 1998 for an example in practice, among others). There is other work that is also
exploring the connection between kernels and mixtures, for example Priebe (1994) and
Marchette, et al. (1996).
Finally, much of the work done here required some sort of numerical optimization. Most
of this optimization was done in S-Plus using the function nlminb (Statistical Sciences,
1995) which is based on subroutines from the PORT Mathematical Subroutine Library
(A.T.&T. Bell Laboratories, 1984). In some cases the gradient was also computed and
used in the optimization. As expected, optimizing these criteria is time and cpu intensive,
to the point of making these procedures unattractive in some cases for real-time estimation.



--R


"On Bandwidth Variation in Kernel Estimates - A Square Root Law,"
"Arbitrariness of the Pilot Estimator in Adaptive Kernel Methods,"
"An Alternative Method of Cross-Validation for the Smoothing of Density Estimates,"
"Variable Kernel Estimates of Multivariate Densities,"
"SiZer for Exploration of Structures in Curves,"
"Variable Kernel Estimates: On the Impossibility of Tuning the Parameters,"
"Projection Pursuit Regression,"
"The Influence of Rounding Errors on Some Nonparametric Estimators of a Density and its Derivatives,"
"On Near Neighbor Estimates of a Multivariate Density,"
"On Global Properties of Variable Bandwidth Density Estimators,"
"On the Accuracy of Binned Kernel Density Estimators,"
"Improved Variable Window Kernel Estimates of Probability Densities,"
"Smoothing by Weighted Averaging of Rounded Points,"
"Bandwidth Selection for Local Density Estimation,"
"Variable Kernel Density Estimates and Variable Kernel Density Estimates,"
"A Nonparametric Estimate of a Multivariate Density Function,"
"A Note on the Bias Reduction in Variable Kernel Density Estimates,"
"Multivariate K-Nearest Neighbor Density Esti- mates,"
"Filtered Kernel Density Estimation,"
"Exact Mean Integrated Squared
"The mode tree: A Tool for Visualization of Nonparametric Density Features,"
"Adaptive Mixtures,"
"Empirical Choice of Histogram and Kernel Density Estimators,"
"Outlier Detection from a Mixture Distribution when Training Data are Unlabeled,"
"On Locally Adaptive Density Estimation,"
"Zero-Bias Locally Adaptive Density Estimators,"
Multivariate Density Estimation: Theory
"Plasma Lipids as Collateral Risk Factors in Coronary Heart Disease - A Study of 371 Males with Chest Pain,"
"Feasibility of Multivariate Density Estimates,"
"A Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation,"
"Kernel Density Estimation using the Fast Fourier Transform,"
Density Estimation for Statistics and Data Analysis
"S-PLUS Guide to Statistics and Mathematical Analysis,"
"Variable Kernel Density Estimation,"
"Fast Computation of Multivariate Kernel Estimators,"
"Comparison of Smoothing Parameterizations in Bivariate Kernel Density Estimation,"
Kernel Smoothing
"Optimal Smoothing Parameters For Multivariate Fixed and Adaptive Kernel Methods,"
--TR
On the accuracy of binned kernel density estimators

--CTR
Anders Adamson , Marc Alexa, Anisotropic point set surfaces, Proceedings of the 4th international conference on Computer graphics, virtual reality, visualisation and interaction in Africa, January 25-27, 2006, Cape Town, South Africa
Yaser Sheikh , Mubarak Shah, Bayesian Modeling of Dynamic Scenes for Object Detection, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.11, p.1778-1792, November 2005
binning;smoothing matrix;mean integrated squared error;cross-validation;variable band-width kernel