[Stpp-commits] r21 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 31 11:42:58 CET 2008
Author: gabriele
Date: 2008-10-31 11:42:57 +0100 (Fri, 31 Oct 2008)
New Revision: 21
Added:
pkg/R/rinter.r
pkg/man/rinter.Rd
pkg/man/rpcp.Rd
pkg/man/rpp.Rd
Removed:
pkg/R/pcp.weight.r
pkg/R/rinhib.r
Modified:
pkg/R/animation.r
pkg/R/pcp.larger.region.r
pkg/R/rpcp.r
pkg/R/rpp.r
pkg/R/spatial.inhibition.r
pkg/R/temporal.inhibition.r
pkg/man/animation.Rd
Log:
Modified: pkg/R/animation.r
===================================================================
--- pkg/R/animation.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/animation.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,4 +1,4 @@
-animation <- function(xyt, s.region, t.region, runtime=1, incident="red", prevalent="pink3", pch=19, cex=0.35, plot.s.region=T, scales=T, border.frac=0.05, add=F)
+animation <- function(xyt, s.region, t.region, runtime=1, incident="red", prevalent="pink3", pch=19, cex=0.5, plot.s.region=T, scales=T, border.frac=0.05, add=F)
{
#
# Description:
@@ -38,7 +38,7 @@
xy<-as.matrix(sxyt[,1:2])
tt<-sxyt[,3]
npts<-length(tt)
- T0=length(unique(tt))
+ T0 <- max(t.region)
if (add==F)
{
Modified: pkg/R/pcp.larger.region.r
===================================================================
--- pkg/R/pcp.larger.region.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/pcp.larger.region.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,194 +1,196 @@
-pcp.larger.region <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infecD=TRUE, maxradlarger, ...)
-{
- if (missing(cluster)) cluster <- "uniform"
-
- 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(maxrad)) maxrad <- c(0.05,0.05)
- if (length(maxrad)==1) maxrad=rep(maxrad,2)
- maxrads <- maxrad[1]
- maxradt <- maxrad[2]
-
- if (missing(maxradlarger)) maxradlarger <- maxrad
- maxradlargers <- maxradlarger[1]
- maxradlargert <- maxradlarger[2]
-
- if (length(cluster)==1)
- {
- s.distr <- cluster
- t.distr <- cluster
- }
- else
- {
- s.distr <- cluster[1]
- t.distr <- cluster[2]
- }
-
- t.region <- sort(t.region)
- s.area <- areapl(s.region)
- t.area <- t.region[2]-t.region[1]
-
- if (maxradlargers==0)
- s.larger <- s.region
- else
- {
- s.larger <- rbind(s.region+maxradlargers,s.region-maxradlargers,cbind(s.region[,1]+maxradlargers,s.region[,2]-maxradlargers),cbind(s.region[,1]-maxradlargers,s.region[,2]+maxradlargers))
- M <- chull(s.larger)
- s.larger <- cbind(s.larger[M,1],s.larger[M,2])
- }
- t.larger <- c(t.region[1]-maxradlargert,t.region[2]+maxradlargert)
-
- if (t.larger[1]<0) t.larger[1] <- 0
-
- if (is.null(nparents))
- {
- if (is.null(lambda)) stop("please specify either the number of parents or the intensity of parents process")
- if (is.function(lambda)) stop("please specify the number of parents")
- if (is.numeric(lambda))
- lambda <- lambda*(areapl(s.larger)*(t.larger[2]-t.larger[1]))/(s.area*t.area)
- npar <- rpois(1,lambda)
- }
- else npar <- nparents
-
- if (is.null(lambda)) lambda <- (npar/(s.area*t.area))*(areapl(s.larger)*(t.larger[2]-t.larger[1]))/(s.area*t.area)
-
- if (is.null(npoints))
- {
- if (is.null(mc)) stop("please specify either the number of points to simulate or the mean number of children per parents")
- else
- {
- nchild <- matrix(rpois(npar,lambda=mc),npar,1)
- npts <- sum(nchild)
- }
- }
- else
- {
- npts <- npoints
- if (is.null(mc)) mc <- npts/npar
- nchild <- matrix(rpois(npar,lambda=mc),npar,1)
- nc <- sum(nchild)
- while (nc < npts)
- {
- nchild <- matrix(rpois(npar,lambda=mc),npar,1)
- nc <- sum(nchild)
- }
- }
-
- parpts <- rpp(lambda=lambda,s.region=s.larger,t.region=t.larger,npoints=npar,...)$xyt
-
- pattern.interm <- NULL
- ipt <- 1
- nchild2 <- rep(0,npar)
- for (ipar in 1:npar)
- {
- xpar <- parpts[ipar,1]
- ypar <- parpts[ipar,2]
- zpar <- parpts[ipar,3]
- nc <- nchild[ipar,1]
-
- if (s.distr=="uniform")
- {
- xp <- xpar+runif(nc,min=-maxrads,max=maxrads)
- yp <- ypar+runif(nc,min=-maxrads,max=maxrads)
- }
- if (s.distr=="normal")
- {
- xp <- xpar+rnorm(nc,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(nc,mean=0,sd=maxrads/2)
- }
- if (s.distr=="exponential")
- {
- xp <- xpar+rexp(nc,rate=1/maxrads)
- yp <- ypar+rexp(nc,rate=1/maxrads)
- }
- if (t.distr=="uniform")
- {
- if (infecD==TRUE)
- zp <- zpar+runif(nc,min=-maxradt,max=maxradt)
- else
- zp <- zpar+runif(nc,min=0,max=maxradt)
- }
- if (t.distr=="normal")
- {
- if (infecD==TRUE)
- zp <- zpar+abs(rnorm(nc,mean=0,sd=maxradt/2))
- else
- zp <- zpar+rnorm(nc,mean=0,sd=maxradt/2)
- }
- if (t.distr=="exponential")
- {
- if (infecD==TRUE)
- zp <- zpar+rexp(nc,rate=1/maxradt)
- else
- zp <- zpar+sample(c(-1,1),1)*rexp(nc,rate=1/maxradt)
- }
-
- mask <- ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1))
- nchild2[ipar] <- sum(mask)
- pattern.interm <- rbind(pattern.interm,cbind(xp[mask],yp[mask],zp[mask],rep(ipar,sum(mask))))
- }
-
- ipt <- dim(pattern.interm)[[1]]
- if (ipt > npts)
- {
- samp <- sample(1:ipt,npts)
- pattern.interm <- pattern.interm[samp,]
- }
- while(ipt < npts)
- {
- ipar <- sample(1:npar,1)
- nchild2[ipar] <- nchild2[ipar]+1
- xpar <- parpts[ipar,1]
- ypar <- parpts[ipar,2]
- zpar <- parpts[ipar,3]
-
- if (s.distr=="uniform")
- {
- xp <- xpar+runif(1,min=-maxrads,max=maxrads)
- yp <- ypar+runif(1,min=-maxrads,max=maxrads)
- }
- if (s.distr=="normal")
- {
- xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
- }
- if (s.distr=="exponential")
- {
- xp <- xpar+rexp(1,rate=1/maxrads)
- yp <- ypar+rexp(1,rate=1/maxrads)
- }
- if (t.distr=="uniform")
- {
- if (infecD==TRUE)
- zp <- zpar+runif(1,min=-maxradt,max=maxradt)
- else
- zp <- zpar+runif(1,min=0,max=maxradt)
- }
- if (t.distr=="normal")
- {
- if (infecD==TRUE)
- zp <- zpar+abs(rnorm(1,mean=0,sd=maxradt/2))
- else
- zp <- zpar+rnorm(1,mean=0,sd=maxradt/2)
- }
- if (t.distr=="exponential")
- {
- if (infecD==TRUE)
- zp <- zpar+rexp(1,rate=1/maxradt)
- else
- zp <- zpar+sample(c(-1,1),1)*rexp(1,rate=1/maxradt)
- }
-
- if ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1))
- {
- pattern.interm <- rbind(pattern.interm,cbind(xp,yp,zp,ipar))
- ipt <- ipt+1
- }
- }
-
- pts <- pattern.interm[,1:3]
- invisible(return(list(pts=pts,s.larger=s.larger,t.larger=t.larger,index.child=pattern.interm[,4],nchild=nchild2,parpts=parpts)))
-}
-
-
+pcp.larger.region <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infecD=TRUE, maxradlarger, ...)
+{
+ if (missing(cluster)) cluster <- "uniform"
+
+ 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(maxrad)) maxrad <- c(0.05,0.05)
+ if (length(maxrad)==1) maxrad=rep(maxrad,2)
+ maxrads <- maxrad[1]
+ maxradt <- maxrad[2]
+
+ if (missing(maxradlarger)) maxradlarger <- maxrad
+ maxradlargers <- maxradlarger[1]
+ maxradlargert <- maxradlarger[2]
+
+ if (length(cluster)==1)
+ {
+ s.distr <- cluster
+ t.distr <- cluster
+ }
+ else
+ {
+ s.distr <- cluster[1]
+ t.distr <- cluster[2]
+ }
+
+ t.region <- sort(t.region)
+ s.area <- areapl(s.region)
+ t.area <- t.region[2]-t.region[1]
+
+ if (maxradlargers==0)
+ s.larger <- s.region
+ else
+ {
+ s.larger <- rbind(s.region+maxradlargers,s.region-maxradlargers,cbind(s.region[,1]+maxradlargers,s.region[,2]-maxradlargers),cbind(s.region[,1]-maxradlargers,s.region[,2]+maxradlargers))
+ M <- chull(s.larger)
+ s.larger <- cbind(s.larger[M,1],s.larger[M,2])
+ }
+ t.larger <- c(t.region[1]-maxradlargert,t.region[2]+maxradlargert)
+
+ if (t.larger[1]<0) t.larger[1] <- 0
+
+ if (is.null(nparents))
+ {
+ if (is.null(lambda)) stop("please specify either the number of parents or the intensity of parents process")
+ if (is.function(lambda)) stop("please specify the number of parents")
+ if (is.numeric(lambda))
+ lambda <- lambda*(areapl(s.larger)*(t.larger[2]-t.larger[1]))/(s.area*t.area)
+ npar <- rpois(1,lambda)
+ }
+ else npar <- nparents
+
+ if (is.null(lambda)) lambda <- (npar/(s.area*t.area))*(areapl(s.larger)*(t.larger[2]-t.larger[1]))/(s.area*t.area)
+
+ if (is.null(npoints))
+ {
+ if (is.null(mc)) stop("please specify either the number of points to simulate or the mean number of children per parents")
+ else
+ {
+ nchild <- matrix(rpois(npar,lambda=mc),npar,1)
+ npts <- sum(nchild)
+ }
+ }
+ else
+ {
+ npts <- npoints
+ if (is.null(mc)) mc <- npts/npar
+ nchild <- matrix(rpois(npar,lambda=mc),npar,1)
+ nc <- sum(nchild)
+ while (nc < npts)
+ {
+ nchild <- matrix(rpois(npar,lambda=mc),npar,1)
+ nc <- sum(nchild)
+ }
+ }
+
+ parpts <- rpp(lambda=lambda,s.region=s.larger,t.region=t.larger,npoints=npar,...)$xyt
+
+ pattern.interm <- NULL
+ ipt <- 1
+ nchild2 <- rep(0,npar)
+ for (ipar in 1:npar)
+ {
+ xpar <- parpts[ipar,1]
+ ypar <- parpts[ipar,2]
+ zpar <- parpts[ipar,3]
+ nc <- nchild[ipar,1]
+
+ if (s.distr=="uniform")
+ {
+ xp <- xpar+runif(nc,min=-maxrads,max=maxrads)
+ yp <- ypar+runif(nc,min=-maxrads,max=maxrads)
+ }
+ if (s.distr=="normal")
+ {
+ xp <- xpar+rnorm(nc,mean=0,sd=maxrads/2)
+ yp <- ypar+rnorm(nc,mean=0,sd=maxrads/2)
+ }
+ if (s.distr=="exponential")
+ {
+ xp <- xpar+rexp(nc,rate=1/maxrads)
+ yp <- ypar+rexp(nc,rate=1/maxrads)
+ }
+ if (t.distr=="uniform")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+runif(nc,min=-maxradt,max=maxradt)
+ else
+ zp <- zpar+runif(nc,min=0,max=maxradt)
+ }
+ if (t.distr=="normal")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+abs(rnorm(nc,mean=0,sd=maxradt/2))
+ else
+ zp <- zpar+rnorm(nc,mean=0,sd=maxradt/2)
+ }
+ if (t.distr=="exponential")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+rexp(nc,rate=1/maxradt)
+ else
+ zp <- zpar+sample(c(-1,1),1)*rexp(nc,rate=1/maxradt)
+ }
+
+ mask <- ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1))
+ nchild2[ipar] <- sum(mask)
+ pattern.interm <- rbind(pattern.interm,cbind(xp[mask],yp[mask],zp[mask],rep(ipar,sum(mask))))
+ }
+
+ ipt <- dim(pattern.interm)[[1]]
+ if (ipt > npts)
+ {
+ samp <- sample(1:ipt,npts)
+ pattern.interm <- pattern.interm[samp,]
+ }
+ while(ipt < npts)
+ {
+ ipar <- sample(1:npar,1)
+ nchild2[ipar] <- nchild2[ipar]+1
+ xpar <- parpts[ipar,1]
+ ypar <- parpts[ipar,2]
+ zpar <- parpts[ipar,3]
+
+ if (s.distr=="uniform")
+ {
+ xp <- xpar+runif(1,min=-maxrads,max=maxrads)
+ yp <- ypar+runif(1,min=-maxrads,max=maxrads)
+ }
+ if (s.distr=="normal")
+ {
+ xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
+ yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
+ }
+ if (s.distr=="exponential")
+ {
+ xp <- xpar+rexp(1,rate=1/maxrads)
+ yp <- ypar+rexp(1,rate=1/maxrads)
+ }
+ if (t.distr=="uniform")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+runif(1,min=-maxradt,max=maxradt)
+ else
+ zp <- zpar+runif(1,min=0,max=maxradt)
+ }
+ if (t.distr=="normal")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+abs(rnorm(1,mean=0,sd=maxradt/2))
+ else
+ zp <- zpar+rnorm(1,mean=0,sd=maxradt/2)
+ }
+ if (t.distr=="exponential")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+rexp(1,rate=1/maxradt)
+ else
+ zp <- zpar+sample(c(-1,1),1)*rexp(1,rate=1/maxradt)
+ }
+
+ if ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1))
+ {
+ pattern.interm <- rbind(pattern.interm,cbind(xp,yp,zp,ipar))
+ ipt <- ipt+1
+ }
+ }
+
+ pts <- pattern.interm[,1:3]
+ ott<-order(pts[,3])
+ pts<-pts[ott,]
+ invisible(return(list(pts=pts,s.larger=s.larger,t.larger=t.larger,index.child=pattern.interm[,4],nchild=nchild2,parpts=parpts)))
+}
+
+
Deleted: pkg/R/pcp.weight.r
===================================================================
--- pkg/R/pcp.weight.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/pcp.weight.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,117 +0,0 @@
-pcp.weight <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infecD=TRUE, ...)
-{
- if (missing(cluster)) cluster <- "uniform"
-
- 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(maxrad)) maxrad <- c(0.05,0.05)
- maxrads <- maxrad[1]
- maxradt <- maxrad[2]
-
- t.region <- sort(t.region)
- s.area <- areapl(s.region)
- t.area <- t.region[2]-t.region[1]
-
- if (is.null(nparents))
- {
- if (is.null(lambda)) stop("please specify either the number of parents or the intensity of parents process")
- if (is.function(lambda)) stop("please specify the number of parents")
- npar <- rpois(1,lambda)
- }
- else npar <- nparents
- if (is.null(lambda)) lambda <- npar/(s.area*t.area)
-
- parpts <- rpp(lambda=lambda,s.region=s.region,t.region=t.region,npoints=npar,...)$pts
- npoly <- length(s.region[, 1])
- polyx <- c(s.region[, 1], s.region[1, 1])
- polyy <- c(s.region[, 2], s.region[1, 2])
- binft <- t.region[1]
- bsupt <- t.region[2]
- W <- rep(0,npar)
- cont <- as.numeric(infecD)
-
- weight <- .Fortran("edgecorrect", as.double(parpts[,1]),
- as.double(parpts[,2]), as.double(parpts[,3]),
- as.integer(npar), as.double(polyx),
- as.double(polyy), as.integer(npoly),
- as.double(maxrads), as.double(maxradt),
- as.double(binft), as.double(bsupt),
- as.integer(cont), as.double(W))
- W <- weight[[13]]
-
- if (is.null(npoints))
- {
- if (is.null(mc)) stop("please specify either the number of points to simulate or the mean number of children per parents")
- npts <- round(sum(rpois(npar,lambda=mc)*W),0)
- }
- else npts <- npoints
- pattern.interm <- matrix(0,npts,3)
- nchild <- matrix(0,npar,1)
- ipt <- 1
- while(ipt <= npts)
- {
- ipar <- sample(1:npar,1,prob=W)
- xpar <- parpts[ipar,1]
- ypar <- parpts[ipar,2]
- zpar <- parpts[ipar,3]
- nchild[ipar,1] <- nchild[ipar,1]+1
-
- if (length(cluster)==1)
- {
- s.distr <- cluster
- t.distr <- cluster
- }
- else
- {
- s.distr <- cluster[1]
- t.distr <- cluster[2]
- }
- outcirc <- 1
- while(outcirc == 1)
- {
- if (s.distr=="uniform")
- {
- xp <- xpar+runif(1,min=-maxrads,max=maxrads)
- yp <- ypar+runif(1,min=-maxrads,max=maxrads)
- }
- if (s.distr=="normal")
- {
- xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
- }
- if (s.distr=="exponential")
- {
- xp <- xpar+rexp(1,rate=1/maxrads)
- yp <- ypar+rexp(1,rate=1/maxrads)
- }
- if (t.distr=="uniform")
- {
- if (infecD==TRUE)
- zp <- zpar+runif(1,min=-maxradt,max=maxradt)
- else
- zp <- zpar+runif(1,min=0,max=maxradt)
- }
- if (t.distr=="normal")
- {
- if (infecD==TRUE)
- zp <- zpar+abs(rnorm(1,mean=0,sd=maxradt/2))
- else
- zp <- zpar+rnorm(1,mean=0,sd=maxradt/2)
- }
- if (t.distr=="exponential")
- {
- if (infecD==TRUE)
- zp <- zpar+rexp(1,rate=1/maxradt)
- else
- zp <- zpar+sample(c(-1,1),1)*rexp(1,rate=1/maxradt)
- }
-
- if ((inout(as.points(x=xp,y=yp),poly=s.region)==TRUE) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/maxrads^2) + (((ypar-yp)^2)/maxrads^2)) < 1)) outcirc <- 0
- }
- pattern.interm[ipt,1:3] <- cbind(xp,yp,zp)
- ipt <- ipt+1
- }
- pts <- pattern.interm
-}
-
Deleted: pkg/R/rinhib.r
===================================================================
--- pkg/R/rinhib.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/rinhib.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,80 +0,0 @@
-rinhib <- function(npoints,s.region,t.region,hs="step",ps="min",thetas=0,deltas,ht="step",pt="min",thetat=1,deltat,recent="all",nsim=1,discrete.time=FALSE,replace=FALSE,inhibition=TRUE,...)
- {
- #
- # Simulate an inhibition or a contagious 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.
- #
- # h: a function of the distance between points and theta.
- # If inhibition=TRUE, h is monotone, increasing, and must tend
- # to 1 when the distance tends to infinity. 0<= h(d,theta) <= 1.
- # Else, h is monotone, decreasing, and must tend
- # to 1 when the distance tends to 0. 0<= h(d,theta) <= 1.
- #
- # p: a function among "min", "max", "prod".
- #
- # replace: logical allowing times repetition.
- #
- # npoints: number of points to simulate.
- #
- # nsim: number of simulations to generate. Default is 1.
- #
- #
- # Value:
- # pattern: list containing the points (x,y,t) of the simulated point
- # process.
- #
- ##
- ## E. GABRIEL, 09/03/2007
- ##
- ## last modification: 13/03/2007
- ## 27/03/2007
- ##
- ##
-
- 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 (!(is.function(hs)))
- {
- if ((inhibition==TRUE) && (thetas==0) && (hs=="step") && (npoints * pi * deltas^2/4 > areapl(s.region)))
- stop(paste("s.region is too small to fit", npoints, "points", "at minimum distance", deltas))
-
- if ((inhibition==TRUE) && (thetas==0) && (hs=="step") && ((max(t.region)-min(t.region))/deltat<npoints))
- stop(paste("t.region is too small to fit", npoints, "points", "at minimum time interval", deltat))
- }
-
- pattern <- list()
- ni <- 1
- while(ni<=nsim)
- {
- pattern.interm <- spatial.inhibition(npoints,h=hs,theta=thetas,delta=deltas,p=ps,recent=recent,s.region=s.region,inhibition=inhibition,...)$pts
- times.interm <- temporal.inhibition(npoints,h=ht,theta=thetat,delta=deltat,p=pt,recent=recent,t.region=t.region,inhibition=inhibition,discrete.time=discrete.time,replace=replace,...)$times
-
- if (nsim==1)
- {
- pattern <- cbind(pattern.interm,times.interm)
- ni <- ni+1
- }
- else
- {
- pattern[[ni]] <- cbind(pattern.interm,times.interm)
- ni <- ni+1
- }
- }
- invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region)))
-}
-
-
-
-
Copied: pkg/R/rinter.r (from rev 18, pkg/R/rinhib.r)
===================================================================
--- pkg/R/rinter.r (rev 0)
+++ pkg/R/rinter.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -0,0 +1,74 @@
+rinter <- function(npoints,s.region,t.region,hs="step",gs="min",thetas=0,deltas,ht="step",gt="min",thetat=1,deltat,recent="all",nsim=1,discrete.time=FALSE,replace=FALSE,inhibition=TRUE,...)
+ {
+ #
+ # Simulate an inhibition or a contagious 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.
+ #
+ # hs, ht: a function of the distance between points and theta.
+ # If inhibition=TRUE, h is monotone, increasing, and must tend
+ # to 1 when the distance tends to infinity. 0<= h(d,theta) <= 1.
+ # Else, h is monotone, decreasing, and must tend
+ # to 1 when the distance tends to 0. 0<= h(d,theta) <= 1.
+ #
+ # gs, gt: a function among "min", "max", "prod".
+ #
+ # replace: logical allowing times repetition.
+ #
+ # npoints: number of points to simulate.
+ #
+ # nsim: number of simulations to generate. Default is 1.
+ #
+ #
+ # Value:
+ # xyt: list containing the points (x,y,t) of the simulated point
+ # process.
+ #
+
+
+ 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 (!(is.function(hs)))
+ {
+ if ((inhibition==TRUE) && (thetas==0) && (hs=="step") && (npoints * pi * deltas^2/4 > areapl(s.region)))
+ stop(paste("s.region is too small to fit", npoints, "points", "at minimum distance", deltas))
+
+ if ((inhibition==TRUE) && (thetas==0) && (hs=="step") && ((max(t.region)-min(t.region))/deltat<npoints))
+ stop(paste("t.region is too small to fit", npoints, "points", "at minimum time interval", deltat))
+ }
+
+ pattern <- list()
+ ni <- 1
+ while(ni<=nsim)
+ {
+ pattern.interm <- spatial.inhibition(npoints,h=hs,theta=thetas,delta=deltas,p=gs,recent=recent,s.region=s.region,inhibition=inhibition,...)$pts
+ times.interm <- temporal.inhibition(npoints,h=ht,theta=thetat,delta=deltat,p=gt,recent=recent,t.region=t.region,inhibition=inhibition,discrete.time=discrete.time,replace=replace,...)$times
+
+ if (nsim==1)
+ {
+ pattern <- cbind(pattern.interm,times.interm)
+ ni <- ni+1
+ }
+ else
+ {
+ pattern[[ni]] <- cbind(pattern.interm,times.interm)
+ ni <- ni+1
+ }
+ }
+ invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region)))
+}
+
+
+
+
Property changes on: pkg/R/rinter.r
___________________________________________________________________
Name: svn:mergeinfo
+
Modified: pkg/R/rpcp.r
===================================================================
--- pkg/R/rpcp.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/rpcp.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,4 +1,4 @@
-rpcp <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, contagious=TRUE, edge = "larger.region", ...)
+rpcp <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infectious=TRUE, edge = "larger.region", ...)
{
#
# Simulate a space-time Poisson cluster process in a region D x T.
@@ -49,7 +49,7 @@
# to their parent, such that children lies in the 95% IC
# of the normal distribution.
#
- # contagious: If TRUE (default), corresponds to contagious diseases
+ # infectious: If TRUE (default), corresponds to infectious diseases
# (times of children are always greater than parent's time).
#
# edge: specify the edge correction to use, "weight", "larger.region",
@@ -58,16 +58,10 @@
# ...: additional parameters of the intensity of the parent process.
#
# Value:
- # pts: matrix (or list of matrix if nsim>1) containing the points (x,y,t)
+ # xyt: matrix (or list of matrix if nsim>1) containing the points (x,y,t)
# of the simulated point process.
#
- ##
- ## E. GABRIEL, 09/10/2006
- ##
- ## last modification:
- ##
-# library(splancs)
if (missing(cluster)) cluster <- "uniform"
@@ -98,14 +92,11 @@
while(ni<=nsim)
{
- if (edge=="weight")
- pattern.interm <- pcp.weight(s.region=s.region, t.region=t.region, nparents=nparents, npoints=npoints, lambda=lambda, mc=mc, cluster=cluster, maxrad=maxrad, infecD=contagious, ...)
-
if (edge=="larger.region")
- pattern.interm <- pcp.larger.region(s.region=s.region, t.region=t.region, nparents=nparents, npoints=npoints, lambda=lambda, mc=mc, cluster=cluster, maxrad=maxrad, infecD=contagious, ...)$pts
+ pattern.interm <- pcp.larger.region(s.region=s.region, t.region=t.region, nparents=nparents, npoints=npoints, lambda=lambda, mc=mc, cluster=cluster, maxrad=maxrad, infecD=infectious, ...)$pts
if (edge=="without")
- pattern.interm <- pcp.larger.region(s.region=s.region, t.region=t.region, nparents=nparents, npoints=npoints, lambda=lambda, mc=mc, cluster=cluster, maxrad=maxrad, infecD=contagious, maxradlarger=c(0,0), ...)$pts
+ pattern.interm <- pcp.larger.region(s.region=s.region, t.region=t.region, nparents=nparents, npoints=npoints, lambda=lambda, mc=mc, cluster=cluster, maxrad=maxrad, infecD=infectious, maxradlarger=c(0,0), ...)$pts
if (nsim==1)
pattern <- pattern.interm
Modified: pkg/R/rpp.r
===================================================================
--- pkg/R/rpp.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/rpp.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,4 +1,4 @@
-rpp <- function(lambda, s.region, t.region, npoints=NULL, nsim=1, replace=T, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, Lambda=NULL, mut=NULL, ...)
+rpp <- function(lambda, s.region, t.region, npoints=NULL, nsim=1, replace=TRUE, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, Lambda=NULL, mut=NULL, ...)
{
#
# Simulate a space-time Poisson process in a region D x T.
@@ -8,7 +8,11 @@
# Arguments:
#
# lambda: Spatio-temporal intensity of the Poisson process.
- # Either a single positive number or a function(x,y,t,...)
+ # If 'lambda' is a single positive number, the function
+ # generates realisations of a homogeneous Poisson process,
+ # whilst if 'lambda' is a function of the form
+ # lambda(x,y,t,...) or a character it generates
+ # realisations of an inhomogeneous Poisson process.
#
# s.region: two columns matrix specifying polygonal region containing
# all data locations. If s.region is missing, the unit square
@@ -34,8 +38,13 @@
# lmax: upper bound for the value of lambda(x,y,t), if lambda
# is a function.
#
+ # Lambda: matrix of spatial intensity if 'lambda' is a character.
+ #
+ # mut: vector of temporal intensity if 'lambda' is a character.
+ #
+ #
# Value:
- # pts: matrix (or list if nsim>1) containing the points (x,y,t)
+ # xyt: matrix (or list if nsim>1) containing the points (x,y,t)
# of the simulated point process.
# Lambda: an array of the intensity surface at each time.
#
@@ -43,14 +52,7 @@
# is:
# p(s,t)=lambda(s,t)/(max_{(s_i,t_i), i=1,...,ngrid} lambda(s_i,t_i)).
#
- ##
- ## E. GABRIEL, 28/12/2005
- ##
- ## last modification: 06/10/2006
- ##
-# library(splancs)
-
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)
Modified: pkg/R/spatial.inhibition.r
===================================================================
--- pkg/R/spatial.inhibition.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/spatial.inhibition.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,219 +1,170 @@
-spatial.inhibition <- function(npoints,h,theta,delta,p,recent="all",s.region,inhibition=TRUE,...)
- {
- #
- # Simulate an inhibition or a contagious spatial point process in a region D
- #
- # Requires Splancs.
- #
- # Arguments:
- #
- # s.region: two columns matrix specifying polygonal region containing
- # all data locations. If poly is missing, the unit square is
- # considered.
- #
- # h: a function of the distance between points.
- # If inhibition=TRUE, h is monotone, increasing, and must tend
- # to 1 when the distance tends to infinity. 0<= h(d,theta) <= 1.
- # Else, h is monotone, decreasing, and must tend
- # to 1 when the distance tends to 0. 0<= h(d,theta) <= 1.
- #
- # p: a function among "min", "max", "prod".
- #
- # npoints: number of points to simulate.
- #
- #
- # Value:
- # pattern: list containing the points (x,y) of the simulated point
- # process.
- #
- ##
- ## E. GABRIEL, 26/03/2007
- ##
- ## last modification: 21/10/2008
- ##
- ##
-
- if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
-
- if ((inhibition==TRUE) && (npoints * pi * delta[1]^2/4 > areapl(s.region)))
- stop(paste("s.region is too small to fit", npoints, "points", "at minimum distance", delta[1]))
-
- if (!(is.function(h)))
- {
- models <- c("step","exponential","gaussian")
- if (sum(h==models)==0)
- {
- message <- paste("this model is not implemented, please choose among: ",paste(models,"",sep=" ",collapse="and "))
- stop(message)
- }
- if (h=="step")
- {
- hk <- function(d,theta,delta)
- {
- res <- rep(1,length(d))
- res[d<=delta] <- theta
- return(res)
- }
- }
-
- if (h=="exponential")
- {
- hk <- function(d,theta,delta)
- {
- res <- exp(d*theta)/max(exp(d*theta))
- return(res)
- }
- }
- if (h=="gaussian")
- {
- hk <- function(d,theta,delta)
- {
- res <- exp((d^2)*theta)/max(exp((d^2)*theta))
- return(res)
- }
- }
- }
-
- if (inhibition==FALSE)
- {
- if (!(is.function(h)))
- {
- hh <- hk
- hk <- function(d,theta,delta)
- {
- res <- hh(max(d)-d,theta,max(d)-delta)
- return(res)
- }
- }
- else
- {
- hh <- h
- hk <- function(d,theta,delta)
- {
- res <- h(max(d)-d,theta,max(d)-delta,...)
- return(res)
- }
- }
- }
-
- if (p=="min")
- {
- pk <- function(d,h,recent,theta,delta,...)
- {
- if (recent=="all")
- res <- min(h(d=d,theta=theta,delta=delta,...))
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(d))
- res <- min(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta,...))
- else
- res <- min(h(d=d,theta=theta,delta=delta,...))
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
-
- if (p=="max")
- {
- pk <- function(d,h,recent,theta,delta,...)
- {
- if (recent=="all")
- res <- max(h(d=d,theta=theta,delta=delta,...))
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(d))
- res <- max(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta,...))
- else
- res <- max(h(d=d,theta=theta,delta=delta,...))
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
-
- if (p=="prod")
- {
- pk <- function(d,h,recent,theta,delta,...)
- {
- if (is.numeric(recent) && recent<=length(d))
- res <- prod(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta,...))
- else
- res <- prod(h(d=d,theta=theta,delta=delta,...))
- return(res)
- }
- }
-
- xy <- csr(n=1,poly=s.region)
- npts <- 1
- pattern.interm <- cbind(x=xy[1],y=xy[2])
- if (inhibition==TRUE)
- {
- while(npts < npoints)
- {
- prob <- runif(1)
-
- xy <- csr(n=1,poly=s.region)
- if (all((sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2)) > delta))
- umax <- 1
- else
- {
- if (is.function(h))
- umax <- pk(d=sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2),h,recent,theta,delta,...)
- else
- {
- h <- hk
- umax <- pk(d=sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2),h,recent,theta,delta,...)
- }
- }
- if (prob<umax)
- {
- pattern.interm <- rbind(pattern.interm,c(xy[1],xy[2]))
- npts <- npts+1
- }
- }
- }
- else
- {
- while(npts < npoints)
- {
- prob <- runif(1)
-
- continue <- FALSE
- while(continue==FALSE)
- {
- xy <- csr(n=1,poly=s.region)
- if (sqrt((xy[1] - pattern.interm[npts,1])^2 + (xy[2] - pattern.interm[npts,2])^2) < delta)
- umax <- 1
- else
- {
- if (is.function(h))
- umax <- pk(d=sqrt((xy[1] - pattern.interm[npts,1])^2 + (xy[2] - pattern.interm[npts,2])^2),h,recent,theta,delta,...)
- else
- {
- h <- hk
- umax <- pk(d=sqrt((xy[1] - pattern.interm[npts,1])^2 + (xy[2] - pattern.interm[npts,2])^2),h,recent,theta,delta,...)
- }
- }
- if (prob < umax)
- {
- pattern.interm <- rbind(pattern.interm,c(xy[1],xy[2]))
- npts <- npts+1
- continue <- TRUE
- }
- }
- }
- }
-
- invisible(return(list(pts=pattern.interm,s.region=s.region)))
-}
-
-
+spatial.inhibition <- function(npoints,h,theta,delta,p,recent="all",s.region,inhibition=TRUE,...)
+ {
+ #
+ # Simulate an inhibition or a contagious spatial point process in a region D
+ #
+ # Requires Splancs.
+ #
+ # Arguments:
+ #
+ # s.region: two columns matrix specifying polygonal region containing
+ # all data locations. If poly is missing, the unit square is
+ # considered.
+ #
+ # h: a function of the distance between points.
+ # If inhibition=TRUE, h is monotone, increasing, and must tend
+ # to 1 when the distance tends to infinity. 0<= h(d,theta,delta) <= 1.
+ # Else, h is monotone, decreasing, and must tend
+ # to 1 when the distance tends to 0. 0<= h(d,theta,delta) <= 1.
+ #
+ # p: a function among "min", "max", "prod".
+ #
+ # npoints: number of points to simulate.
+ #
+ #
+ # Value:
+ # pattern: list containing the points (x,y) of the simulated point
+ # process.
+ #
+ ##
+ ## E. GABRIEL, 26/03/2007
+ ##
+ ## last modification: 21/10/2008
+ ##
+ ##
+
+ if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
+
+ if (!(is.function(h)))
+ {
+ models <- c("step","gaussian")
+ if (sum(h==models)==0)
+ {
+ message <- paste("this model is not implemented, please choose among: ",paste(models,"",sep=" ",collapse="and "))
+ stop(message)
+ }
+ if (h=="step")
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- rep(1,length(d))
+ if (inhibition==TRUE) res[d<=delta] <- theta
+ else res[d>=delta] <- theta
+ return(res)
+ }
+ }
+ if (h=="gaussian")
+ {
+ hk <- function(d,theta,delta)
+ {
+ if (inhibition==TRUE)
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<=delta) res=c(res,0)
+ if (d[i]>(delta+theta/2)) res=c(res,1)
+ if (d[i]>delta & d[i]<=(delta+theta/2)) res=c(res,exp(-((d[i]-delta-theta/2)^2)/(2*(theta/8)^2)))
+ }
+ }
+ else
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<delta) res=c(res,1)
+ else res=c(res,exp(-((d[i]-delta)^2)/(2*(theta/8)^2)))
+ }
+ }
+ return(res)
+ }
+ }
+ }
+ else
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- h(d,theta,delta)
+ return(res)
+ }
+ }
+
+ pk <- function(d,h,recent,theta,delta)
+ {
+ if (recent=="all")
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(d))
+ {
+ if (p=="min") res <- min(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ }
+ else
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ }
+ else
+ stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+
+ xy <- csr(n=1,poly=s.region)
+ npts <- 1
+ pattern.interm <- cbind(x=xy[1],y=xy[2])
+ if (inhibition==TRUE)
+ {
+ while(npts < npoints)
+ {
+ prob <- runif(1)
+ xy <- csr(n=1,poly=s.region)
+ if (all((sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2)) > delta))
+ umax <- 1
+ else
+ umax <- pk(d=sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2),hk,recent,theta,delta)
+
+ if (prob<umax)
+ {
+ pattern.interm <- rbind(pattern.interm,c(xy[1],xy[2]))
+ npts <- npts+1
+ }
+ }
+ }
+ else
+ {
+ while(npts < npoints)
+ {
+ prob <- runif(1)
+ continue <- FALSE
+ while(continue==FALSE)
+ {
+ xy <- csr(n=1,poly=s.region)
+ if (all(sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2) < delta))
+ umax <- 1
+ else
+ umax <- pk(d=sqrt((xy[1] - pattern.interm[,1])^2 + (xy[2] - pattern.interm[,2])^2),hk,recent,theta,delta)
+
+ if (prob < umax)
+ {
+ pattern.interm <- rbind(pattern.interm,c(xy[1],xy[2]))
+ npts <- npts+1
+ continue <- TRUE
+ }
+ }
+ }
+ }
+ invisible(return(list(pts=pattern.interm,s.region=s.region)))
+}
+
+
Modified: pkg/R/temporal.inhibition.r
===================================================================
--- pkg/R/temporal.inhibition.r 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/R/temporal.inhibition.r 2008-10-31 10:42:57 UTC (rev 21)
@@ -1,235 +1,184 @@
-temporal.inhibition <- function(npoints,h,theta,delta,p,recent="all",t.region,discrete.time=FALSE,replace=FALSE,inhibition=TRUE,...)
- {
- #
- # Simulate an inhibition or a contagious temporal point process in T
- #
- # Requires Splancs.
- #
- # Arguments:
- #
- # t.region: vector containing the minimum and maximum values of
- # the time interval.
- #
- # h: a function of the distance between times and theta.
- # If inhibition=TRUE, h is monotone, increasing, and must tend
- # to 1 when the distance tends to infinity. 0<= h(d,theta) <= 1.
- # Else, h is monotone, decreasing, and must tend
- # to 1 when the distance tends to 0. 0<= h(d,theta) <= 1.
- #
- # p: a function among "min", "max", "prod".
- #
- # replace: logical allowing times repetition.
- #
- # npoints: number of points to simulate.
- #
- #
- # Value:
- # pattern: list containing times t of the simulated process.
- #
- ##
- ## E. GABRIEL, 26/03/2007
- ##
- ## last modification: 21/10/2008
- ##
- ##
-
- if (missing(t.region)) t.region <- c(0,1)
-
- if (inhibition==TRUE)
- {
- if ((max(t.region)-min(t.region))/delta<npoints)
- stop(paste("t.region is too small to fit", npoints, "points", "at minimum time interval", delta))
- }
-
- if (!(is.function(h)))
- {
- models <- c("step","exponential","gaussian")
- if (sum(h==models)==0)
- {
- message <- paste("this model is not implemented, please choose among: ",paste(models,"",sep=" ",collapse="and "))
- stop(message)
- }
- if (h=="step")
- {
- hk <- function(d,theta,delta)
- {
- res <- rep(1,length(d))
- res[d<=delta] <- theta
- return(res)
- }
- }
-
- if (h=="exponential")
- {
- hk <- function(d,theta,delta)
- {
- res <- exp(d*theta)/max(exp(d*theta))
- return(res)
- }
- }
- if (h=="gaussian")
- {
- hk <- function(d,theta,delta)
- {
- res <- exp((d^2)*theta)/max(exp((d^2)*theta))
- return(res)
- }
- }
- }
-
- if (inhibition==FALSE)
- {
- if (!(is.function(h)))
- {
- hh <- hk
- hk <- function(d,theta,delta)
- {
- res <- hh(max(d)-d,theta,max(d)-delta)
- return(res)
- }
- }
- else
- {
- hh <- h
- hk <- function(d,theta,delta)
- {
- res <- h(max(d)-d,theta,max(d)-delta,...)
- return(res)
- }
- }
- }
- if (p=="min")
- {
- pk <- function(d,h,recent,theta,delta,...)
- {
- if (recent=="all")
- res <- min(h(d=d,theta=theta,delta=delta,...))
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(d))
- res <- min(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta,...))
- else
- res <- min(h(d=d,theta=theta,delta=delta,...))
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
-
- if (p=="max")
- {
- pk <- function(d,h,recent,theta,delta,...)
- {
- if (recent=="all")
- res <- max(h(d=d,theta=theta,delta=delta,...))
- else
- {
- if (is.numeric(recent))
- {
- if(recent<=length(d))
- res <- max(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta,...))
- else
- res <- max(h(d=d,theta=theta,delta=delta,...))
- }
- else
- stop("'recent' must be numeric")
- }
- return(res)
- }
- }
-
- if (p=="prod")
- {
- pk <- function(d,h,recent,theta,delta,...)
- {
- if (is.numeric(recent) && recent<=length(d))
- res <- prod(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta,...))
- else
- res <- prod(h(d=d,theta=theta,delta=delta,...))
- return(res)
- }
- }
-
- if (discrete.time==FALSE)
- ti <- runif(1,min=t.region[1],max=t.region[1]+delta)
- else
- ti <- sample(floor(t.region[1]):ceiling(t.region[1]+delta),1)
- times <- ti
- npts <- 1
- if (inhibition==TRUE)
- {
- while(npts < npoints)
- {
- if (discrete.time==FALSE)
- ti <- runif(1,min=t.region[1],max=t.region[2])
- else
- ti <- sample(floor(t.region[1]):ceiling(t.region[2]),1)
-
- prob <- runif(1)
-
- if (all(abs(ti - times) > delta))
- umax <- 1
- else
- {
- if (is.function(h))
- umax <- pk(d=abs(ti - times),h,recent,theta,delta,...)
- else
- {
- h <- hk
- umax <- pk(d=abs(ti - times),h,recent,theta,delta,...)
- }
- }
- if (prob<umax)
- {
- times <- c(times,ti)
- npts <- npts+1
- }
- }
- }
- else
- {
- while(npts < npoints)
- {
- prob <- runif(1)
-
- continue <- FALSE
- while(continue==FALSE)
- {
- if (discrete.time==FALSE)
- ti <- runif(1,min=t.region[1],max=t.region[2])
- else
- ti <- sample(floor(t.region[1]):ceiling(t.region[2]),1)
-
- if (abs(ti - times[npts]) < delta)
- umax <- 1
- else
- {
- if (is.function(h))
- umax <- pk(d=abs(ti - times[npts]),h,recent,theta,delta,...)
- else
- {
- h <- hk
- umax <- pk(d=abs(ti - times[npts]),h,recent,theta,delta,...)
- }
- }
- if (prob < umax)
- {
- times <- c(times,ti)
- npts <- npts+1
- continue <- TRUE
- }
- }
- }
- }
-
- samp <- sample(1:npoints,npoints,replace=replace)
- times <- sort(times[samp])
-
- invisible(return(list(times=times,t.region=t.region)))
-}
-
-
+temporal.inhibition <- function(npoints,h,theta,delta,p,recent="all",t.region,discrete.time=FALSE,replace=FALSE,inhibition=TRUE)
+ {
+ #
+ # Simulate an inhibition or a contagious temporal point process in T
+ #
+ # Requires Splancs.
+ #
+ # Arguments:
+ #
+ # t.region: vector containing the minimum and maximum values of
+ # the time interval.
+ #
+ # h: a function of the distance between times and theta.
+ # If inhibition=TRUE, h is monotone, increasing, and must tend
+ # to 1 when the distance tends to infinity. 0<= h(d,theta,delta) <= 1.
+ # Else, h is monotone, decreasing, and must tend
+ # to 1 when the distance tends to 0. 0<= h(d,theta,delta) <= 1.
+ #
+ # p: a function among "min", "max", "prod".
+ #
+ # replace: logical allowing times repetition.
+ #
+ # npoints: number of points to simulate.
+ #
+ #
+ # Value:
+ # pattern: list containing times t of the simulated process.
+ #
+ ##
+ ## E. GABRIEL, 26/03/2007
+ ##
+ ## last modification: 31/10/2008
+ ##
+ ##
+
+ if (missing(t.region)) t.region <- c(0,1)
+
+ if (!(is.function(h)))
+ {
+ models <- c("step","gaussian")
+ if (sum(h==models)==0)
+ {
+ message <- paste("this model is not implemented, please choose among: ",paste(models,"",sep=" ",collapse="and "))
+ stop(message)
+ }
+ if (h=="step")
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- rep(1,length(d))
+ if (inhibition==TRUE) res[d<=delta] <- theta
+ else res[d>=delta] <- theta
+ return(res)
+ }
+ }
+ if (h=="gaussian")
+ {
+ hk <- function(d,theta,delta)
+ {
+ if (inhibition==TRUE)
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<=delta) res=c(res,0)
+ if (d[i]>(delta+theta/2)) res=c(res,1)
+ if (d[i]>delta & d[i]<=(delta+theta/2)) res=c(res,exp(-((d[i]-delta-theta/2)^2)/(2*(theta/8)^2)))
+ }
+ }
+ else
+ {
+ res=NULL
+ for(i in 1:length(d))
+ {
+ if (d[i]<delta) res=c(res,1)
+ else res=c(res,exp(-((d[i]-delta)^2)/(2*(theta/8)^2)))
+ }
+ }
+ return(res)
+ }
+ }
+ }
+ else
+ {
+ hk <- function(d,theta,delta)
+ {
+ res <- h(d,theta,delta)
+ return(res)
+ }
+ }
+
+ pk <- function(d,h,recent,theta,delta)
+ {
+ if (recent=="all")
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ else
+ {
+ if (is.numeric(recent))
+ {
+ if(recent<=length(d))
+ {
+ if (p=="min") res <- min(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d[(length(d)-recent+1):length(d)],theta=theta,delta=delta))
+ }
+ else
+ {
+ if (p=="min") res <- min(h(d=d,theta=theta,delta=delta))
+ if (p=="max") res <- max(h(d=d,theta=theta,delta=delta))
+ if (p=="prod") res <- prod(h(d=d,theta=theta,delta=delta))
+ }
+ }
+ else stop("'recent' must be numeric")
+ }
+ return(res)
+ }
+
+ if (discrete.time==FALSE)
+ ti <- runif(1,min=t.region[1],max=t.region[1]+delta)
+ else
+ ti <- sample(floor(t.region[1]):ceiling(t.region[1]+delta),1)
+ times <- ti
+ npts <- 1
+ if (inhibition==TRUE)
+ {
+ while(npts < npoints)
+ {
+ if (discrete.time==FALSE)
+ ti <- runif(1,min=t.region[1],max=t.region[2])
+ else
+ ti <- sample(floor(t.region[1]):ceiling(t.region[2]),1)
+
+ prob <- runif(1)
+
+ if (all(abs(ti - times) > delta))
+ umax <- 1
+ else
+ umax <- pk(d=abs(ti - times),hk,recent,theta,delta)
+ if (prob<umax)
+ {
+ times <- c(times,ti)
+ npts <- npts+1
+ }
+ }
+ }
+ else
+ {
+ while(npts < npoints)
+ {
+ prob <- runif(1)
+
+ continue <- FALSE
+ while(continue==FALSE)
+ {
+ if (discrete.time==FALSE)
+ ti <- runif(1,min=t.region[1],max=t.region[2])
+ else
+ ti <- sample(floor(t.region[1]):ceiling(t.region[2]),1)
+
+ if (abs(ti - times[npts]) < delta)
+ umax <- 1
+ else
+ umax <- pk(d=abs(ti - times),hk,recent,theta,delta)
+ if (prob < umax)
+ {
+ times <- c(times,ti)
+ npts <- npts+1
+ continue <- TRUE
+ }
+ }
+ }
+ }
+
+ samp <- sample(1:npoints,npoints,replace=replace)
+ times <- sort(times[samp])
+
+ invisible(return(list(times=times,t.region=t.region)))
+}
+
+
Modified: pkg/man/animation.Rd
===================================================================
--- pkg/man/animation.Rd 2008-10-30 15:33:19 UTC (rev 20)
+++ pkg/man/animation.Rd 2008-10-31 10:42:57 UTC (rev 21)
@@ -7,7 +7,7 @@
\usage{
animation(xyt, s.region, t.region, runtime=1, incident="red", prevalent="pink3", pch=19, cex=0.25, plot.s.region=T, scales=T, border.frac=0.05, add=F)
- }
+}
\arguments{
\item{xyt}{data-matrix containing the (x,y,t)-coordinates. }
Added: pkg/man/rinter.Rd
===================================================================
--- pkg/man/rinter.Rd (rev 0)
+++ pkg/man/rinter.Rd 2008-10-31 10:42:57 UTC (rev 21)
@@ -0,0 +1,176 @@
+\name{rinter}
+\alias{rinter}
+\title{Generate interaction point patterns}
+
+\description{
+Generate one (or several) realisation(s) of the inhibition or contagious process
+in a region S x T.
+}
+
+\usage{
+ rinter(npoints, s.region, t.region, hs="step", gs="min", thetas=0, deltas,
+ ht="step", gt="min", thetat=1, deltat, recent="all", nsim=1,
+ discrete.time=FALSE, replace=FALSE, inhibition=TRUE)
+}
+
+\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{hs, ht}{function which depends on the distance between points
+ and \code{theta}. Can be chosen among "\code{step}" and "\code{gaussian}"
+ or can refer to a user defined function which only depend on d, theta, and delta
+ (see details).
+ If \code{inhibition}=TRUE, \code{h} is monotone, increasing, and must tend
+ to 1 when the distance tends to infinity. 0 \eqn{\leq}{<=} h(d,theta)
+ \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{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{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)
+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{
+# simple inhibition process
+inh1 = rinter(npoints=200,thetas=0,deltas=0.05,thetat=0,deltat=0.001,inhibition=TRUE)
+stani(inh1$xyt)
+
+# inhibition process using hs and ht defined by the user
+hs = function(d,theta,delta,mus=0.1)
+{
+ res=NULL
+ a=(1-theta)/mus
+ b=theta-a*delta
+ for(i in 1:length(d))
+ {
+ if (d[i]<=delta) res=c(res,theta)
+ if (d[i]>(delta+mus)) res=c(res,1)
+ if (d[i]>delta & d[i]<=(delta+mus)) res=c(res,a*d[i]+b)
+ }
+ return(res)
+}
+ht = function(d,theta,delta,mut=0.3)
+{
+ res=NULL
+ a=(1-theta)/mut
+ b=theta-a*delta
+ for(i in 1:length(d))
+ {
+ if (d[i]<=delta) res=c(res,theta)
+ if (d[i]>(delta+mut)) res=c(res,1)
+ if (d[i]>delta & d[i]<=(delta+mut)) res=c(res,a*d[i]+b)
+ }
+ return(res)
+}
+
+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)
+
+# 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)
+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/rpcp.Rd
===================================================================
--- pkg/man/rpcp.Rd (rev 0)
+++ pkg/man/rpcp.Rd 2008-10-31 10:42:57 UTC (rev 21)
@@ -0,0 +1,82 @@
+\name{rpcp}
+\alias{rpcp}
+\title{Generate Poisson cluster point patterns}
+
+\description{
+Generate one (or several) realisation(s) of the Poisson cluster process in a region
+S x T.
+}
+
+\usage{
+ rpcp(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL,
+ mc=NULL, nsim=1, cluster="uniform", maxrad, infectious=TRUE,
+ edge = "larger.region", ...)
+}
+
+\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{nparents}{number of parents. If NULL, \code{nparents} is from a
+ Poisson distribution with intensity \code{lambda}.}
+ \item{npoints}{number of points to simulate. If NULL (default), the
+ number of points is from a Poisson distribution with mean the double integral
+ of the intensity over \code{s.region} and \code{t.region}.}
+ \item{lambda}{intensity of the parent process. Can be either a numeric
+ value, a function, or a character (see \code{\link{rpp}}).
+ If \code{NULL}, it is constant and equal to \code{nparents} / volume
+ of the domain.}
+ \item{mc}{average number of children per parent. It is used when
+ \code{npoints} is NULL.}
+ \item{nsim}{number of simulations to generate. }
+ \item{cluster}{distribution of children: ``uniform'', ``normal'' and
+ ``exponential'' are currently implemented.
+ Either a single value if the distribution in space and time is the
+ same, or a vector of length 2, giving first the spatial distribution of
+ children and then the temporal distribution.}
+ \item{maxrad}{2-vector giving the maximum spatial and temporal
+ variation of children. If missing, \code{maxrad}=(0.05,0.05).
+ The maximum spatial variation is the maximum distance between parent
+ and children (radius of a circle centred at the parent). For a normal
+ distribution of children, \code{maxrad} corresponds to the 2 * standard
+ deviation of location of children relative to their parent, such that children
+ lies in the 95 \% IC of the normal distribution.
+ The maximum temporal variation is the maximum time
+ separating parent and children.}
+ \item{infectious}{If TRUE, offspring's times are always greater than
+ parent's time).}
+ \item{edge}{specify the edge correction to use "larger.region" or "without".}
+ \item{...}{additional parameters of the intensity of the parent process.}
+}
+
+\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{
+# homogeneous Poisson distribution of parents
+
+data(northcumbria)
+pcp1 <- rpcp(nparents=50, npoints=500, s.region=northcumbria, t.region=c(1,365),cluster=c("normal","exponential"), maxrad=c(5000,5))
+animation(pcp1$xyt, s.region=pcp1$s.region, t.region=pcp1$t.region,runtime=5)
+
+# inhomogeneous Poisson distribution of parents
+
+lbda <- function(x,y,t,a){a*exp(-4*y) * exp(-2*t)}
+pcp2 <- rpcp(nparents=50, npoints=500, cluster="normal", lambda=lbda, a=4000/((1-exp(-4))*(1-exp(-2))))
+stani(pcp2$xyt)
+}
Added: pkg/man/rpp.Rd
===================================================================
--- pkg/man/rpp.Rd (rev 0)
+++ pkg/man/rpp.Rd 2008-10-31 10:42:57 UTC (rev 21)
@@ -0,0 +1,91 @@
+\name{rpp}
+\alias{rpp}
+\title{Generate Poisson point patterns}
+
+\description{
+Generate one (or several) realisation(s) of the (homogeneous or inhomogeneous) Poisson process in a region S x T.
+}
+
+\usage{
+rpp(lambda, s.region, t.region, npoints=NULL, nsim=1, replace=TRUE,
+ discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL,
+ Lambda=NULL, mut=NULL,...)
+}
+
+\arguments{
+ \item{lambda}{Spatio-temporal intensity of the Poisson process.
+ If \code{lambda} is a single positive number, the function generates
+ realisations of a homogeneous Poisson process, whilst if \code{lambda}
+ is a function of the form \eqn{\lambda(x,y,t,\dots)}{lambda(x,y,t,...)}
+ or a character it generates realisations of an inhomogeneous Poisson process.}
+ \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{replace}{logical allowing times repeat.}
+ \item{npoints}{number of points to simulate. If \code{NULL}, the
+ 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{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. }
+ \item{lmax}{upper bound for the value of \eqn{\lambda(x,y,t)}{lambda(x,y,t)}, if
+ \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.}
+}
+
+\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{t.index}{vector of times index.}
+\item{Lambda}{nx x ny x nt array of the intensity surface at each time.}
+\item{s.region, t.region, lambda}{parameters passed in argument.}
+}
+
+\author{
+Edith Gabriel <edith.gabriel at univ-avignon.fr> and Peter J Diggle.
+}
+
+\seealso{
+ \code{\link{animation}} and \code{\link{stani}} for plotting space-time point patterns.
+ }
+
+\examples{
+# Homogeneous Poisson process
+# ---------------------------
+hpp1 <- rpp(lambda=200,replace=FALSE)
+stani(hpp1$xyt)
+
+# fixed number of points, discrete time, with time repeat.
+data(northcumbria)
+hpp2 = rpp(npoints=500, s.region=northcumbria, t.region=c(1,1000), discrete.time=TRUE)
+polymap(northcumbria)
+animation(hpp2$xyt, s.region=hpp2$s.region, t.region=hpp2$t.region, runtime=10, add=TRUE)
+
+# Inhomogeneous Poisson process
+# -----------------------------
+
+# intensity defined by a function
+lbda1 = function(x,y,t,a){a*exp(-4*y) * exp(-2*t)}
+ipp1 = rpp(lambda=lbda1, npoints=400, a=3200/((1-exp(-4))*(1-exp(-2))))
+stani(ipp1$xyt)
+
+# intensity defined by a 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=100, ny=100)
+Lt = dim(fmd)[1]*density(fmd[,3], n=200)$y
+ipp2 = rpp(lambda="m", Lambda=Ls$z, mut=Lt, s.region=northcumbria,
+t.region=c(1,200), discrete.time=TRUE)
+image(Ls$x, Ls$y, Ls$z, col=grey((1000:1)/1000)); polygon(northcumbria)
+animation(ipp2$xyt, add=TRUE, cex=0.7, runtime=15)
+}
\ No newline at end of file
More information about the Stpp-commits
mailing list