J. R. Statist. Soc.

B (1991)
53,IVo. 3,pp. 683-690

A Reliable Data-based Bandwidth Selection Method for Kernel Density

Estimation
By S. J. SHEATHER and M. C. JONESt
University ofNew South Wales, Sydney, Australia IBMResearch Division, Yorktown Heights, USA

[Received August 1989. Final revision July 1990)

Downloaded from [Link] by guest on 17 May 2026

SUMMARY
We present a new method for data-based selection of the bandwidth in kernel density
estimation which has excellent properties. It improves on a recent procedure of Park and
Marron (which itself is a good method) in various ways. First, the new method has superior
theoretical performance; second, it also has a computational advantage; third, the new
method has reliably good performance for smooth densities in simulations, performance
that is second to none in the existing literature. These methods are based on choosing the
bandwidth to (approximately) minimize good quality estimates of the mean integrated
squared error. The key to the success of the current procedure is the reintroduction of a non-
stochastic term which was previously omitted together with use of the bandwidth to reduce
bias in estimation without inflating variance.
Keywords: ADAPTIVE CHOICE; BIAS REDUCTION; FUNCTIONAL ESTIMATION; SMOOTHING;
SQUARED ERROR LOSS FUNCTIONS

1. INTRODUCTION
There is currently much interest in the problem of providing good data-based
procedures for selecting the smoothing parameter-which controls the degree of
smoothing applied to the data-employed in statistical curve estimation techniques.
This is particularly so for nonparametric probability density function estimation by
the kernel method, which is described in Silverman (1986), for example. In this
context, we call the smoothing parameter the bandwidth and denote it by h.
An important recent paper in this area is Park and Marron (1990). Park and
Marron studied an estimator, the performance of which-in both theory and
simulations-proved to be clearly superior to other methods currently popular in the
literature, such as the well-known 'least squares cross-validation', for estimating
smooth densities. The work of the current paper makes further advances to the
methodology over and above those already made by Park and Marron (1990).
The improvements that we make are on several levels. First, on the theoretical side,
we improve the asymptotic rate of convergence of the estimated bandwidth to its
theoretical (but practically unavailable) optimum value. Also, and this point is
important because the rate improvement is only slight, we note that the constant
coefficient of the leading term in the expansion of the performance measure we use is
considerably reduced. Both Park and Marron's (1990) bandwidth selection procedure
and our favoured one require the numerical solution of an equation; it appears that

sAddress for correspondence: Department of Statistics, The Open University, Walton Hall, Milton Keynes, MK7
6AA, UK.

© 1991 Royal Statistical Society 0035-9246/91/53683 $2.00

684 SHEATHER AND JONES [No.3,
our procedure is easier to compute because the function that we need to zero is the
better behaved of the two. Finally, and most importantly, it must be stressed that the
bandwidth estimator that we recommend has a practical performance second to none
in the existing literature on the subject. Although the simulation study presented here
is rather small, we can be confident of our claim of excellent practical results because
of extensive recent, but as yet unpublished, simulations of J. S. Marron which
confirm the reliably good performance of our proposed procedure.

Downloaded from [Link] by guest on 17 May 2026

2. PARK AND MARRON'S h SELECTOR
The usual kernel density estimate fh of a univariate density f based on a random
sample Xl> ... , X; of size n is
n
fix) = n :' L; h- I K{h-I(x-XJ}. (1)
i=1

The bandwidth h has already been introduced in Section 1; the function K is the kernel
function which we take to be a symmetric probability density. All the data-based h
selection procedures discussed in this paper are based on choosing h to (at least
approximately) minimize a kernel-based estimate of mean integrated squared error
(MISE) via the first two terms of its usual asymptotic expansion (AMISE) valid as
n -. 00 and h=h(n) -. 0:

AMISE(h) = (nh)-I R(K) + ih 4ai R(f") (2)

(e.g. Silverman (1986), section 3.3). Here, the notation follows the convention R (g) =
Jg 2(X) dx and a; = Jx 2 g(x) dx for appropriate functions g, and the quantities
involved are assumed to exist and be finite. Each objective function is thus ofthe form
1/;(h) = (nh)-I R(K) + ih 4ai Sea) (3)
where Sea) is a kernel-based estimate of R(f"), using some appropriate bandwidth a;
Sea) is discussed below. Note that if a did not depend on h, the minimization of 1/;
could be performed analytically to give
h= [R(K)/{aiS(a)}]1/5 n-1/5, (4)
an estimate of the usual expression for the asymptotically optimal bandwith, h., say.
In its turn, h. is an approximation to h o, the exact minimizer of MISE(h).
Differing in a negligible way from Park and Marron (1990), their h estimator takes
Sea) to be

SND(a) = {n(n-I)}-l a - 5 L;L; LiV{a-I(Xi-Xj )} (5)

i*j

(this derives from the estimate n -I kl~V(Xi) of R(f") and the subscript ND, standing
for 'no diagonals', refers to the fact that the double sum does not include i = j terms).
Notice that, importantly, a is another bandwidth differing from h, and also that Lis
another symmetric density not necessarily K. Taking a to be an estimate of the
asymptotically optimal bandwidth al> say, for estimating R(f") leads to Park and
1991] DATA-BASED BANDWIDTH SELECTION METHOD 685
Marron's (1990) 'plug-in' estimator of ho, which develops ideas in Hall (1980) and
Sheather (1983, 1986). From Hall and Marron (1987),
al = CI(L) C 2(f)n -2113
= C 3(L ) C4(f)h~0I13 (6)
where

and

Downloaded from [Link] by guest on 17 May 2026

C2(f ) = {R(j)/R 2(f"')}1/13 = R- 2I13(f") Cij)·
Now, al in turn depends on an unknown functional of j. However, at this second
stage, it turns out to be sufficiently good to estimate this functional less well, in fact
using a scale modelforj, i.e. writej =g~, say, whereg>.(x) = A. -I gl(A.-Ix), gl is a fixed
density, such as the (suitably standardized) normal, and the scale parameter is
estimated robustly, giving A, say. Park and Marron's (1990) published algorithm,
with which we compare the estimators to follow, is completed by considering equation
(6) to give a general relationship between a and h, namely a = a(h) =
C3(L) C4(g»h lOl13 , i.e. a is taken to depend ~n h; then fil , say, is given by that h, found
numerically, that solves equation (4) with h = h, S = SNO and a = a(h).
Theorem 3.3 of Park and Marron (1990) describes the asymptotic performance of
iii. They show that, for sufficiently smoothj,
Ii/h o = 1 + Op(n -4/13). (7)
Moreover, the mean squared relative error (MSRE) of iiI is given by
sg(lil / h o- l )2 = 100- 1 R -2(f"){18 R(L iv ) R(f)}4/13
x {(11 R (f"')} 18/13( Q4 + 4Q - 9/9)n - 8113 (8)
where Q = C4( g>.)/C 4(f ) (essentially as in Park and Marron (1990)). Note that the
above rate of convergence is the best exhibited by Park and Marron and that, by and
large, this superior performance carries over to (small sample) simulation results as
well. Indeed, it is difficult to find another example of an automatic bandwidth
selection procedure with as consistently good a simulation performance as iiI in the
literature.

3. AN IMPROVED h SELECTOR
Given that the estimates of R(f") studied in Jones and Sheather (1991) improve,
theoretically, on those of Hall and Marron (1987) on which Park and Marron's h
selector is based, it is now natural to replace SNO by Jones and Sheather's So in the
above. The latter estimate of R (f") is given by equation (5) with i = j terms added in
(the subscript D means 'diagonals in'). The difference between So and SNO is a non-
stochastic term which contributes a positive amount to the bias in estimating R (f");
the trick is then to recognize that the bias due to the smoothing is negative and to use
the bandwidth a to (approximately) cancel the 'diagonal' term with the leading
smoothing bias term. The value of a that this leads to is a2, say, which, from equation
(4) of Jones and Sheather (1991), is
686 SHEATHER AND JONES [No.3,
cxz = DI(L) R -1/7(f"')n -1/7 (9)
where
D,(L) = {2 L iv(0)/a1J1/7.
By analogy with Park and Marron's (1990) algorithm, our first h, denoted by hzs ,
say, solves the equation
h = [R(K)/{a1So(cxz(h»}F/5n-1I5,
where, noting that CXz = CI h~/7 for appropriate CI, we have written CXz = cxz(h) = CI h 5/7,

Downloaded from [Link] by guest on 17 May 2026

where CI estimates CI. It is tempting simply to use a scale model to estimate CI but this is
not quite good enough here. The inconsistency of such an estimate (unless g>-. is fortui-
tously the true f) means that the leading bias terms do not cancel out sufficiently well
and a non-negligible bias term is reintroduced. To obtain the necessary sufficiently
small bias, it turns out that we must estimate R (ffll) by some Tsuch that T = R(f"') +
op(n -1114) (Jones and Sheather, 1991). Thus, any of the consistent estimators of
R (f"') discussed by Hall and Marron (1987) or Jones and Sheather (1991) will suffice.
There are, however, various alternative options, each of which is investigated in the
simulation study of Section 4. First, note that a computationally simpler and
essentially asymptotically equivalent alternative estimator of ho arises by leaving
CXz = czn -117, for appropriate Cz as in equation (9), estimating Cz and using equation (4)
to give a direct formula for Ii. Call this hzp, say. A third approach uses essentially the
same device as does hzs in that we write CXz in terms of h, but uses it in the minimization
of expression (3) rather than the solution of equation (4); this gives hZM • Asymptotic
analysis (not given) actually yields a slightly different optimal cx in the minimization
case which we utilize in practice. Although it is not theoretically necessary to do so, we
find it best in practice to use optimal bandwidths in diagonals-in formulae for
estimating functionals at the second stage, and estimate only the scale in these, using
the same scale model ideas as in Section 2. We give full implementation details only
for what proves to be the most successful of these options in Section 5. At times, we
refer to any of the above estimators as hz•
At the (slight) expense of introducing a third stage in the above estimation
procedure, we have obtained estimators hz of h owith better theoretical properties than
hI. This was proved by Jones and Sheather (1991). Their result 2 gives that, for flittle
smoother than that necessary for equations (7) and (8),
(10)
and that the MSRE of hz is
.5;f'(hz/ h o- l )z = 25- 1 x 2- 217 R -z(f") R(LiV)
x R (f){ a1 R (f''')/ L iv (O)} 9/7n - 5/7 (11)

for hzs and hzp, while hZM also has convergence rate (10) but proves to have an inferior
constant coefficient of MSRE to that in approximation (11) (not given). Of course
although the Op(n- 5/ 14) term appearing in equation (10) guarantees hz's superiority
over hI-which has a corresponding Op(n -4113) term in equation (7)-for sufficiently
large n, it is not so clear what the small sample repercussions are. However, the very
similarity of these rates makes it meaningful to compare constant multipliers in
1991] DATA-BASED BANDWIDTH SELECTION METHOD 687
MSREs. Denoting these by 11;, i = 1,2, corresponding to h;, i = 1,2, respectively (h2 =
h2P or h2 = h2S here only), numerical calculation in the special case/is normal (and
also taking K and L to be normal) gives 112 :::::: 0.503111 in the case Q= 1. (Getting g
'wrong' increases 111> although discarding the scale model trick for a further kernel
estimation step would retain 111 as the appropriate constant). Such an improvement in
the constant augurs well for h2 ' s performance, too.
Reverting to diagonals-in estimates So instead of SNO was originally motivated by
algorithmic considerations in numerically solving equations like equation (4) for Park
and Marron's (1990) h selector. When using SNO' with a normal L, it was typical for
SNO < 0 for some small values of h. When SNO changes sign, the function in the

Downloaded from [Link] by guest on 17 May 2026

equation being zeroed has a discontinuity and changes sign. This is particularly
troublesome to many root finding procedures and can mislead them into reporting
this discontinuity as a solution. Such problems do not exist if we use So, since it has the
great advantage that it is necessarily always positive. Plots of the relevant functions
based on SNO (full curve) and So (broken curve) are given for a typical sample of size
n = 100 from the standard normal distribution in Fig. 1; the figure illustrates this
computational advantage of our proposal well. Moreover, although the h2S function
appears always to have the same shape as this one, in many cases the Park-Marron
function is even more badly behaved than this, with more (discontinuous) zero
crossings.

0
<D
ci

"ci

0
N
ci

0
ci

0
N
9

0
-e-
9

0
<D
9

0
co
9
0.0 0.2 0.4 0.6 0.8
Bandwidth, h

Fig. 1. Functions that must be zeroed: - - , for Park and Marron's (1990) method; --------, for fi2S
688 SHEATHER AND JONES [No.3,
4. SIMULAnON RESULTS
The simulation study reported here is rather brief. We compare all three hz
proposals of Section 3 with each other and with Park and Marron's (1990) procedure,
hi' We took the options of setting both K and L to be c/J, the standard normal density,
of using the interquartile range as Awherever needed and of setting gl(x) = 1.349 x
c/J (1.349x). For those hs not given by a direct formula, we used the Newton-Raphson
method to solve the required equations numerically.
We generated W = 500 realizations of data sets of size n = 50 and n = 100 from
each of four test densities. These densities are c/J, the normal mean mixture fl(x)
t t

Downloaded from [Link] by guest on 17 May 2026

= tc/J(x+ t) + c/J(x- t) and the normal variance mixtures fz(x) = c/J(x) +
t .J1Oc/J(.JlOx)andf3(x) = t c/J(x) + 5 c/J(lOx).
In Fig. 2 are shown approximate confidence intervals for the quantities R M =
100 if {MISE(h)/MISE(h o) - I}, where h stands for any of the bandwidths of
interest. MISE(h) denotes plugging the value of h into the MISE formula for fixed h
and its use is best justified by observing that R M is essentially proportional to
if (h/h o _1)z. It is fairly widely accepted that assessing the worth of any lfiis sensibly
done by comparing h with its 'target' ho (see Jones (1991) for reasons). R M , however,
gives us a feel for the practical worthwhileness of any bandwidth differences. Our
approximate confidence intervals are precisely the pivoted confidence intervals of
Park and Marron (1990) which we denote by (RM / (1 + V), RM / (1 - V). Here, RM is
the obvious estimate of R M obtained from the simulations and V = 1.96(2/ W)112 =
0.124.
There are some interesting conclusions to be drawn from Fig. 2. Recall first that hi
is considered to exhibit good performance by the current standards of the literature.
Away from the normal distribution, however, hi is often considerably worse than the
hzs. Of the hzs, hzss is never dominated by hZM ; hzp does better than hzs for the

A
hI
I H l--I
£ ~
l--I I--

A
h2P H H H ~

A
h25 l--I H f-l l--I

A
h2M l--I f-l 1--1 1--1

I
0 5 10 15 20 25 0 5 10 1520250 5 10 15 2025 0 5 10 15 2025
A
hI H 1--1 "-100 H f-l

pH H H H

A
h25 H H H H

A
h2M H H H H
I

Fig. 2. Approximate 95010 confidence intervals for R M for each of four data-based bandwidths, for two
sample sizes from each of four different underlying distributions (definitions of all these quantities are in
the text)
1991] DATA-BASED BANDWIDTH SELECTION METHOD 689
normal and (normal-like)fz distributions, worse for fl and is considerably inferior to
hzs for data from f3. On this evidence, the consistently good performance of hzs
suggests it as the method of choice. Further, on J. S. Marron's extensive simulation
evidence, hzs continues to perform very well, and better than other methods tried,
over a wide range of smooth density shapes.

5. THE Ii OF CHOICE
In Section 4, we saw that the bandwidth selection procedure resulting in hzs has

Downloaded from [Link] by guest on 17 May 2026

much to recommend it. For the reader's convenience, here we give full details of the
formulae for this bandwidth selector.
The bandwidth hzs , for use in equation (1), is the solution to the equation
(12)
where
n n
SD(a) = {n(n-I)}-l a - 5 ~ ~ <f>iv{a-I(X;-Xj)}.
;=1 j=1

From manipulation of equation (9), we have

otz{h) = 1.357{SD(a)/TD(b)}1I7h517.

Here, the constant is D 1(<f»/ R 117(<f» and the term in brackets is the second stage
estimate of R(f")/R(flll);
n n
TD(b) = - {n(n -I)} -lb- 7 ~ ~ <f>vi{b-I(X;-X)}.
;=1 j=1

For this estimate the bandwidths a and b are given by a normal scale model estimate of
equation (9) and of the corresponding formula for estimating R(flll) in Jones and
Sheather (1991) respectively to be
a = O.920An- 1I7 and b = O.9I2An- 1I9 ,
where Ais the sample interquartile range.
We successfully use the Newton-Raphson method to solve equation (12). A
Fortran subroutine is available on request from the first author.

ACKNOWLEDGEMENTS
The authors are very grateful to Peter Hall and especially Steve Marron for many
helpful comments and discussions. We also acknowledge Rob Hyndman for
programming assistance in an earlier version of this work. The editorial process has
been greatly beneficial in improving the standard of presentation of this material.
Some of M. C. Jones's work was supported by a Mathematical Sciences Research
Centre Visiting Fellowship at the Australian National University, Canberra,
Australia.
690 SHEATHER AND JONES [No.3,
REFERENCES
Hall, P. (1980) Objective methods for the estimation of window size in the nonparametric estimation of a
density. Unpublished.
Hall, P. and Marron, J. S. (1987) Estimation of integrated squared density derivatives. Statist. Probab.
Lett., 6, 109-115.
Jones, M. C. (1991) The roles of ISE and MISE in density estimation. Statist. Probab. Lett., to be
published.
Jones, M. C. and Sheather, S. J. (1991) Using non-stochastic terms to advantage in kernel-based
estimation of integrated squared density derivatives. Statist. Probab. Lett., to be published.
Park, B. U. and Marron, J. S. (1990) Comparison of data-driven bandwidth selectors. J. Am. Statist.

Downloaded from [Link] by guest on 17 May 2026

Ass., 85, 66-72.
Sheather, S. J. (1983) A data-based algorithm for choosing the window width when estimating the
density at a point. Comput. Statist. Data Anal., 1, 229-238.
--(1986) An improved data-based algorithm for choosing the window width when estimating the
density at a point. Comput. Statist. Data Anal., 4, 61-65.
Silverman, B. W. (1986) Density Estimation/or Statistics and Data Analysis. London: Chapman and
Hall.