2.4 Bandwidth selection
As we saw in the previous sections, the kde critically depends on the bandwidth employed. The purpose of this section is to introduce objective and automatic bandwidth selectors that attempt to minimize the estimation error of the target density
The first step is to define a global, rather than local, error criterion. The Integrated Squared Error (ISE),
is the squared distance between the kde and the target density. The ISE is a random quantity, since it depends directly on the sample As a consequence, looking for an optimal-ISE bandwidth is a hard task, since the optimality is dependent on the sample itself and not only on the population and To avoid this problem, it is usual to compute the Mean Integrated Squared Error (MISE):
The MISE is convenient due to its mathematical tractability and its natural relation to the MSE. There are, however, other error criteria that present attractive properties, such as the Mean Integrated Absolute Error (MIAE):
The MIAE, unlike the MISE, has the appeal of being invariant with respect to monotone transformations of the density. For example, if is the density of and then the change of variables37 gives
Despite this attractive invariance property, the analysis of MIAE is substantially more complicated than the MISE. We refer to Devroye and Györfi (1985) for a comprehensive treatment of absolute value metrics for kde.
Once the MISE is set as the error criterion to be minimized, our aim is to find
For that purpose, we need an explicit expression of the MISE which we can attempt to minimize. The following asymptotic expansion for the MISE solves this issue.
Corollary 2.2 Under A1–A3,
Therefore, when
Exercise 2.13 Prove Corollary 2.2 by integrating the MSE of
The dominating part of the MISE is denoted by AMISE, which stands for Asymptotic MISE:
The closed-form expression of the AMISE allows obtaining a bandwidth that minimizes this error.
Corollary 2.3 The bandwidth that minimizes the AMISE is
The optimal AMISE is:
Exercise 2.14 Prove Corollary 2.3 by solving
The AMISE-optimal order deserves some further inspection. It can be seen in Section 3.2 in Scott (2015) that the AMISE-optimal order for the histogram of Section 2.1.1 (not the moving histogram) is Two aspects are interesting when comparing this result with (2.24):
- First, the MISE of the histogram is asymptotically larger than the MISE of the kde.38 This is a quantification of the quite apparent visual improvement of the kde over the histogram.
- Second, appears instead of evidencing that the histogram is affected by how fast varies and not only by the curvature of the target density
Unfortunately, the AMISE bandwidth depends on which measures the curvature of the density. As a consequence, it can not be readily applied in practice, as is unknown! In the next subsection we will see how to plug-in estimates for
2.4.1 Plug-in rules
A simple solution to estimate is to assume that is the density of a and then plug-in the form of the curvature for such density,39
While doing so, we approximate the curvature of an arbitrary density by means of the curvature of a normal and we have that
Interestingly, the bandwidth is directly proportional to the standard deviation of the target density. Replacing by an estimate yields the normal scale bandwidth selector, which we denote to emphasize its randomness:
The estimate can be chosen as the standard deviation or, in order to avoid the effects of potential outliers, as the standardized interquantile range
or as
When combined with a normal kernel, for which and this particularization of featuring (2.25) gives the famous rule-of-thumb for bandwidth selection:
is implemented in R through the function bw.nrd
.40
# Data
set.seed(667478)
n <- 100
x <- rnorm(n)
# Rule-of-thumb
bw.nrd(x = x)
## [1] 0.4040319
# bwd.nrd employs 1.34 as an approximation for diff(qnorm(c(0.25, 0.75)))
# Same as
iqr <- diff(quantile(x, c(0.25, 0.75))) / diff(qnorm(c(0.25, 0.75)))
1.06 * n^(-1/5) * min(sd(x), iqr)
## [1] 0.4040319
The previous selector is an example of a zero-stage plug-in selector, a terminology inspired by the fact that the scalar was estimated by plugging a parametric assumption directly, without attempting a nonparametric estimation first. Another possibility could have been to estimate nonparametrically and then to plug-in the estimate into Let’s explore this possibility in more detail next.
First, note the useful equality
This equality follows by an iterative application of integration by parts. For example, for take and It gives
under the assumption that the derivatives vanish at infinity. Applying again integration by parts with and gives the result. This simple derivation has an important consequence: for estimating the functionals it suffices to estimate the functionals
for In particular,
Thanks to the expression (2.26), a possible way to estimate nonparametrically is
where is the -th derivative of a kde with bandwidth and kernel i.e.,
Note that and can be different from and respectively. It turns out that estimating involves the adequate selection of a bandwidth The agenda is analogous to the one for but now taking into account that both and are scalar quantities:
Under certain regularity assumptions,41 the asymptotic bias and variance of are obtained. With them, we can compute the asymptotic expansion of the MSE42 and obtain the Asymptotic Mean Squared Error AMSE:
Note: is the highest integer such that In these notes we have restricted the exposition to the case for the kernels but there are theoretical gains if one allows high-order kernels with vanishing even moments larger than for estimating
Obtain the AMSE-optimal bandwidth:
The order of the optimal AMSE is
which shows that a parametric-like rate of convergence can be achieved with high-order kernels. If we consider and then
The result above has a major problem: it depends on ! Thus, if we want to estimate by we will need to estimate But will depend on and so on! The solution to this convoluted problem is to stop estimating the functional after a given number of stages, hence the terminology -stage plug-in selector. At the -th stage, the functional inside the AMSE-optimal bandwidth for estimating nonparametrically43 is computed assuming that the density is a 44 for which
Typically, two stages are considered a good trade-off between bias (mitigated when increases) and variance (increases with ) of the plug-in selector. This is the method proposed by Sheather and Jones (1991), where they consider and yielding what we call the Direct Plug-In (DPI). The algorithm is:
Estimate using where is given in (2.25).
Estimate using from (2.27), where
Estimate using from (2.27), where
The selected bandwidth is
Remark. The derivatives for the normal kernel can be obtained using the (probabilists’) Hermite polynomials: For these are: and Hermite polynomials satisfy the recurrence relation for
is implemented in R through the function bw.SJ
(use method = "dpi"
). An alternative and faster implementation is ks::hpi
, which also accounts for more flexibility and has a somehow more complete documentation.
# Data
set.seed(672641)
x <- rnorm(100)
# DPI selector
bw.SJ(x = x, method = "dpi")
## [1] 0.5006905
# Similar to
ks::hpi(x) # Default is two-stage
## [1] 0.4999456
Exercise 2.15 Apply and inspect the kde of airquality$Ozone
using the DPI selector (both bw.SJ
and ks::hpi
). What can you conclude from the estimate?
2.4.2 Cross-validation
We now turn our attention to a different philosophy of bandwidth estimation. Instead of trying to minimize the AMISE by plugging estimates for the unknown curvature term, we directly attempt to minimize the MISE. The idea is to use the sample twice: one for computing the kde and other for evaluating its performance on estimating To avoid the clear dependence on the sample, we do this evaluation in a cross-validatory way: the data used for computing the kde is not used for its evaluation.
We begin by expanding the square in the MISE expression:
Since the last term does not depend on minimizing is equivalent to minimizing
This quantity is unknown, but it can be estimated unbiasedly by
where is the leave-one-out kde and is based on the sample with the removed:
Exercise 2.17 (Exercise 3.3 in Wand and Jones (1995)) Prove that
The motivation for (2.29) is the following. The first term is unbiased by design.45 The second arises from approximating by Monte Carlo from the sample in other words, by replacing with This gives
and, in order to mitigate the dependence of the sample, we replace with above. In that way, we use the sample for estimating the integral involving but for each we compute the kde on the remaining points.
The Least Squares Cross-Validation (LSCV) selector, also denoted Unbiased Cross-Validation (UCV) selector, is defined as
Numerical optimization is required for obtaining contrary to the previous plug-in selectors, and there is little control on the shape of the objective function. This will also be the case for the subsequent bandwidth selectors. The following remark warns about the dangers of numerical optimization in this context.
Remark. Numerical optimization of the LSCV function can be challenging. In practice, several local minima are possible, and the roughness of the objective function can vary notably depending on and As a consequence, optimization routines may get trapped in spurious solutions. To be on the safe side, it is always advisable to check (when possible) the solution by plotting for a range of or to perform a search in a bandwidth grid:
is implemented in R through the function bw.ucv
. This function uses R’s optimize
, which is quite sensitive to the selection of the search interval.46 Therefore, some care is needed; that is why the bw.ucv.mod
function is presented below.
# Data
set.seed(123456)
x <- rnorm(100)
# UCV gives a warning
bw.ucv(x = x)
## [1] 0.4499177
# Extend search interval
bw.ucv(x = x, lower = 0.01, upper = 1)
## [1] 0.5482419
# bw.ucv.mod replaces the optimization routine of bw.ucv with an exhaustive
# search on "h_grid" (chosen adaptatively from the sample) and optionally
# plots the LSCV curve with "plot_cv"
bw.ucv.mod <- function(x, nb = 1000L,
h_grid = 10^seq(-3, log10(1.2 * sd(x) *
length(x)^(-1/5)), l = 200),
plot_cv = FALSE) {
if ((n <- length(x)) < 2L)
stop("need at least 2 data points")
n <- as.integer(n)
if (is.na(n))
stop("invalid length(x)")
if (!is.numeric(x))
stop("invalid 'x'")
nb <- as.integer(nb)
if (is.na(nb) || nb <= 0L)
stop("invalid 'nb'")
storage.mode(x) <- "double"
hmax <- 1.144 * sqrt(var(x)) * n^(-1/5)
Z <- .Call(stats:::C_bw_den, nb, x)
d <- Z[[1L]]
cnt <- Z[[2L]]
fucv <- function(h) .Call(stats:::C_bw_ucv, n, d, cnt, h)
## Original
# h <- optimize(fucv, c(lower, upper), tol = tol)$minimum
# if (h < lower + tol | h > upper - tol)
# warning("minimum occurred at one end of the range")
## Modification
obj <- sapply(h_grid, function(h) fucv(h))
h <- h_grid[which.min(obj)]
if (h %in% range(h_grid))
warning("minimum occurred at one end of h_grid")
if (plot_cv) {
plot(h_grid, obj, type = "o")
rug(h_grid)
abline(v = h, col = 2, lwd = 2)
}
h
}
# Compute the bandwidth and plot the LSCV curve
bw.ucv.mod(x = x, plot_cv = TRUE, h_grid = 10^seq(-1.25, 0.5, l = 200))
## [1] 0.5431561
# We can compare with the default bw.ucv output
abline(v = bw.ucv(x = x), col = 3)

Figure 2.9: LSCV curve evaluated for a grid of bandwidths, with a clear global minimum corresponding to
The next cross-validation selector is based on Biased Cross-Validation (BCV). The BCV selector presents a hybrid strategy that combines plug-in and cross-validation ideas. It starts by considering the AMISE expression in (2.23)
and then plugs-in an estimate for based on a modification of The modification is
a leave-out-diagonals estimate of It is designed to reduce the bias of since (Scott and Terrell 1987). Plugging (2.30) into the AMISE expression yields the BCV objective function and the BCV bandwidth selector:
where stands for the smallest local minimizer of The consideration of the local minimum is because, by design, as for fixed : at a faster rate than Therefore, when minimizing some care is required, as one is actually interested in obtaining the smallest local minimizer. Consequently, bandwidth grids with an upper extreme that is too large47 are to be avoided, as these will miss the local minimum in favor of the global one at
The most appealing property of is that it has a considerably smaller variance than This reduction in variance comes at the price of an increased bias, which tends to make larger than
is implemented in R through the function bw.bcv
. Again, bw.bcv
uses R’s optimize
so the bw.bcv.mod
function is presented to have better guarantees on finding the first local minimum. Quite some care is needed with the range of bandwidth grid, though, to avoid the global minimum for large bandwidths.
# Data
set.seed(123456)
x <- rnorm(100)
# BCV gives a warning
bw.bcv(x = x)
## [1] 0.4500924
# Extend search interval
args(bw.bcv)
## function (x, nb = 1000L, lower = 0.1 * hmax, upper = hmax, tol = 0.1 *
## lower)
## NULL
bw.bcv(x = x, lower = 0.01, upper = 1)
## [1] 0.5070129
# bw.bcv.mod replaces the optimization routine of bw.bcv with an exhaustive
# search on "h_grid" (chosen adaptatively from the sample) and optionally
# plots the BCV curve with "plot_cv"
bw.bcv.mod <- function(x, nb = 1000L,
h_grid = 10^seq(-3, log10(1.2 * sd(x) *
length(x)^(-1/5)), l = 200),
plot_cv = FALSE) {
if ((n <- length(x)) < 2L)
stop("need at least 2 data points")
n <- as.integer(n)
if (is.na(n))
stop("invalid length(x)")
if (!is.numeric(x))
stop("invalid 'x'")
nb <- as.integer(nb)
if (is.na(nb) || nb <= 0L)
stop("invalid 'nb'")
storage.mode(x) <- "double"
hmax <- 1.144 * sqrt(var(x)) * n^(-1/5)
Z <- .Call(stats:::C_bw_den, nb, x)
d <- Z[[1L]]
cnt <- Z[[2L]]
fbcv <- function(h) .Call(stats:::C_bw_bcv, n, d, cnt, h)
## Original code
# h <- optimize(fbcv, c(lower, upper), tol = tol)$minimum
# if (h < lower + tol | h > upper - tol)
# warning("minimum occurred at one end of the range")
## Modification
obj <- sapply(h_grid, function(h) fbcv(h))
h <- h_grid[which.min(obj)]
if (h %in% range(h_grid))
warning("minimum occurred at one end of h_grid")
if (plot_cv) {
plot(h_grid, obj, type = "o")
rug(h_grid)
abline(v = h, col = 2, lwd = 2)
}
h
}
# Compute the bandwidth and plot the BCV curve
bw.bcv.mod(x = x, plot_cv = TRUE, h_grid = 10^seq(-1.25, 0.5, l = 200))
## [1] 0.5111433
# We can compare with the default bw.bcv output
abline(v = bw.bcv(x = x), col = 3)

Figure 2.10: BCV curve evaluated for a grid of bandwidths, with a clear local minimum corresponding to Observe the decreasing pattern when
2.4.3 Comparison of bandwidth selectors
Next, we state some insights from the convergence results of the DPI, LSCV, and BCV selectors. All of them are based on results of the kind
where depends on and only, and measures how variable the selector is. The rate serves to quantify how fast48 the relative error decreases (the larger the the faster the convergence).
Under certain regularity conditions, we have:
and Both cross-validation selectors have a slow rate of convergence.49 Inspection of the variances of both selectors reveals that, for the normal kernel, Therefore, LSCV is considerably more variable than BCV.
Thus, the DPI selector has a convergence rate much faster than the cross-validation selectors. There is an appealing explanation for this phenomenon. Recall that minimizes the slightly modified version of given by
and
is a leave-out-diagonals estimate of Despite being different from it serves for building a DPI analogous to BCV that points towards the precise fact that drags down the performance of BCV. The modified version of the DPI minimizes
where is independent of The two methods differ in the way is chosen: BCV sets and the modified DPI looks for the best in terms of the It can be seen that whereas the used in BCV is asymptotically This suboptimality on the choice of is the reason of the asymptotic deficiency of BCV.
We now focus on exploring the empirical performance of bandwidth selectors. The workhorse for doing that is simulation. A popular collection of simulation scenarios was given by Marron and Wand (1992) and is conveniently available through the package nor1mix
. This collection is formed by normal -mixtures of the form
where and Densities of this form are especially attractive since they allow for arbitrary flexibility and, if the normal kernel is employed, they allow for explicit and exact MISE expressions:
These exact MISE expressions are highly convenient for computing as defined in (2.22), by minimizing (2.33) numerically.
# Available models
?nor1mix::MarronWand
# Simulating
samp <- nor1mix::rnorMix(n = 500, obj = nor1mix::MW.nm9)
# MW object in the second argument
hist(samp, freq = FALSE)
# Density evaluation
x <- seq(-4, 4, length.out = 400)
lines(x, nor1mix::dnorMix(x = x, obj = nor1mix::MW.nm9), col = 2)
# Plot a MW object directly
# A normal with the same mean and variance is plotted in dashed lines --
# you can remove it with argument p.norm = FALSE
par(mfrow = c(2, 2))
plot(nor1mix::MW.nm5)
plot(nor1mix::MW.nm7)
plot(nor1mix::MW.nm10)
plot(nor1mix::MW.nm12)
lines(nor1mix::MW.nm7, col = 2) # Also possible
Exercise 2.18 Implement the using (2.22) and (2.33) for model nor1mix::MW.nm6
. Then, compute by Monte Carlo the densities of and Compare them for Describe in detail the results and the major takeaways in relation to the properties of the bandwidth selectors.
Exercise 2.19 Compare the MISE and AMISE criteria in three densities in nor1mix
of your choice. To that purpose, code (2.33) and the AMISE expression for the normal kernel, and compare the two error curves. Compare them for adding a vertical line to represent the and bandwidths. Describe in detail the results and the major takeaways in relation to using as a proxy for .
A key practical issue that emerges after discussing several bandwidth selectors is the following:
Which bandwidth selector is the most adequate for a given dataset?
Unfortunately, there is no simple and universal answer to this question. There are, however, a series of useful facts and suggestions:
- Trying several selectors and inspecting the results may help to determine which one is estimating the density better.
- The DPI selector has a convergence rate much faster than the cross-validation selectors. Therefore, in theory it is expected to perform better than LSCV and BCV. For this reason, it tends to be among the preferred bandwidth selectors in the literature.
- Cross-validatory selectors may be better suited for highly non-normal and rough densities, in which plug-in selectors may end up oversmoothing.
- LSCV tends to be considerably more variable than BCV.
- The RT is a quick, simple, and inexpensive selector. However, it tends to give bandwidths that are too large for non-normal-like data.
Figure 2.11 presents a visualization of the performance of the kde with different bandwidth selectors, carried out in the family of mixtures by Marron and Wand (1992).
Figure 2.11: Performance comparison of bandwidth selectors. The RT, DPI, LSCV, and BCV are computed for each sample for a normal mixture density. For each sample, the animation computes the ISEs of the selectors and sorts them from best to worst. Changing the scenarios gives insight into the adequacy of each selector to hard- and simple-to-estimate densities. Application available here.
References
For defining we consider the transformed kde seen in Section 2.5.1.↩︎
Since ↩︎
We only use a parametric assumption for estimating the curvature of in present not for directly estimating itself.↩︎
Not to confuse with
bw.nrd0
!↩︎See Section 3.5 in Wand and Jones (1995) for full details.↩︎
Recall there is no “I”, since we are estimating a scalar, not a function.↩︎
For the rule-of-thumb selector, and we estimate parametrically. If we estimate parametrically, which is then required for estimating nonparametrically.↩︎
Hence, what we are doing is “sweeping under nonparametric carpets” this parametric estimate on the characteristic of the density, ultimately required for computing the bandwidth selector.↩︎
Why so? Because is unbiased for estimating ↩︎
Long intervals containing the solution may lead to unsatisfactory termination of the search; short intervals might not contain the minimum.↩︎
The precise point at which a bandwidth is “too big” can be formalized with the maximal smoothing principle, as elaborated in Section 3.2.2 in Wand and Jones (1995).↩︎
Recall that another way of writing (2.31) is as the relative error being asymptotically distributed as a ↩︎
Compare it with the rate of the CLT!↩︎