rNeymanScott: Simulate Neyman-Scott Process
In spatstat.random: Random Generation Functionality for the 'spatstat' Family

rNeymanScott

R Documentation

Simulate Neyman-Scott Process

Description

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

Usage

 rNeymanScott(kappa, expand, rcluster, win = unit.square(),
              ..., nsim=1, drop=TRUE,
              nonempty=TRUE, saveparents=TRUE,
              kappamax=NULL, mumax=NULL)

Arguments

`kappa`	Intensity of the Poisson process of cluster centres. A single positive number, a function, or a pixel image.
`expand`	Size of the expansion of the simulation window for generating parent points. A single non-negative number.
`rcluster`	A function which generates random clusters, or other data specifying the random cluster mechanism. See Details.
`win`	Window in which to simulate the pattern. An object of class `"owin"` or something acceptable to `as.owin`.
`...`	Arguments passed to `rcluster`.
`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.
`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.
`saveparents`	Logical value indicating whether to save the locations of the parent points as an attribute.
`kappamax`	Optional. Upper bound on the values of `kappa` when `kappa` is a function or pixel image.
`mumax`	Optional. Upper bound on the values of `mu` when `mu=rcluster[[1]]` is a function or pixel image.

Details

This algorithm generates a realisation of the general Neyman-Scott process, with the cluster mechanism given by the function rcluster.

First, the algorithm generates a Poisson point process of “parent” points with intensity kappa in an expanded window as explained below. 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, and the restriction of this pattern to the window win is then returned as the result of rNeymanScott.

The expanded window consists of as.rectangle(win) extended by the amount expand in each direction. The size of the expansion is saved in the attribute "expand" and may be extracted by attr(X, "expand") where X is the generated point pattern.

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

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 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 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

and \rolf

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.

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))

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