rlgcp: Generate log-Gaussian Cox point patterns

View source: R/rlgcp.R

rlgcpR Documentation

Generate log-Gaussian Cox point patterns

Description

Generate one (or several) realisation(s) of the log-Gaussian cox process in a region S\times T.

Usage

 rlgcp(s.region, t.region, replace=TRUE, npoints=NULL, nsim=1, 
 nx=100, ny=100, nt=100,separable=TRUE,model="exponential",
 param=c(1,1,1,1,1,1), scale=c(1,1),var.grf=1,mean.grf=0,
 lmax=NULL,discrete.time=FALSE,exact=FALSE,anisotropy=FALSE,ani.pars=NULL)

Arguments

s.region

Two-column matrix specifying polygonal region containing all data locations. If s.region is missing, the unit square is considered.

t.region

Vector containing the minimum and maximum values of the time interval. If t.region is missing, the interval [0,1] is considered.

npoints

Number of points to simulate. If NULL, the number of points is from a Poisson distribution with mean the double integral of \lambda(s,t) over s.region and t.region.

nsim

number of simulations to generate. Default is 1.

separable

Logical. If TRUE, the covariance function of the Gaussian random field is separable.

model

Vector of length 1 or 2 specifying the model(s) of covariance of the Gaussian random field. If separable=TRUE and model is of length 2, then the elements of model define the spatial and temporal covariances respectively. If separable=TRUE and model is of length 1, then the spatial and temporal covariances belongs to the same class of covariances, among "matern", "exponential", "stable", "cauchy" and "wave" (see Details). If separable=FALSE, model must be of length 1 and is either "gneiting" or "cesare" (see Details).

param

(\alpha_1,\alpha_2,\alpha_3,\alpha_4,\alpha_5,\alpha_6). Vector of parameters of the covariance function (see Details).

scale

Vector of length 2 defining the spatial and temporal scale.

var.grf

Variance of the Gaussian random field.

mean.grf

Mean of the Gaussian random field.

replace

Logical allowing times repeat.

nx, ny, nt

Define the size of the 3-D grid on which the intensity is evaluated.

lmax

Upper bound for the value of \lambda(x,y,t).

discrete.time

If TRUE, times belong to N, otherwise belong to R^+.

exact

logical allowing exact simulation of Gaussian random fields (see manual for details).

anisotropy

If TRUE, simulate an anisotropic point pattern. Currently only implemented for separable covariance functions.

ani.pars

Vector of length 2, the anisotropy angle and the anisotropy ratio, respectively (see details).

Details

We implemented stationary, isotropic spatio-temporal covariance functions.

Separable covariance functions

c(h,t) = c_s(\| h \|) \, c_t(|t|) , h \in S, t \in T

The purely spatial and purely temporal covariance functions can be:

  • Exponential: c(r) = \exp(-r),

  • Stable: c(r) = \exp(-r^\alpha), \alpha \in [0,2],

  • Cauchy: c(r) = (1+r^2)^{-\alpha}, \alpha >0,

  • Wave: c(r) = \frac{\sin(r)}{r} if r>0, c(0)=1,

  • Matern: c(r) = \frac{(\alpha r)^\nu}{2^{\nu-1}\Gamma(\nu)} {\cal K}_{\nu}(\alpha r), \nu > 0 and \alpha > 0.

    {\cal K}_{\nu} is the modified Bessel function of second kind:

    {\cal K}_{\nu}(x) = \frac{\pi}{2} \frac{I_{-\nu}(x) - I_{\nu}(x)}{\sin(\pi \nu)},

    with I_{\nu}(x) = \left( \frac{x}{2} \right)^{\nu} \sum_{k=0}^\infty \frac{1}{k! \Gamma(\nu+k+1)} \left( \frac{x}{2} \right)^{2k}.

    The parameters \alpha_1 and \alpha_2 correspond to the parameters of the spatial and temporal covariance respectively. For the Matern model, the parameters \alpha_1, \alpha_3 and \alpha_2, \alpha_4 correspond to the parameters \nu, \alpha of the spatial and temporal covariance.

Non-separable covariance functions

The spatio-temporal covariance function can be:

  • Gneiting: c(h,t) = \psi (t^2/\beta_2)^{-\alpha_6} \phi \left(\frac{h^2/ \beta_1}{\psi (t^2/\beta_2)} \right), \beta_1, \beta_2 >0 are spatial and temporal scales respectively,

    • If \alpha_2=1, \phi(r) is the Stable model.

    • if \alpha_2=2, \phi(r) is the Cauchy model.

    • If \alpha_2=3, \phi(r) is the Matern model.

    • If \alpha_5=1, \psi(r) = (r^{\alpha_3}+ 1)^{\alpha_4},

    • If \alpha_5=2, \psi(r) = (\alpha_4^{-1} r^{\alpha_3} +1)/(r^{\alpha_3}+1),

    • If \alpha_5=3, \psi(r) = \log(r^{\alpha_3} + \alpha_4)/\log {\alpha_4},

    The parameter \alpha_1 is the respective parameter for the model of \phi(\cdot), \alpha_3 \in (0,1], \alpha_4 \in (0,1] and \alpha_6 \geq 2.

  • cesare: c(h,t) = \left( 1 + (h/\beta_1)^{\alpha_1} + (t/\beta_2)^{\alpha_2} \right)^{-\alpha_3}, \beta_1, \beta_2 >0, \alpha_1, \alpha_2 \in [1,2] and \alpha_3 \geq 3/2.

We also implemented anisotropic Log-Gaussian Cox processes. We considered geometric spatial anisotropy (see Moller and Toftaker, 2014). In this case the covariance function is elliptical and anisotropy is characterized by two parameters: the anisotropy angle 0 \leq \theta < \pi and the anisotropy ratio 0 < \delta \leq 1 of the minor axis 2 \omega \delta and the major axis 2 \omega.

C(h,t)=C_0\left( \sqrt{h \Sigma^{-1} h'},t \right), \ h \in R^2.

Value

A list containing:

xyt

Matrix (or list of matrices if nsim>1) containing the points (x,y,t) of the simulated point pattern. xyt (or any element of the list if nsim>1) is an object of the class stpp.

s.region, t.region

parameters passed in argument.

Lambda

nx \times ny \times nt array (or list of array if nsim>1) of the intensity.

Author(s)

Edith Gabriel <edith.gabriel@inrae.fr>, Peter J Diggle.

References

Chan, G. and Wood A. (1997). An algorithm for simulating stationary Gaussian random fields. Applied Statistics, Algorithm Section, 46, 171–181.

Chan, G. and Wood A. (1999). Simulation of stationary Gaussian vector fields. Statistics and Computing, 9, 265–268.

Gneiting T. (2002). Nonseparable, stationary covariance functions for space-time data. Journal of the American Statistical Association, 97, 590–600.

Moller J. and Toftaker H. (2014). Geometric anisotropic spatial point pattern analysis and Cox processes. Scandinavian Journal of Statistics, 41, 414–435.

See Also

plot.stpp, animation and stan for plotting space-time point patterns.

Examples


# non separable covariance function: 
lgcp1 <- rlgcp(npoints=200, nx=50, ny=50, nt=50, separable=FALSE, 
model="gneiting", param=c(1,1,1,1,1,2), var.grf=1, mean.grf=0)
N <- lgcp1$Lambda[,,1];for(j in 2:(dim(lgcp1$Lambda)[3])){N <-
N+lgcp1$Lambda[,,j]}
image(N,col=grey((1000:1)/1000));box()
animation(lgcp1$xyt, cex=0.8, runtime=10, add=TRUE, prevalent="orange")

# separable covariance function: 
lgcp2 <- rlgcp(npoints=200, nx=50, ny=50, nt=50, separable=TRUE, 
model="exponential", param=c(1,1,1,1,1,2), var.grf=2, mean.grf=-0.5*2)
N <- lgcp2$Lambda[,,1];for(j in 2:(dim(lgcp2$Lambda)[3])){N <-
N+lgcp2$Lambda[,,j]}
image(N,col=grey((1000:1)/1000));box()
animation(lgcp2$xyt, cex=0.8, pch=20, runtime=10, add=TRUE,
prevalent="orange")

# anisotropic

sigma2=0.5
simlgcp <- rlgcp(npoints=500,nx=250, ny=200, nt=50,separable=TRUE,
s.region=matrix(c(0,2,2,0,0,0,0.5,0.5),byrow=FALSE,ncol=2), model="exponential", 
param=c(1,1,1,1,1,2), var.grf=sigma2, mean.grf=-0.5*sigma2,anisotropy = TRUE,
ani.pars = c(pi/4,0.1))

N <- simlgcp$Lambda[,,1];for(j in 2:dim(simlgcp$Lambda)[3]){N <- N+simlgcp$Lambda[,,j]}
image(x=simlgcp$grid[[1]]$x,y=simlgcp$grid[[1]]$y,z=N,col=grey((1000:1)/1000));box()
points(simlgcp$xyt[,1:2],pch=19,cex=0.25,col=2)



stpp documentation built on June 28, 2024, 9:11 a.m.