BOOT.spattemp.RdBandwidth selection for standalone spatiotemporal density/intensity based on bootstrap estimation of the MISE, providing an isotropic scalar spatial bandwidth and a scalar temporal bandwidth.
BOOT.spattemp(pp, tt = NULL, tlim = NULL, eta = NULL, nu = NULL,
sedge = c("uniform", "none"), tedge = sedge, ref.density = NULL,
sres = 64, tres = sres, start = NULL, verbose = TRUE)An object of class ppp giving the spatial coordinates of the observations to be smoothed. Possibly marked with the time of each event; see argument tt.
A numeric vector of equal length to the number of points in pp, giving the time corresponding to each spatial observation. If unsupplied, the function attempts to use the values in the marks attribute of the ppp.object in pp.
A numeric vector of length 2 giving the limits of the temporal domain over which to smooth. If supplied, all times in tt must fall within this interval (equality with limits allowed). If unsupplied, the function simply uses the range of the observed temporal values.
Fixed scalar bandwidth to use for the spatial margin of 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'.
Fixed scalar bandwidth to use for the temporal margin of the reference density estimate; if NULL it is calculated from tt using the univariate version of Terrell's (1990) oversmoothing principle. Ignored if ref.density is supplied. See `Details'.
Character string dictating spatial edge correction. "uniform" (default) corrects based on evaluation grid coordinate. Setting sedge="none" requests no edge correction.
As sedge, for temporal edge correction.
Optional. An object of class stden giving the reference density from which data is assumed to originate in the bootstrap. Must be spatially edge-corrected if sedge = "uniform".
Numeric value > 0. Resolution of the [sres \(\times\) sres] evaluation grid in the spatial margin.
Numeric value > 0. Resolution of the evaluation points in the temporal margin as defined by the tlim interval. If unsupplied, the density is evaluated at integer values between tlim[1] and tlim[2].
Optional positive numeric vector of length 2 giving starting values for the internal call to optim, in the order of (<spatial bandwidth>, <temporal bandwidth>).
Logical value indicating whether to print a function progress bar to the console during evaluation.
For a spatiotemporal kernel density estimate \(\hat{f}\) defined on \(W x T \in R^3\), the mean integrated squared error (MISE) is given by \(E[\int_W \int_T (\hat{f}(x,t) - f(x,t))^2 dt dx]\), where \(f\) is the corresponding true density. Given observed spatiotemporal locations \(X\) (arguments pp and tt) of \(n\) observations, this function finds the scalar spatial bandwidth \(h\) and scalar temporal bandwidth \(\lambda\) that jointly minimise
$$E^*[\int_W \int_T (\hat{f}^*(x,t) - \hat{f}(x,t))^2 dt dx],$$
where \(\hat{f}(x,t)\) is a density estimate of \(X\) constructed with `reference' bandwidths \(\eta\) (spatial; argument eta) and \(\nu\) (temporal; argument nu); \(\hat{f}^*(x,t)\) is a density estimate using bandwidths \(h\) and \(\lambda\) of \(n\) observations \(X^*\) generated from \(\hat{f}(x,t)\). The notation \(E^*\) denotes expectation with respect to the distribution of the \(X^*\). The user may optionally supply ref.density as an object of class stden, which must be evaluated on the same spatial and temporal domains \(W\) and \(T\) as the data (arguments pp, tt, and tlim). In this case, the reference bandwidths are extracted from this object, and eta and nu are ignored.
This function is based on an extension of the theory of Taylor (1989) to the spatiotemporal domain and to cope with the inclusion of edge-correction factors. No resampling is necessary due to the theoretical properties of the Gaussian kernel.
A numeric vector of length 2 giving the jointly optimised spatial and temporal bandwidths (named h and lambda respectively).
Taylor, C.C. (1989) Bootstrap choice of the smoothing parameter in kernel density estimation, Biometrika, 76, 705-712.
Bootstrapping for spatiotemporal bandwidth selection for spatiotemporal data is very computationally demanding. Keeping verbose = TRUE offers an indication of the computational burden by printing each pair of bandwidths at each iteration of the optimisation routine. The `Examples' section also offers some rough indications of evaluation times on this author's local machine.
# \donttest{
data(burk) # Burkitt's Uganda lymphoma data
burkcas <- burk$cases
#~85 secs
hlam1 <- BOOT.spattemp(burkcas)
#> Initialising...Done.
#> Optimising...
#> h = 11.24394 ; lambda = 560.6947
#> h = 67.31342 ; lambda = 560.6947
#> h = 11.24394 ; lambda = 616.7642
#> h = 39.27868 ; lambda = 574.7121
#> h = 25.26131 ; lambda = 581.7208
#> h = 18.25263 ; lambda = 585.2251
#> h = 18.25263 ; lambda = 641.2946
#> h = 16.50046 ; lambda = 621.1447
#> h = 9.491773 ; lambda = 652.6837
#> h = 16.06241 ; lambda = 602.0898
#> h = 21.31893 ; lambda = 606.4702
#> h = 13.76269 ; lambda = 614.1907
#> h = 13.32465 ; lambda = 595.1359
#> h = 14.1186 ; lambda = 601.6381
#> h = 11.81888 ; lambda = 613.739
#> h = 15.00153 ; lambda = 605.0021
#> h = 12.87976 ; lambda = 610.8267
#> h = 14.47109 ; lambda = 606.4582
#> h = 14.11518 ; lambda = 619.0109
#> h = 14.11347 ; lambda = 627.6973
#> h = 14.82186 ; lambda = 619.9648
#> h = 14.55707 ; lambda = 618.5213
#> h = 14.19945 ; lambda = 639.7604
#> h = 14.06363 ; lambda = 656.4115
#> h = 13.62003 ; lambda = 665.5875
#> h = 13.1515 ; lambda = 689.1206
#> h = 13.10167 ; lambda = 717.8347
#> h = 12.59577 ; lambda = 762.9034
#> h = 11.68364 ; lambda = 795.6125
#> h = 13.46863 ; lambda = 691.2117
#> h = 12.9129 ; lambda = 764.9946
#> h = 12.7936 ; lambda = 802.9316
#> h = 11.92073 ; lambda = 874.6233
#> h = 12.30771 ; lambda = 828.7704
#> h = 12.50553 ; lambda = 868.7986
#> h = 12.46042 ; lambda = 921.7462
#> h = 12.94631 ; lambda = 895.9074
#> h = 13.2656 ; lambda = 929.4758
#> h = 12.93242 ; lambda = 1048.29
#> h = 13.00184 ; lambda = 1170.97
#> h = 13.73761 ; lambda = 1056.02
#> h = 13.41831 ; lambda = 1022.452
#> h = 13.08513 ; lambda = 1141.266
#> h = 13.13025 ; lambda = 1088.319
#> h = 12.64436 ; lambda = 1114.157
#> h = 12.25739 ; lambda = 1160.01
#> h = 12.44654 ; lambda = 1074.129
#> h = 12.61746 ; lambda = 1077.677
#> h = 12.3294 ; lambda = 1143.544
#> h = 12.78167 ; lambda = 1072.104
#> h = 12.80857 ; lambda = 1108.584
#> h = 12.76079 ; lambda = 1100.858
#> h = 12.62348 ; lambda = 1142.911
#> h = 12.74212 ; lambda = 1089.806
#> h = 12.66303 ; lambda = 1125.209
#> h = 12.72235 ; lambda = 1098.656
#> h = 12.60592 ; lambda = 1111.956
#> h = 12.72207 ; lambda = 1103.632
#> h = 12.80006 ; lambda = 1088.131
#> h = 12.68329 ; lambda = 1107.651
#> h = 12.68301 ; lambda = 1112.627
#> h = 12.69284 ; lambda = 1109.134
#> h = 12.65406 ; lambda = 1113.153
#> h = 12.70507 ; lambda = 1106.012
#> h = 12.69551 ; lambda = 1104.529
#> h = 12.69484 ; lambda = 1105.68
#> h = 12.67306 ; lambda = 1107.319
#> h = 12.68106 ; lambda = 1106.992
#> h = 12.66951 ; lambda = 1108.963
#> h = 12.68851 ; lambda = 1106.501
#> h = 12.68629 ; lambda = 1105.842
#> h = 12.68554 ; lambda = 1106.294
#> h = 12.69298 ; lambda = 1105.803
#> h = 12.68404 ; lambda = 1106.695
#> Done.
#~75 secs. Widen time limits, reduce ref. bw.
hlam2 <- BOOT.spattemp(burkcas,tlim=c(400,5800),eta=8,nu=450)
#> Initialising...Done.
#> Optimising...
#> h = 8 ; lambda = 450
#> h = 53 ; lambda = 450
#> h = 8 ; lambda = 495
#> h = 30.5 ; lambda = 461.25
#> h = 19.25 ; lambda = 466.875
#> h = 13.625 ; lambda = 469.6875
#> h = 13.625 ; lambda = 514.6875
#> h = 16.4375 ; lambda = 547.0312
#> h = 19.25 ; lambda = 489.375
#> h = 10.8125 ; lambda = 493.5938
#> h = 10.8125 ; lambda = 538.5938
#> h = 9.40625 ; lambda = 573.0469
#> h = 8 ; lambda = 517.5
#> h = 12.21875 ; lambda = 515.3906
#> h = 9.40625 ; lambda = 516.7969
#> h = 11.51562 ; lambda = 515.7422
#> h = 11.51562 ; lambda = 560.7422
#> h = 11.33984 ; lambda = 543.9551
#> h = 10.8125 ; lambda = 583.5938
#> h = 10.46094 ; lambda = 617.5195
#> h = 9.757812 ; lambda = 595.3711
#> h = 10.19727 ; lambda = 586.7139
#> h = 9.845703 ; lambda = 665.6396
#> h = 9.362305 ; lambda = 729.1626
#> h = 9.625977 ; lambda = 759.9683
#> h = 9.340332 ; lambda = 846.5955
#> h = 8.527344 ; lambda = 871.6113
#> h = 9.977539 ; lambda = 681.0425
#> h = 10.24121 ; lambda = 711.8481
#> h = 10.02148 ; lambda = 716.1768
#> h = 9.669922 ; lambda = 795.1025
#> h = 9.746826 ; lambda = 766.5875
#> h = 9.351318 ; lambda = 810.379
#> h = 9.853943 ; lambda = 739.7273
#> h = 9.51886 ; lambda = 786.8285
#> h = 9.770172 ; lambda = 751.5026
#> h = 9.891022 ; lambda = 758.1219
#> h = 9.692238 ; lambda = 759.5067
#> h = 9.715584 ; lambda = 744.4217
#> h = 9.723394 ; lambda = 749.9632
#> h = 9.801329 ; lambda = 741.9591
#> h = 9.719511 ; lambda = 755.1198
#> h = 9.672733 ; lambda = 753.5804
#> h = 9.745812 ; lambda = 752.022
#> h = 9.741928 ; lambda = 757.1786
#> h = 9.728028 ; lambda = 751.7671
#> h = 9.75433 ; lambda = 748.6693
#> h = 9.728215 ; lambda = 753.5072
#> h = 9.746 ; lambda = 753.7622
#> h = 9.732521 ; lambda = 752.2658
#> h = 9.750118 ; lambda = 750.7807
#> h = 9.733691 ; lambda = 752.8256
#> h = 9.7204 ; lambda = 753.0693
#> h = 9.739459 ; lambda = 752.2839
#> h = 9.740629 ; lambda = 752.8436
#> h = 9.738602 ; lambda = 752.6992
#> Done.
#~150 secs. Increase ref. bw., custom starting vals
hlam3 <- BOOT.spattemp(burkcas,eta=20,nu=800,start=c(7,400))
#> Initialising...Done.
#> Optimising...
#> h = 7 ; lambda = 400
#> h = 47 ; lambda = 400
#> h = 7 ; lambda = 440
#> h = 47 ; lambda = 440
#> h = 37 ; lambda = 430
#> h = 77 ; lambda = 390
#> h = 59.5 ; lambda = 402.5
#> h = 67 ; lambda = 420
#> h = 52 ; lambda = 405
#> h = 42 ; lambda = 415
#> h = 57 ; lambda = 410
#> h = 22 ; lambda = 435
#> h = 4.5 ; lambda = 447.5
#> h = 17 ; lambda = 450
#> h = 23.25 ; lambda = 441.25
#> h = 8.25 ; lambda = 446.25
#> h = 29.8125 ; lambda = 434.0625
#> h = 15.4375 ; lambda = 442.1875
#> h = 19.03125 ; lambda = 440.1562
#> h = 17.78125 ; lambda = 433.9062
#> h = 21.88281 ; lambda = 439.4141
#> h = 18.91406 ; lambda = 444.5703
#> h = 17.37109 ; lambda = 449.3555
#> h = 16.0625 ; lambda = 445.3125
#> h = 20.42773 ; lambda = 440.8887
#> h = 20.31055 ; lambda = 445.3027
#> h = 20.9502 ; lambda = 447.876
#> h = 21.82422 ; lambda = 441.6211
#> h = 19.6416 ; lambda = 443.833
#> h = 21.09668 ; lambda = 442.3584
#> h = 20.00537 ; lambda = 443.4644
#> h = 19.88818 ; lambda = 447.8784
#> h = 20.02307 ; lambda = 446.131
#> h = 20.19336 ; lambda = 449.7168
#> h = 20.28735 ; lambda = 452.843
#> h = 20.70972 ; lambda = 450.2673
#> h = 20.50433 ; lambda = 449.6701
#> h = 20.48114 ; lambda = 457.2104
#> h = 20.56644 ; lambda = 463.1642
#> h = 20.34946 ; lambda = 466.3371
#> h = 20.27202 ; lambda = 474.6706
#> h = 20.5511 ; lambda = 484.9918
#> h = 20.68298 ; lambda = 501.0662
#> h = 20.38856 ; lambda = 512.5727
#> h = 20.29962 ; lambda = 537.2769
#> h = 20.71058 ; lambda = 563.6725
#> h = 20.92985 ; lambda = 608.1734
#> h = 20.5465 ; lambda = 644.3841
#> h = 20.47826 ; lambda = 716.043
#> h = 21.10849 ; lambda = 786.9395
#> h = 20.90627 ; lambda = 724.5239
#> h = 20.45467 ; lambda = 832.3935
#> h = 20.21709 ; lambda = 944.5035
#> h = 19.78907 ; lambda = 936.0226
#> h = 19.23047 ; lambda = 1041.772
#> h = 18.96929 ; lambda = 1270.232
#> h = 18.21481 ; lambda = 1547.327
#> h = 17.22819 ; lambda = 1644.596
#> h = 15.73375 ; lambda = 1994.642
#> h = 16.21254 ; lambda = 2150.151
#> h = 18.47598 ; lambda = 1318.867
#> h = 17.48937 ; lambda = 1416.135
#> h = 17.12664 ; lambda = 1350.539
#> h = 15.87885 ; lambda = 1676.268
#> h = 16.52813 ; lambda = 1586.918
#> h = 16.42658 ; lambda = 1292.861
#> h = 16.62698 ; lambda = 1380.795
#> h = 16.02848 ; lambda = 1617.173
#> h = 16.30302 ; lambda = 1550.515
#> h = 16.40187 ; lambda = 1344.392
#> h = 16.43343 ; lambda = 1405.023
#> h = 16.10947 ; lambda = 1574.743
#> h = 16.4976 ; lambda = 1429.282
#> h = 16.62802 ; lambda = 1283.791
#> h = 16.38427 ; lambda = 1483.834
#> h = 16.44844 ; lambda = 1508.092
#> h = 16.43718 ; lambda = 1430.791
#> h = 16.32385 ; lambda = 1485.342
#> h = 16.23697 ; lambda = 1513.373
#> h = 16.27093 ; lambda = 1538.386
#> h = 16.39562 ; lambda = 1457.689
#> h = 16.3352 ; lambda = 1459.198
#> h = 16.31067 ; lambda = 1446.88
#> h = 16.2389 ; lambda = 1474.533
#> h = 16.16053 ; lambda = 1482.955
#> h = 16.14735 ; lambda = 1444.493
#> h = 16.05911 ; lambda = 1424.068
#> h = 15.99722 ; lambda = 1480.568
#> h = 15.84049 ; lambda = 1497.412
#> h = 15.82731 ; lambda = 1458.949
#> h = 15.91062 ; lambda = 1464.951
#> h = 15.60376 ; lambda = 1517.87
#> h = 15.73966 ; lambda = 1499.525
#> h = 15.80978 ; lambda = 1467.065
#> h = 15.79443 ; lambda = 1451.891
#> h = 15.63882 ; lambda = 1501.639
#> h = 15.70677 ; lambda = 1492.467
#> h = 15.7769 ; lambda = 1460.006
#> h = 15.74897 ; lambda = 1489.646
#> h = 15.64595 ; lambda = 1515.048
#> h = 15.76883 ; lambda = 1479.061
#> h = 15.81102 ; lambda = 1476.239
#> h = 15.78496 ; lambda = 1480.296
#> h = 15.80482 ; lambda = 1469.711
#> h = 15.76293 ; lambda = 1484.662
#> h = 15.77906 ; lambda = 1485.898
#> h = 15.77139 ; lambda = 1480.77
#> h = 15.74936 ; lambda = 1485.136
#> h = 15.75826 ; lambda = 1483.926
#> h = 15.76671 ; lambda = 1480.034
#> h = 15.76388 ; lambda = 1483.505
#> h = 15.75075 ; lambda = 1486.661
#> h = 15.76623 ; lambda = 1482.243
#> h = 15.77184 ; lambda = 1481.822
#> h = 15.76845 ; lambda = 1482.348
#> Done.
rbind(hlam1,hlam2,hlam3)
#> h lambda
#> hlam1 12.685536 1106.2945
#> hlam2 9.738602 752.6992
#> hlam3 15.766226 1482.2426
# }