BOOT.density function

Bootstrap bandwidth for a spatial kernel density estimate

Isotropic fixed or global (for adaptive) bandwidth selection for a standalone 2D density based on bootstrap estimation of the MISE.

BOOT.density(pp, hlim = NULL, eta = NULL, type = c("fixed", "adaptive"),
  hp = NULL, edge = c("uniform", "none"), ref.density = NULL,
  resolution = 64, rmdiag = TRUE, sim.adapt = list(N = 50, B = 100,
  dimz = 64, objective = FALSE), parallelise = NA, verbose = TRUE, ...)

Arguments

• pp: An object of class ppp giving the observed 2D data to be smoothed.
• hlim: An optional vector of length 2 giving the limits of the optimisation routine with respect to the bandwidth. If NULL, the function attempts to choose this automatically.
• eta: Fixed scalar bandwidth to use for the reference density estimate; if NULL it is calculated as the oversmoothing bandwidth of pp using OS. Ignored if ref.density is supplied. See `Details'.
• type: A character string indicating selection type. Either "fixed" (default) for selection of a constant bandwidth for the fixed-bandwidth estimator based on theory extended from results in Taylor (1989); or "adaptive" for selection of the global bandwidth for an adaptive kernel density. See `Details'.
• hp: Pilot bandwidth used for adaptive estimates in the bootstrap; see the argument of the same tag in bivariate.density. Ignored when type = "fixed" or when ref.density is supplied.
• edge: Character string dictating edge correction for the bootstrapped estimates. "uniform" (default) corrects based on evaluation grid coordinate. Setting edge="none" requests no edge correction.
• ref.density: Optional. An object of class bivden giving the reference density from which data will be generated. Based on theory, this must be a fixed-bandwidth estimate if type = "fixed"; see Details'. Must be edge-corrected if edge = "uniform"`.
• resolution: Spatial grid size; the optimisation will be based on a [resolution $x$ resolution] density estimate.
• rmdiag: Logical control value for removal of mirrored evaluation points as suggested by Taylor (1989) in the theoretical expression of the fixed-bandwidth MISE estimate. See Details'. Ignored when type = "adaptive"`
• sim.adapt: List of control values for bootstrap simulation in the adaptive case; see Details'. Ignored when type = "fixed"`.
• parallelise: Optional numeric argument to reduce computation time by invoking parallel processing, by giving the number of CPU cores to use in either evaluation (fixed) or in the actual bootstrap replicate generation (adaptive). Experimental. Test your system first using parallel::detectCores() to identify the number of cores available to you.
• verbose: Logical value indicating whether to print function progress during execution.
• ...: Optional arguments controlling scaling to be passed to multiscale.density for the adaptive bootstrap; ignored when type = "fixed".

Details

For a 2D kernel density estimate $\hat{f}$ defined on $W \in R^2$, the mean integrated squared error (MISE) is given by $E[\int_W (\hat{f}(x) - f(x))^2 dx]$, where $f$ is the corresponding true density. Given an observed data set $X$ (argument pp) of $n$ observations, this function finds the bandwidth $h$ that minimises

where $\hat{f}(x)$ is a density estimate of $X$ constructed with reference' bandwidth $\eta$ (argument etaorref.density`), and $\hat{f}^*(x)$ is a density estimate using bandwidth $h$ of $n$ observations $X^*$ generated from $\hat{f}(x)$. The notation $E^*$ denotes expectation with respect to the distribution of the $X^*$.

• Fixed: When type = "fixed", the function assumes you want to select a constant bandwidth for use with the fixed-bandwith density estimator. This implementation is based on extending the remarkable results of Taylor (1989) (see also Sain et al., 1994), who demonstrates that when the Gaussian kernel is being used, we can find the optimal $h$ with respect to the aforementioned bootstrap-estimated MISE without any actual resampling. This implementation extends these results to the bivariate setting, and allows for edge-correction of both the reference and bootstrap densities.

    * Taylor (1989) does not distinguish between the reference bandwidth $\eta$ and the target of optimisation, $h$, thus allowing the reference bandwidth to vary alongside the target in the optimisation. This is not optimal, and this function always assumes a static reference bandwidth. Hall et al. (1992) indicate that a generous amount of smoothing is to be preferred in the reference density (hence the default `eta` set using `OS`).
    * If `ref.density` is supplied, it must be a fixed-bandwidth density estimate as an object of class `bivden` for validity of the theory. Edge-correction must be present if `edge = "uniform"`; and it must be evaluated on the same spatial domain as dictated by `Window(pp)` and `resolution`. If unsupplied, the function internally computes an appropriate fixed-bandwidth density estimate using `eta` as the reference bandwidth.
    * Finally, Taylor (1989) argues it is preferable to avoid summation at identical evaluation grid points in the expression for the optimal bandwidth, which is performed when `rmdiag = TRUE`. Setting `rmdiag = FALSE` disables this correction.
  

• Adaptive: When type = "adaptive", the function assumes you want to select a global bandwidth (argument h0 in bivariate.density) for use in 2D adaptive kernel density estimation.

    * An expression similar to Taylor (1989) is not possible for the adaptive estimator. Thus, in the adaptive setting, the optimal bootstrap bandwidth is calculated by brute force as was performed in Davies and Baddeley (2018) by taking advantage of the multiscale estimation theory implemented in `multiscale.density`. The value that minimises an interpolating cubic spline of the estimated MISE on bandwidth is identified as the optimal global bandwidth.
    * The user can pass either a fixed or adaptive `bivden` object to `ref.density`. If this is the case, `hp` is ignored and the pilot bandwidth for each iteration of the bootstrap in estimation of the $\hat{f}^*(x)$ uses `ref.density$hp` (if `ref.density` is adaptive) or `ref.density$h0` (if `ref.density` is fixed). When `ref.density` is unsupplied, the function uses a fixed-bandwidth kernel estimate with bandwidth `eta` as the reference density, and if additionally `hp` is unsupplied, the same value `eta` is used for the constant pilot bandwidth.
    * Control over the bootstrap is achieved with four optional named arguments passed as a list to `sim.adapt`. `N` controls the number of bootstrap iterates per bandwidth; `B` controls the resolution of the sequence of bandwidths trialled (i.e. between `hlim[1]` and `hlim[2]`); `dimz` specifies the resolution of the bandwidth axis in the trivariate convolution evaluated by `multiscale.density`; and `objective` specifies whether to return the set of estimated MISEs for all bandwidths (nice to plot), or merely the optimal bandwidth (see `Value').
    * The `...` are intended for any relevant optional arguments to be passed to the internal call to `multiscale.density`, such as `gamma.scale` or `trim`.
  

Returns

The optimal fixed or global (for adaptive) scalar bandwidth. If simargs$objective = TRUE for the adaptive bootstrap, the return object is instead a $[$simargs$B $x 2]$ matrix, with the first column giving the trialled bandwidth and the second giving the corresponding value of the estimated bootstrap MISE.

References

Davies, T.M. and Baddeley A. (2018), Fast computation of spatially adaptive kernel estimates, Statistics and Computing, 28 (4), 937-956.

Hall, P., Marron, J.S. and Park, B.U. (1992) Smoothed cross-validation, Probability Theory and Related Fields, 92 , 1-20.

Sain, S.R., Baggerly, K.A. and Scott, D.W. (1994) Cross-validation of multivariate densities, Journal of the American Statistical Association, 89 , 807-817.

Taylor, C.C. (1989) Bootstrap choice of the smoothing parameter in kernel density estimation, Biometrika, 76 , 705-712.

Author(s)

T.M. Davies

Warning

Even with the implemented computational tricks, bootstrapping for bandwidth selection for spatial data is still computationally demanding, especially for adaptive kernel estimates. The user can reduce this time by keeping the evaluation grid at modest resolutions, and experimenting with parallelising the internal loops via parallelise. The `Examples' section offers some rough indications of evaluation times on this author's local machine.

See Also

bivariate.density, OS, multiscale.density

Examples

data(pbc)

## Fixed bandwidth selection ##
BOOT.density(pbc) # ~20 secs
BOOT.density(pbc,eta=OS(pbc)/2) # halve default reference bandwidth
BOOT.density(pbc,eta=OS(pbc)*2) # double default reference bandwidth

# supplying pre-defined reference density as fixed-bandwidth 'bivden' object
pbcfix <- bivariate.density(pbc,h0=2.5,resolution=64)
system.time(hfix <- BOOT.density(pbc,ref.density=pbcfix))
hfix

## Global (for adaptive) bandwidth selection ##
# ~200 secs next line; use 'parallelise' for speedup
system.time(hada <- BOOT.density(pbc,type="adaptive")) # minimal usage for adaptive bootstrap
hada

# ~80 secs next line. Set custom h limits; increase reference bandwidth;
#    set custom pilot bandwidth; return objective function; use 'parallelise' for speedup
system.time(hada <- BOOT.density(pbc,hlim=c(0.9,8),eta=3.5,type="adaptive",
                                 hp=OS(pbc)/2,sim.adapt=list(objective=TRUE)))
hada[which.min(hada[,2]),1]
plot(hada);abline(v=hada[which.min(hada[,2]),1],col=2)