rPSNCP: Simulate Product Shot-noise Cox Process

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

Simulate Product Shot-noise Cox Process

Description

Generate a random multitype point pattern, a realisation of the product shot-noise Cox process.

Usage

 rPSNCP(lambda=rep(100, 4), kappa=rep(25, 4), omega=rep(0.03, 4), 
        alpha=matrix(runif(16, -1, 3), nrow=4, ncol=4), 
        kernels=NULL, nu.ker=NULL, win=owin(), nsim=1, drop=TRUE,
        ...,
        cnames=NULL, epsth=0.001)

Arguments

lambda	List of intensities of component processes. Either a numeric vector determining the constant (homogeneous) intensities or a list of pixel images (objects of class "im") determining the (inhomogeneous) intensity functions of component processes. The length of lambda determines the number of component processes.
kappa	Numeric vector of intensities of the Poisson process of cluster centres for component processes. Must have the same size as lambda.
omega	Numeric vector of bandwidths of cluster dispersal kernels for component processes. Must have the same size as lambda and kappa.
alpha	Matrix of interaction parameters. Square numeric matrix with the same number of rows and columns as the length of lambda, kappa and omega. All entries of alpha must be greater than -1.
kernels	Vector of character string determining the cluster dispersal kernels of component processes. Implemented kernels are Gaussian kernel ("Thomas") with bandwidth omega, Variance-Gamma (Bessel) kernel ("VarGamma") with bandwidth omega and shape parameter nu.ker and Cauchy kernel ("Cauchy") with bandwidth omega. Must have the same length as lambda, kappa and omega.
nu.ker	Numeric vector of bandwidths of shape parameters for Variance-Gamma kernels.
win	Window in which to simulate the pattern. An object of class "owin".
nsim	Number of simulated realisations to be generated.
cnames	Optional vector of character strings giving the names of the component processes.
...	Optional arguments passed to as.mask to determine the pixel array geometry. See as.mask.
epsth	Numerical threshold to determine the maximum interaction range for cluster kernels.
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.

Details

This function generates a realisation of a product shot-noise Cox process (PSNCP). This is a multitype (multivariate) Cox point process in which each element of the multivariate random intensity \Lambda(u) of the process is obtained by

\Lambda_i(u) = \lambda_i(u) S_i(u) \prod_{j \neq i} E_{ji}(u)

where \lambda_i(u) is the intensity of component i of the process,

S_i(u) = \frac{1}{\kappa_{i}} \sum_{v \in \Phi_i} k_{i}(u - v)

is the shot-noise random field for component i and

E_{ji}(u) = \exp(-\kappa_{j} \alpha_{ji} / k_{j}(0)) \prod_{v \in \Phi_{j}} {1 + \alpha_{ji} \frac{k_j(u-v)}{k_j(0)}}

is a product field controlling impulses from the parent Poisson process \Phi_j with constant intensity \kappa_j of component process j on \Lambda_i(u). Here k_i(u) is an isotropic kernel (probability density) function on R^2 with bandwidth \omega_i and shape parameter \nu_i, and \alpha_{ji}>-1 is the interaction parameter.

Value

A point pattern (an object of class "ppp") if nsim=1, or a list of point patterns if nsim > 1. Each point pattern is multitype (it carries a vector of marks which is a factor).

Author(s)

Abdollah Jalilian. Modified by \spatstatAuthors.

References

Jalilian, A., Guan, Y., Mateu, J. and Waagepetersen, R. (2015) Multivariate product-shot-noise Cox point process models. Biometrics 71(4), 1022–1033.

See Also

rmpoispp, rThomas, rVarGamma, rCauchy, rNeymanScott

Examples

  online <- interactive()
  # Example 1: homogeneous components
  lambda <- c(250, 300, 180, 400)
  kappa <- c(30, 25, 20, 25)
  omega <- c(0.02, 0.025, 0.03, 0.02)
  alpha <- matrix(runif(16, -1, 1), nrow=4, ncol=4)
  if(!online) {
     lambda <- lambda[1:3]/10
     kappa  <- kappa[1:3]
     omega  <- omega[1:3]
     alpha  <- alpha[1:3, 1:3]
  }
  X <- rPSNCP(lambda, kappa, omega, alpha)
  if(online) {
    plot(X)
    plot(split(X))
  }

  #Example 2: inhomogeneous components
  z1 <- scaletointerval.im(bei.extra$elev, from=0, to=1)
  z2 <- scaletointerval.im(bei.extra$grad, from=0, to=1)
  if(!online) {
    ## reduce resolution to reduce check time
    z1 <- as.im(z1, dimyx=c(40,80))
    z2 <- as.im(z2, dimyx=c(40,80))
  } 
  lambda <- list(
         exp(-8 + 1.5 * z1 + 0.5 * z2),
         exp(-7.25 + 1 * z1  - 1.5 * z2),
         exp(-6 - 1.5 * z1 + 0.5 * z2),
         exp(-7.5 + 2 * z1 - 3 * z2))
  kappa <- c(35, 30, 20, 25) / (1000 * 500)
  omega <- c(15, 35, 40, 25)
  alpha <- matrix(runif(16, -1, 1), nrow=4, ncol=4)
  if(!online) {
     lambda <- lapply(lambda[1:3], "/", e2=10)
     kappa  <- kappa[1:3]
     omega  <- omega[1:3]
     alpha  <- alpha[1:3, 1:3]
  } else {
     sapply(lambda, integral)
  }
  X <- rPSNCP(lambda, kappa, omega, alpha, win = Window(bei), dimyx=dim(z1))
  if(online) {
    plot(X)
    plot(split(X), cex=0.5)
  }

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