2

I'm reading the documentation on the Kest function (page 731) from the spatstat package.

For the argument $r$, that is the "Vector of values for the argument $r$ at which $K(r)$ should be evaluated", it states that:

Users are advised not to specify this argument; there is a sensible default.

I'm interested in knowing how this argument is specified. I went to look at the source for this function but I'm not an expert on R (I write in Python) and I can't decipher how it is calculated.

I know that the piece of code that determines this is here:

  rmaxdefault <- rmax %orifnull% rmax.rule("K", W, lambda)
  if(is.infinite(rmaxdefault)) rmaxdefault <- diameter(W)
  breaks <- handle.r.b.args(r, breaks, W, rmaxdefault=rmaxdefault)
  r <- breaks$r
  rmax <- breaks$max

  (...)

  # recommended range of r values
  alim <- c(0, min(rmax, rmaxdefault))

Assuming rmax=NULL (the default) and a rectangular window, could you explain how the final vector alim is obtained?

Gabriel
  • 3,072
  • 1
  • 22
  • 49
  • Just realized that this has been asked before over at SO: https://stackoverflow.com/q/14330415/1391441 – Gabriel Jul 07 '20 at 14:42
  • 2
    It's not necessarily a problem that this is a duplicate of a question *on another site*. If you want to delete this, you can. Alternatively, you can answer it yourself & reference the SO thread. – gung - Reinstate Monica Jul 07 '20 at 14:46

1 Answers1

3

The answer is in spatstat's FAQ:

  1. How are the r values determined in Kest(X) ?

The default r values for Kest are computed as follows:

  1. The maximum r value is computed by the function 'rmax.rule', rmax \<- rmax.rule("K", W, lambda) where W is the window containing the data, and lambda is the average density of points per unit area. Currently this rule takes the minimum of
    • Ripley’s rule of thumb: rmax = one quarter of the smallest side of the enclosing rectangle
    • large sample rule: rmax = sqrt(1000/(pi \* lambda))
  2. r values are equally spaced from 0 to rmax with step value eps. If eps is not specified, then eps = rmax/512 so that there are 513 values or 512 intervals.

You can always override the r values if you need to.

This means that, for a window of (length,width)=[0,1]x[0,1] (ie: A=1) and assuming that the number of points is N , the array of radii values is obtained as follows (in pseudo-code):

rho = N / A  # density
thumb = A / 4 # rule of thumb
large = sqrt(1000 / (pi * rho)) # rule of large numbers
rmax = min(thumb, large)

Hence for N<~5000 the array of radii will be an equal-spaced segmentation of the range [0, 0.25] in 512 intervals.

This was also asked over at StackOverflow.

Gabriel
  • 3,072
  • 1
  • 22
  • 49