R: Simulate Neyman-Scott Point Process with Variance Gamma...

rVarGamma {spatstat.random}

R Documentation

Simulate Neyman-Scott Point Process with Variance Gamma cluster kernel

Description

Generate a random point pattern, a simulated realisation of the Neyman-Scott process with Variance Gamma (Bessel) cluster kernel.

Usage

 rVarGamma(kappa, scale, mu, 
              nu,
              win = square(1),
              nsim=1, drop=TRUE,
              ...,
              algorithm=c("BKBC", "naive"),
              nonempty=TRUE, 
              thresh = 0.001,
              poisthresh=1e-6,
              expand = NULL,
              saveparents=FALSE, saveLambda=FALSE,
              kappamax=NULL, mumax=NULL)

Arguments

`kappa`	Intensity of the Poisson process of cluster centres. A single positive number, a function, or a pixel image.
`scale`	Scale parameter for cluster kernel. Determines the size of clusters. A single positive number, in the same units as the spatial coordinates.
`mu`	Mean number of points per cluster (a single positive number) or reference intensity for the cluster points (a function or a pixel image).
`nu`	Shape parameter for the cluster kernel. A number greater than -1.
`win`	Window in which to simulate the pattern. An object of class `"owin"` or something acceptable to `as.owin`.
`nsim`	Number of simulated realisations to be generated.
`drop`	Logical. If `nsim=1` and `drop=TRUE` (the default), the result will be a point pattern, rather than a list containing a point pattern.
`...`	Passed to `clusterfield` to control the image resolution when `saveLambda=TRUE`, and to `clusterradius` when `expand` is missing or `NULL`.
`algorithm`	String (partially matched) specifying the simulation algorithm. See Details.
`nonempty`	Logical. If `TRUE` (the default), a more efficient algorithm is used, in which parents are generated conditionally on having at least one offspring point. If `FALSE`, parents are generated even if they have no offspring. Both choices are valid; the default is recommended unless you need to simulate all the parent points for some other purpose.
`thresh`	Threshold relative to the cluster kernel value at the origin (parent location) determining when the cluster kernel will be treated as zero for simulation purposes. Will be overridden by argument `expand` if that is given.
`poisthresh`	Numerical threshold below which the model will be treated as a Poisson process. See Details.
`expand`	Window expansion distance. A single number. The distance by which the original window will be expanded in order to generate parent points. Has a sensible default, determined by calling `clusterradius` with the numeric threshold value given in `thresh`.
`saveparents`	Logical value indicating whether to save the locations of the parent points as an attribute.
`saveLambda`	Logical. If `TRUE` then the random intensity corresponding to the simulated parent points will also be calculated and saved, and returns as an attribute of the point pattern.
`kappamax`	Optional. Numerical value which is an upper bound for the values of `kappa`, when `kappa` is a pixel image or a function.
`mumax`	Optional. Numerical value which is an upper bound for the values of `mu`, when `mu` is a pixel image or a function.

Details

This algorithm generates a realisation of the Neyman-Scott process with Variance Gamma (Bessel) cluster kernel, inside the window win.

The process is constructed by first generating a Poisson point process of “parent” points with intensity kappa. Then each parent point is replaced by a random cluster of points, the number of points in each cluster being random with a Poisson (mu) distribution, and the points being placed independently and uniformly according to a Variance Gamma kernel.

Note that, for correct simulation of the model, the parent points are not restricted to lie inside the window win; the parent process is effectively the uniform Poisson process on the infinite plane.

The shape of the kernel is determined by the dimensionless index nu. This is the parameter nu^\prime (nu-prime) \nu^\prime = \alpha/2-1 appearing in equation (12) on page 126 of Jalilian et al (2013).

The scale of the kernel is determined by the argument scale, which is the parameter \eta appearing in equations (12) and (13) of Jalilian et al (2013). It is expressed in units of length (the same as the unit of length for the window win).

The algorithm can also generate spatially inhomogeneous versions of the cluster process:

The parent points can be spatially inhomogeneous. If the argument kappa is a function(x,y) or a pixel image (object of class "im"), then it is taken as specifying the intensity function of an inhomogeneous Poisson process that generates the parent points.
The offspring points can be inhomogeneous. If the argument mu is a function(x,y) or a pixel image (object of class "im"), then it is interpreted as the reference density for offspring points, in the sense of Waagepetersen (2006).

If the pair correlation function of the model is very close to that of a Poisson process, deviating by less than poisthresh, then the model is approximately a Poisson process, and will be simulated as a Poisson process with intensity kappa * mu, using rpoispp. This avoids computations that would otherwise require huge amounts of memory.

Value

A point pattern (an object of class "ppp") if nsim=1, or a list of point patterns if nsim > 1.

Additionally, some intermediate results of the simulation are returned as attributes of this point pattern (see rNeymanScott). Furthermore, the simulated intensity function is returned as an attribute "Lambda", if saveLambda=TRUE.

Simulation Algorithm

Two simulation algorithms are implemented.

The naive algorithm generates the cluster process by directly following the description given above. First the window win is expanded by a distance equal to expand. Then the parent points are generated in the expanded window according to a Poisson process with intensity kappa. Then each parent point is replaced by a finite cluster of offspring points as described above. The naive algorithm is used if algorithm="naive" or if nonempty=FALSE.
The BKBC algorithm, proposed by Baddeley and Chang (2023), is a modification of the algorithm of Brix and Kendall (2002). Parents are generated in the infinite plane, subject to the condition that they have at least one offspring point inside the window win. The BKBC algorithm is used when algorithm="BKBC" (the default) and nonempty=TRUE (the default).

The naive algorithm becomes very slow when scale is large, while the BKBC algorithm is uniformly fast (Baddeley and Chang, 2023).

If saveparents=TRUE, then the simulated point pattern will have an attribute "parents" containing the coordinates of the parent points, and an attribute "parentid" mapping each offspring point to its parent.

If nonempty=TRUE (the default), then parents are generated subject to the condition that they have at least one offspring point in the window win. nonempty=FALSE, then parents without offspring will be included; this option is not available in the BKBC algorithm.

Note that if kappa is a pixel image, its domain must be larger than the window win. This is because an offspring point inside win could have its parent point lying outside win. In order to allow this, the naive simulation algorithm first expands the original window win by a distance equal to expand and generates the Poisson process of parent points on this larger window. If kappa is a pixel image, its domain must contain this larger window.

If the pair correlation function of the model is very close to that of a Poisson process, with maximum deviation less than poisthresh, then the model is approximately a Poisson process. This is detected by the naive algorithm which then simulates a Poisson process with intensity kappa * mu, using rpoispp. This avoids computations that would otherwise require huge amounts of memory.

Fitting cluster models to data

The Variance-Gamma cluster model with homogeneous parents (i.e. where kappa is a single number) where the offspring are either homogeneous or inhomogeneous (mu is a single number, a function or pixel image) can be fitted to point pattern data using kppm, or fitted to the inhomogeneous K function using vargamma.estK or vargamma.estpcf.

Currently spatstat does not support fitting the Variance-Gamma cluster process model with inhomogeneous parents.

A Variance-Gamma cluster process model fitted by kppm can be simulated automatically using simulate.kppm (which invokes rVarGamma to perform the simulation).

Warning

The argument nu is the parameter \nu^\prime (nu-prime) \nu^\prime = \alpha/2-1 appearing in equation (12) on page 126 of Jalilian et al (2013). This is different from the parameter called \nu appearing in equation(14) on page 127 of Jalilian et al (2013), defined by \nu = \alpha-1. This has been a frequent source of confusion.

Author(s)

Original algorithm by Abdollah Jalilian and Rasmus Waagepetersen. Adapted for spatstat by Adrian Baddeley Adrian.Baddeley@curtin.edu.au. Brix-Kendall-Baddeley-Chang algorithm implemented by Adrian Baddeley Adrian.Baddeley@curtin.edu.au and Ya-Mei Chang yamei628@gmail.com.

References

Baddeley, A. and Chang, Y.-M. (2023) Robust algorithms for simulating cluster point processes. Journal of Statistical Computation and Simulation. In Press. DOI 10.1080/00949655.2023.2166045.

Brix, A. and Kendall, W.S. (2002) Simulation of cluster point processes without edge effects. Advances in Applied Probability 34, 267–280.

Jalilian, A., Guan, Y. and Waagepetersen, R. (2013) Decomposition of variance for spatial Cox processes. Scandinavian Journal of Statistics 40, 119-137.

Waagepetersen, R. (2007) An estimating function approach to inference for inhomogeneous Neyman-Scott processes. Biometrics 63, 252–258.

Examples

 # homogeneous
 X <- rVarGamma(kappa=5, scale=2, mu=5, nu=-1/4)
 # inhomogeneous
 ff <- function(x,y){ exp(2 - 3 * abs(x)) }
 fmax <- exp(2)
 Z <- as.im(ff, W= owin())
 Y <- rVarGamma(kappa=5, scale=2, mu=Z, nu=0)
 YY <- rVarGamma(kappa=ff, scale=2, mu=3, nu=0, kappamax=fmax)

[Package spatstat.random version 3.3-1 Index]