spatstat.core (version 2.3-0)

rPSNCP: 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 numberic 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. Impleneted 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 Varaince-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.

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

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

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

# NOT RUN {
  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/10
  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, "/", e2=10)
  sapply(lambda, integral)
  X <- rPSNCP(lambda, kappa, omega, alpha, win = bei$window, dimyx=dim(z1))
  if(online) {
    plot(X)
    plot(split(X), cex=0.5)
  }
# }