1

It seems that there are methods for using a set number of points when simulating some ppm objects (for example, through manipulating the rmhcontrol options), but I don't see a similar option for LGCP.

After looking through the documentation and source code, it looks like the number of points for an rLGCP is controlled by an integration of lambda and the area of the window being simulated (per documentation from rpoispp and the call to mu=integrate(lambda)). This is then used as the "mu" parameter in a call to rpois, which determines the number of points used in the simulation (nn=rpois(mu,nsims)). So it seems that the only way currently available to control the number of points directly would be to edit the various function calls and subsequent "nn" variable to take on a user defined constant value as opposed to being a random draw based on lambda. Is that correct?

Alternatively, I was considering a work around using random thinning, such that I'd simulate first, then thin down to the desired number of points.

Eric Robinson
  • 47
  • 4

2 Answers2

0

This is a problem of "conditional simulation" - you want to generate simulated realisations of a point process model, subject to the condition that the number of points is exactly equal to n.

This task is easy for Poisson point process models, and relatively straightforward for Gibbs point process models, which are the two types of models represented by ppm objects.

However, conditional simulation is much harder for a Cox process (including log-Gaussian Cox process LGCP). There are solutions, but they are either mathematically complicated or computationally intensive. So we have not implemented them in spatstat.

(Of course you could just hack the code to make it produce patterns with a specified number of points, but this would not be a valid conditional simulation of the LGCP model).

But if you really desperately need to do conditional simulation, try the following brute-force method. Here model is the fitted model of class kppm, and nfix is the desired number of points in the pattern.

bruteBayes <- function(model, nfix, ...) {
   stopifnot(is.kppm(model))
   ## stability factor: probability of exactly nfix points for Poisson
   logpn0 <- dpois(nfix, lambda=integrate(intensity(model)), log=TRUE)
   ## generate a large number of realisations of the random driving intensity
   Xlist <- simulate(model, nsim=1000, saveLambda=TRUE, ...)
   lamlist <- lapply(Xlist, getElement, name="Lambda")
   ## compute probability of exactly nfix points from each intensity
   mu <- sapply(lamlist, integrate)
   logpn <- dpois(nfix, lambda=mu, log=TRUE)
   ## stabilise
   pna <- exp(logpn - logpn0)
   ## choose one intensity at random, probability proportional to pn
   i <- sample(length(pna), 1, prob=pn)
   lami <- lamlist[[i]]
   ## generate 
   Y <- rpoint(nfix, lami)
   return(Y)
}

This may not work in large problems due to numerical underflow.

Adrian Baddeley
  • 2,534
  • 1
  • 5
  • 8
  • Can you comment (in layman' terms) a bit further about how "hacking" it to simulate a set number of points wouldn't be a valid conditional simulation? I guess I viewed the simulation as a sampling of the underlying intensity, which would be the same regardless of how many "observers" you have sampling the pattern. I could see a case where you try to simulate too many points such that you can't fit them into the observational window and satisfy the underlying density predicted by the intensity. Again... my exploration of this topic is nascent, so appreciate your comments. – Eric Robinson Apr 05 '21 at 13:11
0

update:

conditional simulation for Cox processes (holding the number of points fixed) is now implemented in simulate.kppm in the latest release of the package spatstat.core. The argument n.cond specifies the fixed number of points.

Adrian Baddeley
  • 2,534
  • 1
  • 5
  • 8