rThomas: Simulate Thomas Process

In spatstat.random: Random Generation Functionality for the 'spatstat' Family

Simulate Thomas Process

Description

Generate a random point pattern, a simulated realisation of the Thomas cluster process.

Usage

  rThomas(kappa, scale, mu, win = square(1),
             nsim=1, drop=TRUE,
             ...,
             algorithm=c("BKBC", "naive"),
             nonempty=TRUE, 
             poisthresh=1e-6,
             expand = 4*scale,
             saveparents=FALSE, saveLambda=FALSE,
             kappamax=NULL, mumax=NULL, sigma) 

Arguments

kappa	Intensity of the Poisson process of cluster centres. A single positive number, a function, or a pixel image.
scale	Cluster size. Standard deviation of random displacement (along each coordinate axis) of a point from its cluster centre. A single positive number.
mu	Mean number of points per cluster (a single positive number) or reference intensity for the cluster points (a function or a pixel image).
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.
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 in the window. If FALSE, parents are generated even if they have no offspring in the window. The default is recommended unless you need to simulate all the parent points for some other purpose.
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.
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.
sigma	Deprecated. Equivalent to scale.

Details

This algorithm generates a realisation of the (‘modified’) Thomas process, a special case of the Neyman-Scott process, inside the window win.

In the simplest case, where kappa and mu are single numbers, the cluster process is formed by first generating a uniform Poisson point process of “parent” points with intensity kappa. Then each parent point is replaced by a random cluster of “offspring” points, the number of points per cluster being Poisson (mu) distributed, and their positions being isotropic Gaussian displacements from the cluster parent location. The resulting point pattern is a realisation of the classical “stationary Thomas process” generated inside the window win. This point process has intensity kappa * mu.

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 algorithm can also generate spatially inhomogeneous versions of the Thomas 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 (2007). For a given parent point, the offspring constitute a Poisson process with intensity function equal to mu * f, where f is the Gaussian probability density centred at the parent point. Equivalently we first generate, for each parent point, a Poisson (mumax) random number of offspring (where M is the maximum value of mu) with independent Gaussian displacements from the parent location, and then randomly thin the offspring points, with retention probability mu/M.

• Both the parent points and the offspring points can be spatially inhomogeneous, as described above.

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 simulation algorithm first expands the original window win by a distance 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.

The intensity of the Thomas process is kappa * mu if either kappa or mu is a single number. In the general case the intensity is an integral involving kappa, mu and f.

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 Thomas 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 thomas.estK or thomas.estpcf.

Currently spatstat does not support fitting the Thomas cluster process model with inhomogeneous parents.

A Thomas cluster process model fitted by kppm can be simulated automatically using simulate.kppm (which invokes rThomas to perform the simulation).

Author(s)

\adrian

, \rolf and \yamei.

References

\baddchangclustersim

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

Diggle, P. J., Besag, J. and Gleaves, J. T. (1976) Statistical analysis of spatial point patterns by means of distance methods. Biometrics 32 659–667.

Thomas, M. (1949) A generalisation of Poisson's binomial limit for use in ecology. Biometrika 36, 18–25.

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

See Also

rpoispp, rMatClust, rCauchy, rVarGamma, rNeymanScott, rGaussPoisson.

For fitting the model, see kppm, clusterfit.

Examples

  #homogeneous
  X <- rThomas(10, 0.2, 5)
  #inhomogeneous
  Z <- as.im(function(x,y){ 5 * exp(2 * x - 1) }, owin())
  Y <- rThomas(10, 0.2, Z)

spatstat.random documentation built on Oct. 22, 2023, 1:17 a.m.