Simulate Neyman-Scott Process

Description

Generate a random point pattern, a realisation of the Neyman-Scott cluster process.

Usage

 rNeymanScott(kappa, rmax, rcluster, win = owin(c(0,1),c(0,1)), ..., lmax=NULL)

Arguments

Details

This algorithm generates a realisation of the general Neyman-Scott process, with the cluster mechanism given by the function rcluster. The clusters must have a finite maximum possible radius rmax.

First, the algorithm generates a Poisson point process of “parent” points with intensity kappa. Here kappa may be a single positive number, a function kappa(x,y), or a pixel image object of class "im" (see im.object). See rpoispp for details.

Second, each parent point is replaced by a random cluster of points. These clusters are combined together to yield a single point pattern which is then returned as the result of rNeymanScott.

The argument rcluster specifies the cluster mechanism. It may be either:

• A function which will be called to generate each random cluster (the offspring points of each parent point). The function should expect to be called in the form rcluster(x0,y0,...) for a parent point at a location (x0,y0). The return value of rcluster should specify the coordinates of the points in the cluster; it may be a list containing elements x,y, or a point pattern (object of class "ppp"). If it is a marked point pattern then the result of rNeymanScott will be a marked point pattern.

• A list(mu, f) where mu specifies the mean number of offspring points in each cluster, and f generates the random displacements (vectors pointing from the parent to the offspring). In this case, the number of offspring in a cluster is assumed to have a Poisson distribution, implying that the Neyman-Scott process is also a Cox process. The first element mu should be either a single nonnegative number (interpreted as the mean of the Poisson distribution of cluster size) or a pixel image or a function(x,y) giving a spatially varying mean cluster size (interpreted in the sense of Waagepetersen, 2007). The second element f should be a function that will be called once in the form f(n) to generate n independent and identically distributed displacement vectors (i.e. as if there were a cluster of size n with a parent at the origin (0,0)). The function should return a point pattern (object of class "ppp") or something acceptable to xy.coords that specifies the coordinates of n points.

If required, the intermediate stages of the simulation (the parents and the individual clusters) can also be extracted from the return value of rNeymanScott through the attributes "parents" and "parentid". The attribute "parents" is the point pattern of parent points. The attribute "parentid" is an integer vector specifying the parent for each of the points in the simulated pattern.

Neyman-Scott models where kappa is a single number and rcluster = list(mu,f) can be fitted to data using the function kppm.

Value

The simulated point pattern (an object of class "ppp").

Additionally, some intermediate results of the simulation are returned as attributes of this point pattern: see Details.

Inhomogeneous Neyman-Scott Processes

There are several different ways of specifying a spatially inhomogeneous Neyman-Scott process:

• The point process of parent points can be 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 according to which the parent points are generated.

• The number of points in a typical cluster can be spatially varying. If the argument rcluster is a list of two elements mu, f and the first entry mu is a function(x,y) or a pixel image (object of class "im"), then mu is interpreted as the reference intensity 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(x, y) * g(x-x0, y-y0) where g is the probability density of the offspring displacements generated by the function f.

  Equivalently, clusters are first generated with a constant expected number of points per cluster: the constant is mumax, the maximum of mu. Then the offspring are randomly thinned (see rthin) with spatially-varying retention probabilities given by mu/mumax.

• The entire mechanism for generating a cluster can be dependent on the location of the parent point. If the argument rcluster is a function, then the cluster associated with a parent point at location (x0,y0) will be generated by calling rcluster(x0, y0, ...). The behaviour of this function could depend on the location (x0,y0) in any fashion.

Note that if kappa is an image, the spatial domain covered by this image must be large enough to include the expanded window in which the parent points are to be generated. This expanded window consists of as.rectangle(win) extended by the amount rmax in each direction. This requirement means that win must be small enough so that the expansion of as.rectangle(win) is contained in the spatial domain of kappa. As a result, one may wind up having to simulate the process in a window smaller than what is really desired.

In the first two cases, the intensity of the Neyman-Scott process is equal to kappa * mu if at least one of kappa or mu is a single number, and is otherwise equal to an integral involving kappa, mu and f.

Author(s)

Adrian Baddeley Adrian.Baddeley@uwa.edu.au http://www.maths.uwa.edu.au/~adrian/ and Rolf Turner r.turner@auckland.ac.nz

References

Neyman, J. and Scott, E.L. (1958) A statistical approach to problems of cosmology. Journal of the Royal Statistical Society, Series B 20, 1–43.

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

See Also

rpoispp, rThomas, rGaussPoisson, rMatClust, rCauchy, rVarGamma

Examples

  # each cluster consist of 10 points in a disc of radius 0.2
  nclust <- function(x0, y0, radius, n) {
              return(runifdisc(n, radius, centre=c(x0, y0)))
            }
  plot(rNeymanScott(10, 0.2, nclust, radius=0.2, n=5))

  # multitype Neyman-Scott process (each cluster is a multitype process)
  nclust2 <- function(x0, y0, radius, n, types=c("a", "b")) {
     X <- runifdisc(n, radius, centre=c(x0, y0))
     M <- sample(types, n, replace=TRUE)
     marks(X) <- M
     return(X)
  }
  plot(rNeymanScott(15,0.1,nclust2, radius=0.1, n=5))