rPSNCP {spatstat.random}

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 Adrian Baddeley [e-mail], Rolf Turner [e-mail] and Ege Rubak [e-mail].

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:2]/10
     kappa  <- kappa[1:2]
     omega  <- omega[1:2]
     alpha  <- alpha[1:2, 1:2]
  }
  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:2], "/", e2=10)
     kappa  <- kappa[1:2]
     omega  <- omega[1:2]
     alpha  <- alpha[1:2, 1:2]
  } 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)
  }