rMatClust: Simulate Matern Cluster Process

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

Simulate Matern Cluster Process

Description

Generate a random point pattern, a simulated realisation of the \Matern Cluster Process.

Usage

 rMatClust(kappa, scale, mu, win = square(1),
           nsim=1, drop=TRUE, ...,
           n.cond=NULL, w.cond=NULL,
           algorithm=c("BKBC", "naive"),
           nonempty=TRUE, 
           poisthresh=1e-6, 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	Radius of the clusters. 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.
n.cond	Optional. Integer specifying a fixed number of points. See the section on Conditional Simulation.
w.cond	Optional. Conditioning region. A window (object of class "owin") specifying the region which must contain exactly n.cond points. See the section on Conditional Simulation.
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.
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 \Matern's cluster 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 placed and uniformly inside a disc of radius scale centred on the parent point. The resulting point pattern is a realisation of the classical “stationary \Matern cluster process”. This point process has intensity kappa * mu.

The algorithm can also generate spatially inhomogeneous versions of the \Matern 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 (2007). For a given parent point, the offspring constitute a Poisson process with intensity function equal to mu/(pi * scale^2) inside the disc of radius scale centred on the parent point, and zero intensity outside this disc. Equivalently we first generate, for each parent point, a Poisson (M) random number of offspring (where M is the maximum value of mu) placed independently and uniformly in the disc of radius scale centred on 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 inhomogeneous, as described above.

The intensity of the \Matern cluster 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 scale.

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 scale. 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 scale 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 \Matern cluster process 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 matclust.estK or matclust.estpcf.

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

A fitted \Matern cluster process model can be simulated automatically using simulate.kppm (which invokes rMatClust to perform the simulation).

Conditional Simulation

If n.cond is specified, it should be a single integer. Simulation will be conditional on the event that the pattern contains exactly n.cond points (or contains exactly n.cond points inside the region w.cond if it is given).

Conditional simulation uses the rejection algorithm described in Section 6.2 of Moller, Syversveen and Waagepetersen (1998). There is a maximum number of proposals which will be attempted. Consequently the return value may contain fewer than nsim point patterns.

Author(s)

\adrian

, \yamei and \rolf.

References

\baddchangclustersim

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

\Matern

, B. (1960) Spatial Variation. Meddelanden fraan Statens Skogsforskningsinstitut, volume 59, number 5. Statens Skogsforskningsinstitut, Sweden.

\Matern

, B. (1986) Spatial Variation. Lecture Notes in Statistics 36, Springer-Verlag, New York.

\Moller

, J., Syversveen, A. and Waagepetersen, R. (1998) Log Gaussian Cox Processes. Scandinavian Journal of Statistics 25, 451–482.

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

See Also

rpoispp, rThomas, rCauchy, rVarGamma, rNeymanScott, rGaussPoisson.

For fitting the model, see kppm, clusterfit.

Examples

 # homogeneous
 X <- rMatClust(10, 0.05, 4)
 # inhomogeneous
 ff <- function(x,y){ 4 * exp(2 * abs(x) - 1) }
 Z <- as.im(ff, owin())
 Y <- rMatClust(10, 0.05, Z)
 YY <- rMatClust(ff, 0.05, 3)

spatstat.random documentation built on Sept. 30, 2024, 9:46 a.m.