spatstat.explore (version 3.5-2)

densityHeat.ppp: Diffusion Estimate of Point Pattern Intensity

Description

Computes the diffusion estimate of the intensity of a point pattern.

Usage

# S3 method for ppp
densityHeat(x, sigma, ..., weights=NULL,
          connect=8, symmetric=FALSE,
          sigmaX=NULL, k=1, show=FALSE, se=FALSE,
          at=c("pixels", "points"),
          leaveoneout = TRUE,
          extrapolate = FALSE, coarsen = TRUE,
          verbose=TRUE, internal=NULL)

Value

Pixel image (object of class "im") giving the estimated intensity of the point process.

If se=TRUE, the result has an attribute "se"

which is another pixel image giving the estimated standard error.

If at="points" then the result is a numeric vector with one entry for each point of x.

Details

This command computes a diffusion kernel estimate of point process intensity from the observed point pattern x.

The function densityHeat is generic, with methods for point patterns in two dimensions (class "ppp") and point patterns on a linear network (class "lpp"). The function densityHeat.ppp described here is the method for class "ppp". Given a two-dimensional point pattern x, it computes a diffusion kernel estimate of the intensity of the point process which generated x.

Diffusion kernel estimates were developed by Botev et al (2010), Barry and McIntyre (2011) and Baddeley et al (2022).

Barry and McIntyre (2011) proposed an estimator for point process intensity based on a random walk on the pixel grid inside the observation window. Baddeley et al (2022) showed that the Barry-McIntyre method is a special case of the diffusion estimator proposed by Botev et al (2010).

The original Barry-McIntyre algorithm assumes a symmetric random walk (i.e. each possible transition has the same probability \(p\)) and requires a square pixel grid (i.e. equal spacing in the \(x\) and \(y\) directions). Their original algorithm is used if symmetric=TRUE. Use the ... arguments to ensure a square grid: for example, the argument eps specifies a square grid with spacing eps units.

The more general algorithm used here (Baddeley et al, 2022) does not require a square grid of pixels. If the pixel grid is not square, and if symmetric=FALSE (the default), then the random walk is not symmetric, in the sense that the probabilities of different jumps will be different, in order to ensure that the smoothing is isotropic.

This implementation also includes two generalizations to the case of adaptive smoothing (Baddeley et al, 2022).

In the first version of adaptive smoothing, the bandwidth is spatially-varying. The argument sigma should be a pixel image (class "im") or a function(x,y) specifying the bandwidth at each spatial location. The smoothing is performed by solving the heat equation with spatially-varying parameters.

In the second version of adaptive smoothing, each data point in x is smoothed using a separate bandwidth. The argument sigmaX should be a numeric vector specifying the bandwidth for each point of x. The smoothing is performed using the lagged arrival algorithm. The argument sigma can be omitted.

If extrapolate=FALSE (the default), calculations are performed using the Euler scheme for the heat equation. If extrapolate=TRUE, the accuracy of the result will be improved by applying Richardson extrapolation (Baddeley et al, 2022, Section 4). After computing the intensity estimate using the Euler scheme on the desired pixel grid, another estimate is computed using the same method on another pixel grid, and the two estimates are combined by Richardson extrapolation to obtain a more accurate result. The second grid is coarser than the original grid if coarsen=TRUE (the default), and finer than the original grid if coarsen=FALSE. Setting extrapolate=TRUE increases computation time by 35% if coarsen=TRUE and by 400% if coarsen=FALSE.

References

Baddeley, A., Davies, T., Rakshit, S., Nair, G. and McSwiggan, G. (2022) Diffusion smoothing for spatial point patterns. Statistical Science 37 (1) 123--142.

Barry, R.P. and McIntyre, J. (2011) Estimating animal densities and home range in regions with irregular boundaries and holes: a lattice-based alternative to the kernel density estimator. Ecological Modelling 222, 1666--1672.

Botev, Z.I., Grotowski, J.F. and Kroese, D.P. (2010) Kernel density estimation via diffusion. Annals of Statistics 38, 2916--2957.

See Also

density.ppp for the usual kernel estimator, and adaptive.density for the tessellation-based estimator.

Examples

   online <- interactive()
   if(!online) op <- spatstat.options(npixel=32)

   X <- runifpoint(25, letterR)
   Z <- densityHeat(X, 0.2)
   if(online) {
     plot(Z, main="Diffusion estimator")
     plot(X, add=TRUE, pch=16)
     integral(Z) # should equal 25
   }

   Z <- densityHeat(X, 0.2, se=TRUE)
   Zse <- attr(Z, "se")
   if(online) plot(solist(estimate=Z, SE=Zse), main="")

   Zex <- densityHeat(X, 0.2, extrapolate=TRUE)

   ZS <- densityHeat(X, 0.2, symmetric=TRUE, eps=0.125)
   if(online) {
     plot(ZS, main="fixed bandwidth")
     plot(X, add=TRUE, pch=16)
   }

   sig <- function(x,y) { (x-1.5)/10 }
   ZZ <- densityHeat(X, sig)
   if(online) {
     plot(ZZ, main="adaptive (I)")
     plot(X, add=TRUE, pch=16)
   }

   sigX <- sig(X$x, X$y)
   AA <- densityHeat(X, sigmaX=sigX)
   if(online) {
     plot(AA, main="adaptive (II)")
     plot(X, add=TRUE, pch=16)
   }
   if(!online) spatstat.options(op)