[Stpp-commits] r22 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 1 14:52:04 CET 2008
Author: gabriele
Date: 2008-11-01 14:52:03 +0100 (Sat, 01 Nov 2008)
New Revision: 22
Added:
pkg/man/fmd.Rd
pkg/man/northcumbria.Rd
pkg/man/rinfec.Rd
pkg/man/rlgcp.Rd
Modified:
pkg/R/rinfec.r
pkg/R/rlgcp.r
pkg/man/rinter.Rd
pkg/man/rpp.Rd
Log:
Modified: pkg/R/rinfec.r
===================================================================
--- pkg/R/rinfec.r 2008-10-31 10:42:57 UTC (rev 21)
+++ pkg/R/rinfec.r 2008-11-01 13:52:03 UTC (rev 22)
@@ -1,363 +1,364 @@
-rinfec <- function(npoints=NULL,s.region,t.region,nsim=1,alpha,beta,gamma,maxrad,delta,s.distr="exponential",t.distr="uniform",h="step",p="min",recent="all",lambda=NULL,lmax=NULL,nx=100,ny=100,nt=1000,t0,inhibition=FALSE,...)
- {
- #
- # Simulate an infectious point process in a region D x T.
- #
- # Requires Splancs.
- #
- # Arguments:
- #
- # npoints: number of points to simulate.
- #
- # s.region: two columns matrix specifying polygonal region containing
- # all data locations. If poly is missing, the unit square is
- # considered.
- #
- # t.region: vector containing the minimum and maximum values of
- # the time interval.
- #
- # nsim: number of simulations to generate. Default is 1.
- #
- # alpha: numerical value for the latent period.
- #
- # beta: numerical value for the maximum infection rate.
- #
- # gamma: numerical value for the infection period.
- #
- # delta: spatial distance of inhibition.
- #
- # maxrad: single value or 2-vector of spatial and temporal maximum
- # radiation respectively. If single value, the same value is used
- # for space and time.
- #
- # s.distr: spatial distribution. Must be choosen among "uniform", "gaussian",
- # "expoenntial" and "poisson".
- #
- # t.distr: temporal distribution. Must be choosen among "uniform" and "exponential".
- #
- # h: infection rate function which depends on alpha, beta and delta.
- # Must be choosen among "step" and "gaussian".
- #
- # p: Must be choosen among "min", "max" and "prod".
- #
- # recent: If "all" consider all previous events. If is an integer, say N,
- # consider only the N most recent events.
- #
- # lambda: function or matrix define the intensity of a Poisson process if
- # s.distr is Poisson.
- #
- # lmax: upper bound for the value of lambda.
- #
- # nx,ny: define the 2-D grid on which the intensity is evaluated if s.distr
- # is Poisson.
- #
- # nt: used to discretize time to compute the infection rate function.
- #
- # t0: minimum time used to compute the infection rate function.
- # Default in the minimum of t.region.
- #
- # inhibition: Logical. If TRUE, an inhibition process is generated.
- # Otherwise, it is a contagious process.
- #
- #
- # Value:
- # pattern: list containing the points (x,y,t) of the simulated point
- # process.
- #
- ##
- ## E. GABRIEL, 02/04/2007
- ##
- ## last modification: 20/10/2008
- ##
- ##
-
- if (is.null(npoints)) stop("please specify the number of points to generate")
- if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
- if (missing(t.region)) t.region <- c(0,1)
- if (missing(t0)) t0 <- t.region[1]
-
- if (missing(maxrad)) maxrad <- c(0.05,0.05)
- maxrads <- maxrad[1]
- maxradt <- maxrad[2]
-
- if (missing(delta)) delta <- maxrads
-
- if (s.distr=="poisson")
- {
- s.grid <- make.grid(nx,ny,s.region)
- s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
- if (is.function(lambda))
- {
- Lambda <- lambda(as.vector(s.grid$X),as.vector(s.grid$Y),...)
- Lambda <- matrix(Lambda,ncol=ny,byrow=T)
- }
- if (is.matrix(lambda))
- {
- nx <- dim(lambda)[1]
- ny <- dim(lambda)[2]
- Lambda <- lambda
- }
- if (is.null(lambda))
- {
- Lambda <- matrix(npoints/areapl(s.region),ncol=ny,nrow=nx)
- }
- Lambda[s.grid$mask==FALSE] <- NaN
- }
-
- if (h=="step")
- {
- hk <- function(t,t0,alpha,beta,gamma)
- {
- res <- rep(0,length(t))
- mask <- (t <= t0+alpha+gamma) & (t >= t0+alpha)
- res[mask] <- beta
- return(res)
- }
- }
-
- if (h=="gaussian")
- {
- hk <- function(t,t0,alpha,beta,gamma)
- {
- t0 <- alpha+t0
- d <- gamma/8
- top <- beta*sqrt(2*pi)*d
-# res <- top*(1/(sqrt(2*pi)*d))*exp(-(((t-(t0+gamma/2))^2)/(2*d^2)))
- res <- top*dnorm(t,mean=t0+gamma/2,sd=d)
- return(res)
- }
-
- d <- gamma/8
- top <- beta*sqrt(2*pi)*d
- ug <- top*dnorm(0,mean=alpha+gamma/2,sd=d)
- }
-
- if (p=="min")
- {
- pk <- function(mu,recent)
- {
- if (recent=="all")
- res <- min(mu)
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(ITK))
- res <- min(mu[(length(ITK)-recent+1):length(ITK)])
- else
- res <- min(mu)
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
-
- if (p=="max")
- {
- pk <- function(mu,recent)
- {
- if (recent=="all")
- res <- max(mu)
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(ITK))
- res <- max(mu[(length(ITK)-recent+1):length(ITK)])
- else
- res <- max(mu)
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
- if (p=="prod")
- {
- pk <- function(mu,recent)
- {
- if (recent=="all")
- res <- prod(mu)
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(ITK))
- res <- prod(mu[(length(ITK)-recent+1):length(ITK)])
- else
- res <- prod(mu)
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
-
- t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=((t.region[2]-t.region[1])/(nt-1)))
-
- pattern <- list()
- ni <- 1
- while(ni<=nsim)
- {
- xy <- csr(n=1,poly=s.region)
- npts <- 1
- pattern.interm <- cbind(x=xy[1],y=xy[2],t=t0)
- Mu <- rep(0,nt)
- ti <- t0
- ITK <- iplace(t.grid$times,ti,t.grid$tinc)
-
- continue <- FALSE
- while(continue==FALSE)
- {
- while (npts<npoints)
- {
- ht <- hk(t=t.grid$times,t0=ti[npts],alpha=alpha,beta=beta,gamma=gamma)
-# ht <- ht/integrate(hk,lower=t.region[1],upper=t.region[2],subdivisions=nt,t0=t0,alpha=alpha,beta=beta,gamma=gamma)$value
- ht <- ht/(sum(ht)*t.grid$tinc)
-
- mu <- Mu+ht
- if (sum(mu)==0)
- mu2 <- mu
- else
- mu2 <- mu/max(mu)
-
- if (t.distr=="uniform")
- tk <- ti[npts] + runif(1,min=0,max=maxradt)
-
- if (t.distr=="exponential")
- tk <- ti[npts] + rexp(1,rate=1/maxradt)
-# tk <- ti[npts]+rexp(1,rate=1/(sum(mu)*t.grid$tinc))
-
- if (tk>=t.region[2])
- {
- npts <- npoints
- continue <- TRUE
- }
- else
- {
- itk <- iplace(t.grid$times,tk,t.grid$tinc)
-
- if ((mu[itk]==0) || ((h=="gaussian") && (mu[itk]<ug)))
- {
- npts <- npoints
- continue <- TRUE
- }
- else
- {
- uk <- runif(1)
- Mu <- mu
- if (inhibition==FALSE)
- {
- cont <- FALSE
- while(cont==FALSE)
- {
- xpar <- pattern.interm[npts,1]
- ypar <- pattern.interm[npts,2]
-
- out <- TRUE
- while(out==TRUE)
- {
- if (s.distr=="uniform")
- {
- xp <- xpar+runif(1,min=-maxrads,max=maxrads)
- yp <- ypar+runif(1,min=-maxrads,max=maxrads)
- }
- if (s.distr=="gaussian")
- {
- xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
- }
- if (s.distr=="exponential")
- {
- xp <- xpar+sample(c(-1,1),1)*rexp(1,rate=1/maxrads)
- yp <- ypar+sample(c(-1,1),1)*rexp(1,rate=1/maxrads)
- }
- if (s.distr=="poisson")
- {
- if (is.null(lmax))
- lmax <- max(Lambda,na.rm=TRUE)
- retain.eq.F <- FALSE
- while(retain.eq.F==FALSE)
- {
- xy <- csr(poly=s.region,npoints=1)
- xp <- xy[1]
- yp <- xy[2]
- up <- runif(1)
- nix <- iplace(X=s.grid$x,x=xp,xinc=s.grid$xinc)
- niy <- iplace(X=s.grid$y,x=yp,xinc=s.grid$yinc)
- Up <- Lambda[nix,niy]/lmax
- retain <- up <= Up
- if ((retain==FALSE) || is.na(retain))
- retain.eq.F <- FALSE
- else
- retain.eq.F <- TRUE
-
- }
- }
- M <- inout(pts=rbind(c(xp,yp),c(xp,yp)),poly=s.region)
- if ((sqrt((xp - pattern.interm[npts,1])^2 + (yp - pattern.interm[npts,2])^2) < maxrads)) M <- c(M,M)
- if (sum(M)==4)
- out <- FALSE
- }
-
- if (sqrt((xp - pattern.interm[npts,1])^2 + (yp - pattern.interm[npts,2])^2) < delta)
- umax <- 1
- else
- umax <- mu2[itk]
-
- if (uk < umax)
- {
- npts <- npts+1
- ITK <- c(ITK,itk)
- ti <- c(ti,tk)
- pattern.interm <- rbind(pattern.interm,c(xp,yp,tk))
- cont <- TRUE
- }
- }
- }
- else
- {
- xy <- csr(n=1,poly=s.region)
- xp <- xy[1]
- yp <- xy[2]
-
- if (all((sqrt((xp - pattern.interm[,1])^2 + (yp - pattern.interm[,2])^2)) > delta))
- umax <- 1
- else
- umax <- pk(mu=mu2[ITK],recent=recent)
- if (uk < umax)
- {
- npts <- npts+1
- ITK <- c(ITK,itk)
- ti <- c(ti,tk)
- pattern.interm <- rbind(pattern.interm,c(xp,yp,tk))
- }
- }
- }
- }
- }
- continue <- TRUE
- }
- dimnames(pattern.interm) <- list(NULL,c("x","y","t"))
- if (nsim==1)
- {
- pattern <- pattern.interm
- ni <- ni+1
- }
- else
- {
- pattern[[ni]] <- pattern.interm
- ni <- ni+1
- }
- }
- par(cex.lab=1.5,cex.axis=1.5)
- plot(t.grid$times,Mu,type="s",ylim=c(0,max(Mu)),lwd=2,xlab="t",ylab=expression(mu(t)))
- points(ti,Mu[ITK],pch=20,col=2,cex=1.5)
- invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region)))
-}
-
-
+rinfec <- function(npoints, s.region, t.region, nsim=1, alpha, beta, gamma,
+s.distr="exponential", t.distr="uniform", maxrad, delta, h="step", g="min",
+recent=1, lambda=NULL, lmax=NULL, nx=100, ny=100, nt=1000,
+t0, inhibition=FALSE, ...)
+{
+ #
+ # Simulate an infectious point process in a region D x T.
+ #
+ # Requires Splancs.
+ #
+ # Arguments:
+ #
+ # npoints: number of points to simulate.
+ #
+ # s.region: two columns matrix specifying polygonal region containing
+ # all data locations. If poly is missing, the unit square is
+ # considered.
+ #
+ # t.region: vector containing the minimum and maximum values of
+ # the time interval.
+ #
+ # nsim: number of simulations to generate. Default is 1.
+ #
+ # alpha: numerical value for the latent period.
+ #
+ # beta: numerical value for the maximum infection rate.
+ #
+ # gamma: numerical value for the infection period.
+ #
+ # delta: spatial distance of inhibition.
+ #
+ # maxrad: single value or 2-vector of spatial and temporal maximum
+ # radiation respectively. If single value, the same value is used
+ # for space and time.
+ #
+ # s.distr: spatial distribution. Must be choosen among "uniform", "gaussian",
+ # "expoenntial" and "poisson".
+ #
+ # t.distr: temporal distribution. Must be choosen among "uniform" and "exponential".
+ #
+ # h: infection rate function which depends on alpha, beta and delta.
+ # Must be choosen among "step" and "gaussian".
+ #
+ # g: Must be choosen among "min", "max" and "prod".
+ #
+ # recent: If "all" consider all previous events. If is an integer, say N,
+ # consider only the N most recent events.
+ #
+ # lambda: function or matrix define the intensity of a Poisson process if
+ # s.distr is Poisson.
+ #
+ # lmax: upper bound for the value of lambda.
+ #
+ # nx,ny: define the 2-D grid on which the intensity is evaluated if s.distr
+ # is Poisson.
+ #
+ # nt: used to discretize time to compute the infection rate function.
+ #
+ # t0: minimum time used to compute the infection rate function.
+ # Default in the minimum of t.region.
+ #
+ # inhibition: Logical. If TRUE, an inhibition process is generated.
+ # Otherwise, it is a contagious process.
+ #
+ #
+ # Value:
+ # xyt: list containing the points (x,y,t) of the simulated point
+ # process.
+ #
+
+ if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
+ if (missing(t.region)) t.region <- c(0,1)
+ if (missing(t0)) t0 <- t.region[1]
+
+ if (missing(maxrad)) maxrad <- c(0.05,0.05)
+ maxrads <- maxrad[1]
+ maxradt <- maxrad[2]
+
+ if (missing(delta)) delta <- maxrads
+
+ if (s.distr=="poisson")
+ {
+ if (is.matrix(lambda))
+ {
+ nx <- dim(lambda)[1]
+ ny <- dim(lambda)[2]
+ Lambda <- lambda
+ }
+ s.grid <- make.grid(nx,ny,s.region)
+ s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
+ if (is.function(lambda))
+ {
+ Lambda <- lambda(as.vector(s.grid$X),as.vector(s.grid$Y),...)
+ Lambda <- matrix(Lambda,ncol=ny,byrow=T)
+ }
+ if (is.null(lambda))
+ {
+ Lambda <- matrix(npoints/areapl(s.region),ncol=ny,nrow=nx)
+ }
+ Lambda[s.grid$mask==FALSE] <- NaN
+ }
+
+ if (h=="step")
+ {
+ hk <- function(t,t0,alpha,beta,gamma)
+ {
+ res <- rep(0,length(t))
+ mask <- (t <= t0+alpha+gamma) & (t >= t0+alpha)
+ res[mask] <- beta
+ return(res)
+ }
+ }
+
+ if (h=="gaussian")
+ {
+ hk <- function(t,t0,alpha,beta,gamma)
+ {
+ t0 <- alpha+t0
+ d <- gamma/8
+ top <- beta*sqrt(2*pi)*d
+# res <- top*(1/(sqrt(2*pi)*d))*exp(-(((t-(t0+gamma/2))^2)/(2*d^2)))
+ res <- top*dnorm(t,mean=t0+gamma/2,sd=d)
+ return(res)
+ }
+
+ d <- gamma/8
+ top <- beta*sqrt(2*pi)*d
+ ug <- top*dnorm(0,mean=alpha+gamma/2,sd=d)
+ }
+
+ if (g=="min")
+ {
+ pk <- function(mu,recent)
+ {
+ if (recent=="all")
+ res <- min(mu)
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(ITK))
+ res <- min(mu[(length(ITK)-recent+1):length(ITK)])
+ else
+ res <- min(mu)
+ }
+ else
+ stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+ }
+
+ if (g=="max")
+ {
+ pk <- function(mu,recent)
+ {
+ if (recent=="all")
+ res <- max(mu)
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(ITK))
+ res <- max(mu[(length(ITK)-recent+1):length(ITK)])
+ else
+ res <- max(mu)
+ }
+ else
+ stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+ }
+ if (g=="prod")
+ {
+ pk <- function(mu,recent)
+ {
+ if (recent=="all")
+ res <- prod(mu)
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(ITK))
+ res <- prod(mu[(length(ITK)-recent+1):length(ITK)])
+ else
+ res <- prod(mu)
+ }
+ else
+ stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+ }
+
+ t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=((t.region[2]-t.region[1])/(nt-1)))
+
+ pattern <- list()
+ ni <- 1
+ while(ni<=nsim)
+ {
+ xy <- csr(n=1,poly=s.region)
+ npts <- 1
+ pattern.interm <- cbind(x=xy[1],y=xy[2],t=t0)
+ Mu <- rep(0,nt)
+ ti <- t0
+ ITK <- iplace(t.grid$times,ti,t.grid$tinc)
+
+ continue <- FALSE
+ while(continue==FALSE)
+ {
+ while (npts<npoints)
+ {
+ ht <- hk(t=t.grid$times,t0=ti[npts],alpha=alpha,beta=beta,gamma=gamma)
+# ht <- ht/integrate(hk,lower=t.region[1],upper=t.region[2],subdivisions=nt,t0=t0,alpha=alpha,beta=beta,gamma=gamma)$value
+ ht <- ht/(sum(ht)*t.grid$tinc)
+
+ mu <- Mu+ht
+ if (sum(mu)==0)
+ mu2 <- mu
+ else
+ mu2 <- mu/max(mu)
+
+ if (t.distr=="uniform")
+ tk <- ti[npts] + runif(1,min=0,max=maxradt)
+
+ if (t.distr=="exponential")
+ tk <- ti[npts] + rexp(1,rate=1/maxradt)
+# tk <- ti[npts]+rexp(1,rate=1/(sum(mu)*t.grid$tinc))
+
+ if (tk>=t.region[2])
+ {
+ N = npts
+ npts <- npoints
+ continue <- TRUE
+ M = paste("only",N,"points have been generated over",npoints,". Process stopped when getting a time t greater than max(t.region): try a larger t.region or check your parameters",sep=" ")
+ warning(M)
+ }
+ else
+ {
+ itk <- iplace(t.grid$times,tk,t.grid$tinc)
+
+ if ((mu[itk]==0) || ((h=="gaussian") && (mu[itk]<ug)))
+ {
+ N = npts
+ npts <- npoints
+ continue <- TRUE
+ M = paste("only",N,"points have been generated over",npoints,"(process stopped when getting µ(t)=0",sep=" ")
+ warning(M)
+ }
+ else
+ {
+ Mu <- mu
+ if (inhibition==FALSE)
+ {
+ cont <- FALSE
+ while(cont==FALSE)
+ {
+ uk <- runif(1)
+ xpar <- pattern.interm[npts,1]
+ ypar <- pattern.interm[npts,2]
+
+ out <- TRUE
+ while(out==TRUE)
+ {
+ if (s.distr=="uniform")
+ {
+ xp <- xpar+runif(1,min=-maxrads,max=maxrads)
+ yp <- ypar+runif(1,min=-maxrads,max=maxrads)
+ }
+ if (s.distr=="gaussian")
+ {
+ xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
+ yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
+ }
+ if (s.distr=="exponential")
+ {
+ xp <- xpar+sample(c(-1,1),1)*rexp(1,rate=1/maxrads)
+ yp <- ypar+sample(c(-1,1),1)*rexp(1,rate=1/maxrads)
+ }
+ if (s.distr=="poisson")
+ {
+ if (is.null(lmax))
+ lmax <- max(Lambda,na.rm=TRUE)
+ retain.eq.F <- FALSE
+ while(retain.eq.F==FALSE)
+ {
+ xy <- csr(poly=s.region,npoints=1)
+ xp <- xy[1]
+ yp <- xy[2]
+ up <- runif(1)
+ nix <- iplace(X=s.grid$x,x=xp,xinc=s.grid$xinc)
+ niy <- iplace(X=s.grid$y,x=yp,xinc=s.grid$yinc)
+ Up <- Lambda[nix,niy]/lmax
+ retain <- up <= Up
+ if ((retain==FALSE) || is.na(retain))
+ retain.eq.F <- FALSE
+ else
+ retain.eq.F <- TRUE
+
+ }
+ }
+ M <- inout(pts=rbind(c(xp,yp),c(xp,yp)),poly=s.region)
+ if ((sqrt((xp - pattern.interm[npts,1])^2 + (yp - pattern.interm[npts,2])^2) < maxrads)) M <- c(M,M)
+ if (sum(M)==4)
+ out <- FALSE
+ }
+
+ if (all(sqrt((xp - pattern.interm[,1])^2 + (yp - pattern.interm[,2])^2) < delta))
+ umax <- 1
+ else
+ umax <- pk(mu=mu2[ITK],recent=recent)
+# umax <- mu2[itk]
+
+ if (uk < umax)
+ {
+ npts <- npts+1
+ ITK <- c(ITK,itk)
+ ti <- c(ti,tk)
+ pattern.interm <- rbind(pattern.interm,c(xp,yp,tk))
+ cont <- TRUE
+ }
+ }
+ }
+ else
+ {
+ uk <- runif(1)
+ xy <- csr(n=1,poly=s.region)
+ xp <- xy[1]
+ yp <- xy[2]
+
+ if (all((sqrt((xp - pattern.interm[,1])^2 + (yp - pattern.interm[,2])^2)) > delta))
+ umax <- 1
+ else
+ umax <- pk(mu=mu2[ITK],recent=recent)
+ if (uk < umax)
+ {
+ npts <- npts+1
+ ITK <- c(ITK,itk)
+ ti <- c(ti,tk)
+ pattern.interm <- rbind(pattern.interm,c(xp,yp,tk))
+ }
+ }
+ }
+ }
+ }
+ continue <- TRUE
+ }
+ dimnames(pattern.interm) <- list(NULL,c("x","y","t"))
+ if (nsim==1)
+ {
+ pattern <- pattern.interm
+ ni <- ni+1
+ }
+ else
+ {
+ pattern[[ni]] <- pattern.interm
+ ni <- ni+1
+ }
+ }
+ invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region)))
+}
+
+
Modified: pkg/R/rlgcp.r
===================================================================
--- pkg/R/rlgcp.r 2008-10-31 10:42:57 UTC (rev 21)
+++ pkg/R/rlgcp.r 2008-11-01 13:52:03 UTC (rev 22)
@@ -1,246 +1,231 @@
-rlgcp <- function(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,2),scale=c(1,1),var.grf=1,mean.grf=0,lmax=NULL,discrete.time=FALSE,exact=TRUE)
-{
- #
- # Simulate a space-time log-Gaussian point process in a region D x T.
- #
- # Requires Splancs.
- #
- # Arguments:
- #
- # s.region: two columns matrix specifying polygonal region containing
- # all data locations. If poly is missing, the unit square is
- # considered.
- #
- # t.region: vector containing the minimum and maximum values of
- # the time interval.
- #
- # replace: logical allowing times repetition.
- #
- # npoints: number of points to simulate. If NULL (default), the
- # number of points is from a Poisson distribution with
- # mean the double integral of lambda over s.region and
- # t.region.
- #
- # nsim: number of simulations to generate. Default is 1.
- #
- # nx,ny,nt: define the 3-D grid on which the intensity is evaluated.
- #
- #
- # Value:
- # pattern: list containing the points (x,y,t) of the simulated point
- # process.
- # Lambda: an array of the intensity surface at each time.
- #
- # NB: the probability that an event occurs at a location s and a time t
- # is:
- # p(t)=lambda(s,t)/(max_{s_i, i=1,...,ngrid} lambda(s_i,t)).
- #
- ##
- ## E. GABRIEL, 24/01/2006
- ##
- ## last modification: 27/01/2006
- ## 02/02/2007
- ## 13/02/2007
- ##
- ##
-
- if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
- if (missing(t.region)) t.region <- c(0,1)
-
- t.region <- sort(t.region)
- s.area <- areapl(s.region)
- t.area <- t.region[2]-t.region[1]
- tau <- c(start=t.region[1],end=t.region[2],step=(t.region[2]-t.region[1])/(nt-1))
- bpoly <- bbox(s.region)
-
- lambdamax <- lmax
- pattern <- list()
- index.t <- list()
- Lambdafin <- list()
- ni <- 1
-
- s.grid <- make.grid(nx,ny,s.region)
- s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
-
- if (discrete.time==TRUE)
- {
- vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
- if (nt>length(vect))
- {
- nt <- length(vect)
- warning("nt used is less than the one given in argument")
- t.grid <- list(times=vect,tinc=1)
- }
- else
- {
- vect <- round(seq(floor(t.region[1]),ceiling(t.region[2]),length=nt))
- t.grid <- list(times=vect,tinc=round(t.area/(nt-1)))
- }
- }
- else
- t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=(t.area/(nt-1)))
-
- while(ni<=nsim)
- {
-# now <- proc.time()[1]
- S <- gauss3D(nx=nx,ny=ny,nt=nt,xlim=range(s.region[,1]),ylim=range(s.region[,2]),tlim=range(t.region),separable=separable,model=model,param=param,scale=scale,var.grf=var.grf,mean.grf=mean.grf,exact=exact)
-# speed <- proc.time()[1]-now;print(speed)
-
- Lambda <- exp(S)
-
- mut <- rep(0,nt)
- for (it in 1:nt)
- {
- Lambda[,,it][s.grid$mask==FALSE] <- NaN
- mut[it] <- sum(Lambda[,,it],na.rm=TRUE)
- }
-
- if (is.null(npoints))
- {
- en <- sum(Lambda,na.rm=TRUE)*s.grid$xinc*s.grid$yinc
- npoints <- round(rpois(n=1,lambda=en),0)
- }
-
- if (is.null(lambdamax))
- lambdamax <- max(Lambda,na.rm=TRUE)
-# summ <- summary(Lambda)
-# lmax <- summ[[6]] + 0.05 * diff(range(Lambda,na.rm=TRUE))
-
- npts <- round(lambdamax/(s.area*t.area),0)
-## npts <- npoints
- if (npts==0) stop("there is no data to thin")
-
- if ((replace==FALSE) && (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
-
- if (discrete.time==TRUE)
- {
- vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
- times.init <- sample(vect,nt,replace=replace)
- }
- else
- times.init <- runif(nt,min=t.region[1],max=t.region[2])
-
- samp <- sample(1:nt,npts,replace=replace,prob=mut/max(mut,na.rm=TRUE))
- times <- times.init[samp]
-
-# now <- proc.time()[1]
- retain.eq.F <- FALSE
- while(retain.eq.F==FALSE)
- {
- xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
- x <- xy[,1]
- y <- xy[,2]
-
- prob <- NULL
- for(ix in 1:length(x))
- {
- nix <- iplace(X=s.grid$x,x=x[ix],xinc=s.grid$xinc)
- niy <- iplace(X=s.grid$y,x=y[ix],xinc=s.grid$yinc)
- nit <- iplace(X=t.grid$times,x=times[ix],xinc=t.grid$tinc)
- prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
- }
-
- M <- which(is.na(prob))
- if (length(M)!=0)
- {
- x <- x[-M]
- y <- y[-M]
- times <- times[-M]
- prob <- prob[-M]
- npts <- length(x)
- }
-
- u <- runif(npts)
- retain <- u <= prob
- if (sum(retain==FALSE)==length(retain)) retain.eq.F <- FALSE
- else retain.eq.F <- TRUE
- }
-
- x <- x[retain]
- y <- y[retain]
- samp <- samp[retain]
- samp.remain <- (1:nt)[-samp]
- times <- times[retain]
-
- neffec <- length(x)
- if (neffec > npoints)
- {
- retain <- 1:npoints
- x <- x[retain]
- y <- y[retain]
- samp <- samp[retain]
- samp.remain <- (1:nt)[-samp]
- times <- times[retain]
- }
- while(neffec < npoints)
- {
- xy <- as.matrix(csr(poly=s.region,npoints=npoints-neffec))
- if (dim(xy)[2]==1) {wx <- xy[1]; wy <- xy[2]}
- else {wx <- xy[,1]; wy <- xy[,2]}
-
- if (isTRUE(replace))
- wsamp <- sample(1:nt,npoints-neffec,replace=replace,prob=mut/max(mut,na.rm=TRUE))
- else
- wsamp <- sample(samp.remain,npoints-neffec,replace=replace,prob=mut[samp.remain]/max(mut[samp.remain],na.rm=TRUE))
-
- wtimes <- times.init[wsamp]
- prob <- NULL
- for(ix in 1:length(wx))
- {
- nix <- iplace(X=s.grid$x,x=wx[ix],xinc=s.grid$xinc)
- niy <- iplace(X=s.grid$y,x=wy[ix],xinc=s.grid$yinc)
- nit <- iplace(X=t.grid$times,x=wtimes[ix],xinc=t.grid$tinc)
- prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
- }
- M <- which(is.na(prob))
- if (length(M)!=0)
- {
- wx <- wx[-M]
- wy <- wy[-M]
- wtimes <- wtimes[-M]
- prob <- prob[-M]
- }
- if (neffec > 0)
- {
- u <- runif(length(prob))
- retain <- u <= prob
- x <- c(x,wx[retain])
- y <- c(y,wy[retain])
- times <- c(times,wtimes[retain])
- samp <- c(samp,wsamp[retain])
- samp.remain <- (1:nt)[-samp]
- neffec <- length(x)
- }
- }
-# speed <- proc.time()[1]-now;print(speed)
-
- times <- sort(times)
- index.times <- sort(samp)
- pattern.interm <- cbind(x=x,y=y,t=times)
-
- if (nsim==1)
- {
- pattern <- pattern.interm
- index.t <- index.times
- Lambdafin <- Lambda
- }
- else
- {
- pattern[[ni]] <- pattern.interm
- index.t[[ni]] <- index.times
- Lambdafin[[ni]] <- Lambda
- }
- ni <- ni+1
- }
-
- invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region,Lambda=Lambdafin,index.t=index.t)))
-}
-
-
-
-
-
-
-
-
-
+rlgcp <- function(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,2),scale=c(1,1),var.grf=1,mean.grf=0,lmax=NULL,discrete.time=FALSE,exact=FALSE)
+{
+ #
+ # Simulate a space-time log-Gaussian point process in a region D x T.
+ #
+ # Requires Splancs.
+ #
+ # Arguments:
+ #
+ # s.region: two columns matrix specifying polygonal region containing
+ # all data locations. If poly is missing, the unit square is
+ # considered.
+ #
+ # t.region: vector containing the minimum and maximum values of
+ # the time interval.
+ #
+ # replace: logical allowing times repetition.
+ #
+ # npoints: number of points to simulate. If NULL (default), the
+ # number of points is from a Poisson distribution with
+ # mean the double integral of lambda over s.region and
+ # t.region.
+ #
+ # nsim: number of simulations to generate. Default is 1.
+ #
+ # nx,ny,nt: define the 3-D grid on which the intensity is evaluated.
+ #
+ #
+ # Value:
+ # pattern: list containing the points (x,y,t) of the simulated point
+ # process.
+ # Lambda: an array of the intensity surface at each time.
+ #
+ # NB: the probability that an event occurs at a location s and a time t
+ # is:
+ # p(t)=lambda(s,t)/(max_{s_i, i=1,...,ngrid} lambda(s_i,t)).
+ #
+
+ if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
+ if (missing(t.region)) t.region <- c(0,1)
+
+ t.region <- sort(t.region)
+ s.area <- areapl(s.region)
+ t.area <- t.region[2]-t.region[1]
+ tau <- c(start=t.region[1],end=t.region[2],step=(t.region[2]-t.region[1])/(nt-1))
+ bpoly <- bbox(s.region)
+
+ lambdamax <- lmax
+ pattern <- list()
+ index.t <- list()
+ Lambdafin <- list()
+ ni <- 1
+
+ s.grid <- make.grid(nx,ny,s.region)
+ s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
+
+ if (discrete.time==TRUE)
+ {
+ vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
+ if (nt>length(vect))
+ {
+ nt <- length(vect)
+ warning("nt used is less than the one given in argument")
+ t.grid <- list(times=vect,tinc=1)
+ }
+ else
+ {
+ vect <- round(seq(floor(t.region[1]),ceiling(t.region[2]),length=nt))
+ t.grid <- list(times=vect,tinc=round(t.area/(nt-1)))
+ }
+ }
+ else
+ t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=(t.area/(nt-1)))
+
+ while(ni<=nsim)
+ {
+ S <- gauss3D(nx=nx,ny=ny,nt=nt,xlim=range(s.region[,1]),ylim=range(s.region[,2]),tlim=range(t.region),separable=separable,model=model,param=param,scale=scale,var.grf=var.grf,mean.grf=mean.grf,exact=exact)
+
+ Lambda <- exp(S)
+
+ mut <- rep(0,nt)
+ for (it in 1:nt)
+ {
+ Lambda[,,it][s.grid$mask==FALSE] <- NaN
+ mut[it] <- sum(Lambda[,,it],na.rm=TRUE)
+ }
+
+ if (is.null(npoints))
+ {
+ en <- sum(Lambda,na.rm=TRUE)*s.grid$xinc*s.grid$yinc
+ npoints <- round(rpois(n=1,lambda=en),0)
+ }
+
+ if (is.null(lambdamax))
+ lambdamax <- max(Lambda,na.rm=TRUE)
+
+ npts <- round(lambdamax/(s.area*t.area),0)
+ if (npts==0) stop("there is no data to thin")
+
+ if ((replace==FALSE) && (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
+
+ if (discrete.time==TRUE)
+ {
+ vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
+ times.init <- sample(vect,nt,replace=replace)
+ }
+ else
+ times.init <- runif(nt,min=t.region[1],max=t.region[2])
+
+ samp <- sample(1:nt,npts,replace=replace,prob=mut/max(mut,na.rm=TRUE))
+ times <- times.init[samp]
+
+ retain.eq.F <- FALSE
+ while(retain.eq.F==FALSE)
+ {
+ xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
+ x <- xy[,1]
+ y <- xy[,2]
+
+ prob <- NULL
+ for(ix in 1:length(x))
+ {
+ nix <- iplace(X=s.grid$x,x=x[ix],xinc=s.grid$xinc)
+ niy <- iplace(X=s.grid$y,x=y[ix],xinc=s.grid$yinc)
+ nit <- iplace(X=t.grid$times,x=times[ix],xinc=t.grid$tinc)
+ prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
+ }
+
+ M <- which(is.na(prob))
+ if (length(M)!=0)
+ {
+ x <- x[-M]
+ y <- y[-M]
+ times <- times[-M]
+ prob <- prob[-M]
+ npts <- length(x)
+ }
+
+ u <- runif(npts)
+ retain <- u <= prob
+ if (sum(retain==FALSE)==length(retain)) retain.eq.F <- FALSE
+ else retain.eq.F <- TRUE
+ }
+
+ x <- x[retain]
+ y <- y[retain]
+ samp <- samp[retain]
+ samp.remain <- (1:nt)[-samp]
+ times <- times[retain]
+
+ neffec <- length(x)
+ if (neffec > npoints)
+ {
+ retain <- 1:npoints
+ x <- x[retain]
+ y <- y[retain]
+ samp <- samp[retain]
+ samp.remain <- (1:nt)[-samp]
+ times <- times[retain]
+ }
+ while(neffec < npoints)
+ {
+ xy <- as.matrix(csr(poly=s.region,npoints=npoints-neffec))
+ if (dim(xy)[2]==1) {wx <- xy[1]; wy <- xy[2]}
+ else {wx <- xy[,1]; wy <- xy[,2]}
+
+ if (isTRUE(replace))
+ wsamp <- sample(1:nt,npoints-neffec,replace=replace,prob=mut/max(mut,na.rm=TRUE))
+ else
+ wsamp <- sample(samp.remain,npoints-neffec,replace=replace,prob=mut[samp.remain]/max(mut[samp.remain],na.rm=TRUE))
+
+ wtimes <- times.init[wsamp]
+ prob <- NULL
+ for(ix in 1:length(wx))
+ {
+ nix <- iplace(X=s.grid$x,x=wx[ix],xinc=s.grid$xinc)
+ niy <- iplace(X=s.grid$y,x=wy[ix],xinc=s.grid$yinc)
+ nit <- iplace(X=t.grid$times,x=wtimes[ix],xinc=t.grid$tinc)
+ prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
+ }
+ M <- which(is.na(prob))
+ if (length(M)!=0)
+ {
+ wx <- wx[-M]
+ wy <- wy[-M]
+ wtimes <- wtimes[-M]
+ prob <- prob[-M]
+ }
+ if (neffec > 0)
+ {
+ u <- runif(length(prob))
+ retain <- u <= prob
+ x <- c(x,wx[retain])
+ y <- c(y,wy[retain])
+ times <- c(times,wtimes[retain])
+ samp <- c(samp,wsamp[retain])
+ samp.remain <- (1:nt)[-samp]
+ neffec <- length(x)
+ }
+ }
+
+ times <- sort(times)
+ index.times <- sort(samp)
+ pattern.interm <- cbind(x=x,y=y,t=times)
+
+ if (nsim==1)
+ {
+ pattern <- pattern.interm
+ index.t <- index.times
+ Lambdafin <- Lambda
+ }
+ else
+ {
+ pattern[[ni]] <- pattern.interm
+ index.t[[ni]] <- index.times
+ Lambdafin[[ni]] <- Lambda
+ }
+ ni <- ni+1
+ }
+
+ invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region,Lambda=Lambdafin,index.t=index.t)))
+}
+
+
+
+
+
+
+
+
+
Added: pkg/man/fmd.Rd
===================================================================
--- pkg/man/fmd.Rd (rev 0)
+++ pkg/man/fmd.Rd 2008-11-01 13:52:03 UTC (rev 22)
@@ -0,0 +1,21 @@
+\name{fmd}
+\docType{data}
+\alias{fmd}
+\title{2001 food-and-mouth epidemic, north Cumbria (UK)}
+\description{
+This data set gives the spatial locations and reported times of
+food-and-mouth disease in north Cumbria (UK), 2001.
+}
+\usage{data(fmd)}
+\format{A matrix containing (x,y,t) coordinates of the 648 observations.}
+
+\references{
+Diggle, P., Rowlingson, B. and Su, T. (2005). Point process methodology
+for on-line spatio-temporal disease surveillance. Environmetrics, 16, 423--34.
+}
+
+\seealso{
+ \code{\link{northcumbria}} for boundaries of the county of north Cumbria.
+ }
+
+\keyword{datasets}
\ No newline at end of file
Added: pkg/man/northcumbria.Rd
===================================================================
--- pkg/man/northcumbria.Rd (rev 0)
+++ pkg/man/northcumbria.Rd 2008-11-01 13:52:03 UTC (rev 22)
@@ -0,0 +1,15 @@
+\name{northcumbria}
+\docType{data}
+\alias{northcumbria}
+\title{Polygon boundary of north Cumbria}
+\description{
+This data set gives the boundary of the county of north Cumbria (UK).
+}
+\usage{data(northcumbria)}
+\format{A matrix containing (x,y) coordinates of the boundary.}
+\seealso{
+ \code{\link{fmd}} for the space-time pattern of food-and-mouth disease in this
+ county in 2001.
+ }
+
+\keyword{datasets}
\ No newline at end of file
Added: pkg/man/rinfec.Rd
===================================================================
--- pkg/man/rinfec.Rd (rev 0)
+++ pkg/man/rinfec.Rd 2008-11-01 13:52:03 UTC (rev 22)
@@ -0,0 +1,82 @@
+\name{rinfec}
+\alias{rinfec}
+\title{Generate infection point patterns}
+
+\description{
+Generate one (or several) realisation(s) of the infection process
+in a region S x T.
+}
+
+\usage{
+rinfec(npoints, s.region, t.region, nsim=1, alpha, beta, gamma, h="step",
+ s.distr="exponential", t.distr="uniform", maxrad, delta, g="min", recent=1,
+ lambda=NULL, lmax=NULL, nx=100, ny=100, nt=1000, t0, inhibition=FALSE, ...)
+}
+
+\arguments{
+ \item{npoints}{number of points to simulate. }
+ \item{s.region}{two-column matrix specifying polygonal region containing
+ all data locations. If \code{s.region} is missing, the unit square is considered.}
+ \item{t.region}{vector containing the minimum and maximum values of
+ the time interval. If \code{t.region} is missing, the interval [0,1] is considered.}
+ \item{nsim}{number of simulations to generate. Default is 1. }
+ \item{alpha}{numerical value for the latent period.}
+ \item{beta}{numerical value for the maximum infection rate.}
+ \item{gamma}{numerical value for the infection period.}
+ \item{h}{infection rate function which depends on alpha, beta and delta.
+ Must be choosen among "step" and "gaussian".}
+ \item{s.distr}{spatial distribution. Must be choosen among "uniform",
+ "gaussian", "exponential" and "poisson". }
+ \item{t.distr}{temporal distribution. Must be choosen among "uniform" and
+ "exponential".}
+ \item{maxrad}{single value or 2-vector of spatial and temporal maximum
+ radiation respectively. If single value, the same value is used for space and time.}
+ \item{delta}{spatial distance of inhibition/contagion. If missing, the spatial
+ radiation is used}
+ \item{g}{Compute the probability of acceptance of a new point from \code{h}
+ and \code{recent}. Must be choosen among "min", "max" and "prod". }
+ \item{recent}{If ``\code{all}'' consider all previous events. If
+ is an integer, say N, consider only the N most recent events.}
+ \item{lambda}{function or matrix define the intensity of a Poisson process if
+ s.distr is Poisson.}
+ \item{lmax}{upper bound for the value of lambda.}
+ \item{nx, ny}{define the 2-D grid on which the intensity is evaluated if
+ \code{s.distr} is Poisson.}
+ \item{nt}{used to discretize time to compute the infection rate function.}
+ \item{t0}{minimum time used to compute the infection rate function.
+ Default is the minimum of \code{t.region}.}
+ \item{inhibition}{Logical. If TRUE, an inhibition process is generated.
+ Otherwise, it is a contagious process.}
+ \item{...}{additional parameters if \code{lambda} is a function.}
+}
+
+value{
+A list containing:
+\item{xyt}{matrix (or list of matrices if \code{nsim}>1)
+containing the points (x,y,t) of the simulated point pattern.}
+\item{s.region, t.region}{parameters passed in argument.}
+}
+
+\author{
+Edith Gabriel <edith.gabriel at univ-avignon.fr>, Peter J Diggle.
+}
+
+\seealso{
+ \code{\link{animation}} and \code{\link{stani}} for plotting space-time point patterns.
+ }
+
+\examples{
+# inhibition; spatial distribution: uniform
+inf1 = rinfec(npoints=100, alpha=0.2, beta=0.6, gamma=0.5, maxrad=c(0.075,0.5), t.region=c(0,50), s.distr="uniform", t.distr="uniform", h="step", p="min", recent="all", inhibition=TRUE)
+animation(inf1$xyt, cex=0.8, runtime=10)
+
+# contagion; spatial distribution: Poisson with intensity a given matrix
+data(fmd)
+data(northcumbria)
+h = mse2d(as.points(fmd[,1:2]), northcumbria, nsmse=30, range=3000)
+h = h$h[which.min(h$mse)]
+Ls = kernel2d(as.points(fmd[,1:2]), northcumbria, h, nx=50, ny=50)
+inf2 = rinfec(npoints=100, alpha=4, beta=0.6, gamma=20, maxrad=c(12000,20), s.region=northcumbria, t.region=c(1,2000), s.distr="poisson", t.distr="uniform", h="step", p="min", recent=1, lambda=Ls$z, inhibition=FALSE)
+image(Ls$x, Ls$y, Ls$z, col=grey((1000:1)/1000)); polygon(northcumbria,lwd=2)
+animation(inf2$xyt, add=TRUE, cex=0.7, runtime=15)
+}
\ No newline at end of file
Modified: pkg/man/rinter.Rd
===================================================================
--- pkg/man/rinter.Rd 2008-10-31 10:42:57 UTC (rev 21)
+++ pkg/man/rinter.Rd 2008-11-01 13:52:03 UTC (rev 22)
@@ -30,88 +30,21 @@
\eqn{\leq}{<=} 1.
Otherwise, \code{h} is monotone, decreasing, and must tend
to 1 when the distance tends to 0.}
- \item{gs, gt}{Must be choosen among "\code{min}", "\code{max}" and
- "\code{prod}". }
\item{thetas, thetat}{Parameters of \code{hs} and \code{ht} functions.}
\item{deltas, deltat}{Spatial and temporal distance of inhibition.}
+ \item{gs, gt}{Compute the probability of acceptance of a new point from
+ \code{hs} or \code{ht} and \code{recent}. Must be choosen among "\code{min}",
+ "\code{max}" and "\code{prod}". }
\item{recent}{If ``\code{all}'' consider all previous events. If
is an integer, say N, consider only the N most recent events.}
\item{nsim}{number of simulations to generate. Default is 1. }
- \item{discrete.time}{if TRUE, times are discrete, otherwise belong to R+.}
+ \item{discrete.time}{if TRUE, times belong to \eqn{\mathbb{N}}{N},
+ otherwise belong to \eqn{\mathbb{R}^+}{R^+}.}
\item{replace}{Logical. If TRUE allows times repeat.}
\item{inhibition}{Logical. If TRUE, an inhibition process is
generated. Otherwise, it is a contagious process. }
}
-
-#\subsubsection*{Details}
-#
-#Inhibition process ($k$th step):
-#\begin{enumerate}
-#\item Generate uniformly a location $s \in S$ and a time $t \in T$.
-#\item Generate $u_s \sim {\cal U}[0,1]$ and $u_t \sim {\cal U}[0,1]$.
-#\item If $\| s -s_j \| \geq \delta_s$ for all $j=1,\dots,k-1$,
-# then set $p_s = 1$.
-#
-# Otherwise compute $v_s = p_s \left( h_s
-# \left( (\| s -s_j \|)_{j=l,\dots,k-1}, \theta_s, \delta_s
-# \right), r_k \right)$
-# \item If $| t -t_j | \geq \delta_t$ for all $j=1,\dots,k-1$,
-# then set $p_t = 1$.
-#
-# Otherwise compute $v_t = p_t \left( h_t \left( (| t -t_j |)_{j=l,\dots,
-# k-1}, \theta_t, \delta_t \right), r_k \right)$
-# \item If $u_s < v_s$ and $u_t < v_t$, then keep $s$ and $t$.
-#\end{enumerate}
-#The functions $p_s$ and $p_t$ can be choosen among ``min'', ``max''
-#and ``prod'' and depend on the parameter $r_k$. This parameter allows us
-#to consider either all previous events or only the most recent in the
-#conditions 3 et 4.
-#The functions $h_s$ and $h_t$ depend on the distance between locations
-#and times respectively. They are monotone, increasing, tend
-#to 1 when the distance tends to infinity and satisfy $0 \leq h(\cdot) \leq 1$.
-#Currently the following functions are implemented:
-#\begin{itemize}
-#\item step: $h(x) = \lce
-# \begin{array}[l]{ll}
-# 1, \text{ if } x>\delta \\
-# \theta, \text{ otherwise}
-# \end{array} \right., \ \theta \in [0,1]
-# $
-#\item exponential: $h(x) = \exp(\theta x) / \max(\exp(\theta x)), \ \theta
-# \geq 0$,
-#\item gaussian: $h(x) = \exp(\theta x^2) / \max(\exp(\theta x^2)), \ \theta
-# \geq 0$.
-#\end{itemize}
-#Note that the user can call his own functions $h_s$ and $h_t$, provided
-#that they have (spatial or temporal distance between points),
-#\verb#theta# and \verb#delta# (even if not used by the function).
-#
-#{\bf *** should we also allow an user defined $p$ function?}
-#
-#\noindent
-#The simple sequential inhibition process is given using ``min'' for
-#$g$ functions, ``step'' for $h$ functions and 0 for $\theta$ parameter.
-#
-#Contagious process ($k$th step):
-#\begin{enumerate}
-#\item Generate uniformly a location $s \in S$ and a time $t \in T$.
-#\item Generate $u_s \sim {\cal U}[0,1]$ and $u_t \sim {\cal U}[0,1]$.
-#\item If $\| s_{k-1} -s \| < \delta_s$, then set $p_s = 1$.
-# Otherwise, compute $p_s =h_s \left( \| s_{k-1} -s \|, \theta_s,
-# \delta_s \right)$.
-#\item If $| t_{k-1} - t | < \delta_t$, then set $p_t=1$.
-# Otherwise, compute $p_t = h_t \left( | t_{k-1} - t|, \theta_t,
-# \delta_t \right)$.
-#\item If $u_s < p_s$ and $u_t < p_t$, then keep $s$ and $t$.
-#\end{enumerate}
-#The functions $h_s$ and $h_t$ depend on the distance between locations
-#and times respectively. They are monotone, decreasing, tend
-#to 1 when the distance tends to 0 and satisfy $0 \leq h(\cdot) \leq 1$.
-#$h_{cont.} (x,\theta,\delta)=h_{inhib.}(\max(x)-x,\theta,\max(x)-\delta)$.
-#The simple contagious process is given using ``step'' for $h$ functions
-#and 0 for $\theta$ parameters.
-
value{
A list containing:
\item{xyt}{matrix (or list of matrices if \code{nsim}>1)
@@ -159,18 +92,18 @@
}
return(res)
}
+d=seq(0,1,length=100)
+plot(d,hs(d,theta=0.2,delta=0.1,mus=0.1),xlab="",ylab="",type="l",ylim=c(0,1),lwd=2,las=1)
+lines(d,ht(d,theta=0.1,delta=0.05,mut=0.3),col=2,lwd=2)
+legend("bottomright",col=1:2,lty=1,lwd=2,legend=c(expression(h[s]),expression(h[t])),bty="n",cex=2)
-x=seq(0,1,length=100)
-plot(x,hs(x,theta=0.2,delta=0.1,mus=0.1),xlab="",ylab="",type="l",ylim=c(0,1))
-lines(x,ht(x,theta=0.1,delta=0.05,mut=0.3),col=2)
inh2 = rinter(npoints=100, hs=hs, gs="min", thetas=0.2, deltas=0.1, ht=ht, gt="min", thetat=0.1, deltat=0.05, inhibition=TRUE)
-animation(inh2$xyt,runtime=15,cex=0.8)
+animation(inh2$xyt, runtime=15, cex=0.8)
# simple contagious process for given spatial and temporal regions
data(northcumbria)
-cont1 = rinter(npoints=250, s.region=northcumbria, t.region=c(1,200), hs="gaussian",
- gs="min", thetas=1000, deltas=5000, ht="step", gt="min", thetat=0, deltat=10,
- recent=1, inhibition=FALSE)
+cont1 = rinter(npoints=250, s.region=northcumbria, t.region=c(1,200), thetas=0,
+ deltas=5000, thetat=0, deltat=10, recent=1, inhibition=FALSE)
animation(cont1$xyt, s.region=cont1$s.region, t.region=cont1$t.region,
incident="darkgreen", prevalent="green3", runtime=15, cex=0.8)
}
Added: pkg/man/rlgcp.Rd
===================================================================
--- pkg/man/rlgcp.Rd (rev 0)
+++ pkg/man/rlgcp.Rd 2008-11-01 13:52:03 UTC (rev 22)
@@ -0,0 +1,153 @@
+\name{rlgcp}
+\alias{rlgcp}
+\title{Generate log-Gaussian Cox point patterns}
+
+\description{
+Generate one (or several) realisation(s) of the log-Gaussian cox process
+in a region S x T.
+}
+
+\usage{
+ rlgcp(s.region, t.region, npoints=NULL, nsim=1, separable=TRUE, model="exponential",
+ param=c(1,1,1,1,1,2), var.grf=1, mean.grf=0, replace=TRUE, nx=100, ny=100,
+ nt=100, lmax=NULL, discrete.time=FALSE, exact=FALSE)
+}
+
+\arguments{
+ \item{s.region}{two-column matrix specifying polygonal region containing
+ all data locations. If \code{s.region} is missing, the unit square is considered.}
+ \item{t.region}{vector containing the minimum and maximum values of
+ the time interval. If \code{t.region} is missing, the interval [0,1] is considered.}
+ \item{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 \code{s.region} and \code{t.region}. }
+ \item{nsim}{number of simulations to generate. Default is 1.}
+ \item{separable}{Logical. If TRUE, the covariance function of the
+ Gaussian random field is separable. }
+ \item{model}{vector of length 1 or 2 specifying the model(s) of
+ covariance of the Gaussian random field. If \code{separable=TRUE} and
+ \code{model} is of length 2, then the elements of \code{model} define the
+ spatial and temporal covariances respectively.
+ If \code{separable=TRUE} and \code{model} is of length 1, then the
+ spatial and temporal covariances belongs to the same class of covariances,
+ among "exponential", "stable", "cauchy" and "wave" (see Details).
+ If \code{separable=FALSE}, \code{model} must be of length 1 and is either
+ "gneiting" or "cesare" (see Details).}
+ \item{param}{vector of parameters of the covariance function (see
+ Details).}
+ \item{var.grf}{variance of the Gaussian random field.}
+ \item{mean.grf}{mean of the Gaussian random field.}
+ \item{replace}{logical allowing times repeat.}
+ \item{nx,ny,nt}{define the size of the 3-D grid on which the intensity
+ is evaluated.}
+ \item{lmax}{upper bound for the value of lambda(x,y,t).}
+ \item{discrete.time}{if TRUE, times belong to \eqn{\mathbb{N}}{N},
+ otherwise belong to \eqn{\mathbb{R}^+}{R^+}.}
+}
+
+\details{
+We implemented stationary, isotropic spatio-temporal covariance functions.
+
+\emph{Separable covariance functions}
+
+\deqn{c(h,t) = c_s(\| h \|) \, c_t(|t|) , h \in S, t \in T}{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:
+\itemize{
+\item Exponential: \eqn{c(r) = \exp(-r)}{c(r)=exp(-r)},
+\item Stable: \eqn{c(r) = \exp(-r^\alpha)}{c(r)=exp(-r \alpha)},
+\eqn{\alpha \in [0,2]}{\alpha \in [0,2]},
+\item Cauchy: \eqn{c(r) = (1+r^2)^{-\alpha}}{c(r)=(1+r^2)^(-\alpha)},
+\eqn{\alpha >0}{\alpha > 0},
+\item Wave: \eqn{c(r) = \frac{\sin(r)}{r}}{c(r)=sin(r)/r} if \eqn{r>0}{r>0},
+ \eqn{c(0)=1}{c(0)=1},
+\item Matèrn: \eqn{c(r) = \frac{(r/\beta)^\alpha}{2^{\alpha-1}\Gamma(\alpha)}
+ {\cal K}_{\alpha}(r/\beta)}{c(r)={(r/\beta)^\alpha}/{2^(\alpha-1) \Gamma(\alpha)}
+ K_\alpha(r/\beta)}, \eqn{\alpha > 0}{\alpha>0} and \eqn{\beta > 0}{\beta>0}.
+
+\eqn{{\cal K}_{\alpha}}{K_\alpha} is the modified Bessel function of second kind:
+\deqn{{\cal K}_{\alpha}(x) = \frac{\pi}{2} \frac{I_{-\alpha}(x) -
+ I_{\alpha}(x)}{\sin(\pi \alpha)},}{K_\alpha(x) = (\pi/2) * (I_(-\alpha)(x)
+ - I_\alpha(x))/sin(\pi \alpha),} with
+ \eqn{I_{\alpha}(x) = \left( \frac{x}{2} \right)^{\alpha} \sum_{k=0}^\infty
+ \frac{1}{k! \Gamma(\alpha+k+1)} \left( \frac{x}{2} \right)^{2k}}{I_\alpha(x)
+ = (x/2)^\alpha sum_{k=0}^\infty 1/(k! \Gamma(\alpha+k+1)) (x/2)^(2k)}
+The parameters \eqn{alpha_1}{\alpha_1} and \eqn{\alpha_2}{\alpha_2} correspond
+to the parameters of the spatial and temporal covariance respectively.
+}
+
+
+\emph{Non-separable covariance functions}
+
+The spatio-temporal covariance function can be:
+\itemize{
+\item gneiting: \eqn{c(h,t) = \psi (t/\beta_2)^{-\alpha_6} \phi \left(\frac{h/
+ \beta_1}{\psi (t/\beta_2)} \right)}{c(h,t)=\psi(t/\beta_2)^(-\alpha_6)
+ \phi( (h/\beta_1)/\psi(t/\beta_2) )}, \eqn{\beta_1, \beta_2 >0}{\beta_1,
+ \beta_2 >0},
+ \itemize{
+ \item If \eqn{\alpha_2=1}{\alpha_2=1}, \eqn{\phi(r)}{\phi(r)} is the Stable model.
+ \item if \eqn{\alpha_2=2}{\alpha_2=2}, \eqn{\phi(r)}{\phi(r)} is the Cauchy model.
+ \item If \eqn{\alpha_2=3}{\alpha_2=3}, \eqn{\phi(r)}{\phi(r)} is the Matèrn model.
+ \item If \eqn{\alpha_5=1}{\alpha_5=1}, \eqn{\psi^2(r) = (r^{\alpha_3}+
+ 1)^{\alpha_4}}{\psi^2(r) = (r^\alpha_3+1)^\alpha_4},
+ \item If \eqn{\alpha_5=2}{\alpha_5=2}, \eqn{\psi^2(r) = (\alpha_4^{-1} r^{\alpha_3}
+ +1)/(r^{\alpha_3}+1)}{\psi^2(r) = (\alpha_4^(-1) r^\alpha_3 +1)/(r^\alpha_3+1)},
+ \item If \eqn{\alpha_5=3}{\alpha_5=3}, \eqn{\psi^2(r) = -\log(r^{\alpha_3} +
+ 1/{\alpha_4})/\log {\alpha_4}}{\psi^2(r) = -log(r^\alpha_3+1/\alpha_4)/
+ log(\alpha_4)},
+ }
+
+ The parameter \eqn{\alpha_1}{\alpha_1} is the respective parameter for the model of
+ \eqn{\phi(\cdot)}{\phi(.)}, \eqn{\alpha_3 \in (0,2]}{\alpha_3 in (0,2]},
+ \eqn{\alpha_4 \in (0,1]}{\alpha_4 in (0,1]} and \eqn{\alpha_6 \geq 2}{\alpha_6 >= 2}.
+\item cesare: \eqn{c(h,t) = \left( 1 + (h/\beta_1)^{\alpha_1} +
+ (t/\beta_2)^{\alpha_2} \right)^{-\alpha_3}}{c(h,t) = (1 + (h/\beta_1)^\alpha_1 +
+ (t/\beta_2)^\alpha_2)^(-\alpha_3)}, \eqn{\beta_1, \beta_2 >0}{\beta_1,
+ \beta_2 >0}, \eqn{\alpha_1, \alpha_2 \in [1,2]}{\alpha_1, \alpha_2 in [1,2]}
+ and \eqn{\alpha_3 \geq 3/2}{\alpha_3 >= 3/2}.
+ }
+}
+
+\value{
+A list containing:
+\item{xyt}{matrix (or list of matrices if \code{nsim}>1)
+containing the points (x,y,t) of the simulated point pattern.}
+\item{s.region, t.region}{parameters passed in argument.}
+\item{Lambda}{nx * ny * nt array (or list of array if \code{nsim}>1)
+of the intensity.}
+}
+
+\references{
+Chan, G. and Wood A. (1997). An algorithm for simulating stationary
+Gaussian random fields. Applied Statistics, Algorithm Section, 46, 171--181.
+
+Gneiting T. (2002). Nonseparable, stationary covariance functions
+for space-time data. Journal of the American Statistical Association,
+97, 590--600.
+}
+
+\author{
+Edith Gabriel <edith.gabriel at univ-avignon.fr>, Peter J Diggle.
+}
+
+\seealso{
+ \code{\link{animation}} and \code{\link{stani}} 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")
+}
Modified: pkg/man/rpp.Rd
===================================================================
--- pkg/man/rpp.Rd 2008-10-31 10:42:57 UTC (rev 21)
+++ pkg/man/rpp.Rd 2008-11-01 13:52:03 UTC (rev 22)
@@ -29,8 +29,8 @@
number of points is from a
Poisson distribution with mean the double integral of \code{lambda} over
\code{s.region} and \code{t.region}.}
- \item{discrete.time}{if \code{TRUE}, times belong to N,
- otherwise belong to R+. }
+ \item{discrete.time}{if TRUE, times belong to \eqn{\mathbb{N}}{N},
+ otherwise belong to \eqn{\mathbb{R}^+}{R^+}.}
\item{nsim}{number of simulations to generate. Default is 1.}
\item{nx,ny,nt}{define the size of the 3-D grid on which the intensity
is evaluated. }
@@ -38,6 +38,7 @@
\code{lambda} is a function.}
\item{Lambda}{matrix of spatial intensity if \code{lambda} is a character.}
\item{mut}{vector of temporal intensity if \code{lambda} is a character.}
+ \item{...}{additional parameters if \code{lambda} is a function.}
}
\value{
More information about the Stpp-commits
mailing list