[Stpp-commits] r34 - in pkg: R doc man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Aug 31 11:46:17 CEST 2010
Author: gabriele
Date: 2010-08-31 11:46:16 +0200 (Tue, 31 Aug 2010)
New Revision: 34
Modified:
pkg/R/covst.r
pkg/R/plot.stpp.R
pkg/R/rpcp.r
pkg/R/rpp.r
pkg/doc/doc.tex
pkg/doc/source.exemples.R
pkg/man/plot.stpp.Rd
pkg/man/rinfec.Rd
pkg/man/rinter.Rd
pkg/man/rlgcp.Rd
pkg/man/rpcp.Rd
pkg/man/rpp.Rd
pkg/src/covst.f
pkg/src/stpp.dll
Log:
Modified: pkg/R/covst.r
===================================================================
--- pkg/R/covst.r 2010-02-13 14:40:00 UTC (rev 33)
+++ pkg/R/covst.r 2010-08-31 09:46:16 UTC (rev 34)
@@ -36,23 +36,29 @@
if (model=="stable")
{
mods <- 2
- if ((param[1] >2) || (param[1]<0)) stop("parameter of stable model must lie in [0,2]")
+ if ((param[1] >2) || (param[1]<0)) stop("Stable model parameter must lie in [0,2]")
modt <- 2
- if ((param[2] >2) || (param[2]<0)) stop("parameter of stable model must lie in [0,2]")
+ if ((param[2] >2) || (param[2]<0)) stop("Stable model parameter must lie in [0,2]")
}
if (model=="cauchy")
{
mods <- 3
- if (param[1]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ if (param[1]<=0) stop("Cauchy model parameter must be strictly positive")
modt <- 3
- if (param[2]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ if (param[2]<=0) stop("Cauchy model parameter must be strictly positive")
}
if (model=="wave")
{
mods <- 4
modt <- 4
}
- if (model=="matern") stop("You must specify another model for the temporal covariance")
+ if (model=="matern")
+ {
+ mods <- 7
+ if (param[2]<=0 | param[1]<=0) stop("Matérn model parameters must be strictly positive")
+ modt <- 7
+ if (param[3]<=0 | param[4]<=0) stop("Matérn model parameters must be strictly positive")
+ }
}
if (length(model)==2)
{
@@ -67,29 +73,37 @@
if (model[1]=="stable")
{
mods <- 2
- if ((param[1] >2) || (param[1]<0)) stop("parameter of stable model must lie in [0,2]")
+ if ((param[1] >2) || (param[1]<0)) stop("Stable model parameter must lie in [0,2]")
}
if (model[2]=="stable")
{
modt <- 2
- if ((param[2] >2) || (param[2]<0)) stop("parameter of stable model must lie in [0,2]")
+ if ((param[2] >2) || (param[2]<0)) stop("Stable model parameter must lie in [0,2]")
}
if (model[1]=="cauchy")
{
mods <- 3
- if (param[1]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ if (param[1]<=0) stop("Cauchy model parmaeter must be strictly positive")
}
if (model[2]=="cauchy")
{
modt <- 3
- if (param[2]<=0) stop("parameter of cauchy model must be strictly greater than 0")
+ if (param[2]<=0) stop("Cauchy model parameter must be strictly positive")
}
if (model[1]=="wave")
mods <- 4
if (model[2]=="wave")
modt <- 4
if (model[1]=="matern")
+ {
mods <- 7
+ if (param[2]<=0 | param[1]<=0) stop("Matérn model parameters must be strictly positive")
+ }
+ if (model[2]=="matern")
+ {
+ modt <- 7
+ if (param[3]<=0 | param[4]<=0) stop("Matérn model parameters must be strictly positive")
+ }
}
}
if (!(isTRUE(separable)))
@@ -118,6 +132,18 @@
return(model=c(mods,modt,mod))
}
+matern = function (d, scale = 1, alpha = 1, nu = 0.5)
+{
+ if (any(d < 0))
+ stop("distance argument must be nonnegative")
+ d <- d * alpha
+ d[d == 0] <- 1e-10
+ k <- 1/((2^(nu - 1)) * gamma(nu))
+ res <- scale * k * (d^nu) * besselK(d, nu)
+ return(res)
+}
+
+
covst <- function(dist,times,separable=TRUE,model,param=c(1,1,1,1,1,2),sigma2=1,scale=c(1,1),plot=TRUE,nlevels=10)
{
@@ -141,7 +167,13 @@
as.double(param),
as.double(sigma2),
as.double(scale))[[1]]
+ if (model[1]==7) mods=matern(dist,scale=scale[1],alpha=param[2],nu=param[1])
+ if (model[2]==7) modt=matern(times,scale=scale[2],alpha=param[4],nu=param[3])
+ if (model[1]==7 & model[2]==7) gs=mods%*%t(modt)
+ if (model[1]==7 & model[2]!=7) gs=mods*gs
+ if (model[2]==7 & model[1]!=7) gs=gs*modt
+
if (plot==TRUE)
{
# image(dist,times,gs,col=grey(((10*max(length(times),length(dist))):1)/(10*max(length(times),length(dist)))),xlab="h",ylab="t",cex.axis=1.5,cex.lab=2,font=2)
Modified: pkg/R/plot.stpp.R
===================================================================
--- pkg/R/plot.stpp.R 2010-02-13 14:40:00 UTC (rev 33)
+++ pkg/R/plot.stpp.R 2010-08-31 09:46:16 UTC (rev 34)
@@ -1,8 +1,10 @@
-plot.stpp <- function(x, s.region=NULL, t.region=NULL, ...)
+plot.stpp <- function(x, s.region=NULL, t.region=NULL, mark=FALSE, mark.cexmin=0.4, mark.cexmax=1.2, mark.col=1, ...)
{
if (inherits(x,"stpp")==TRUE)
{
+ if (mark==FALSE)
+ {
par(mfrow=c(1,2),pty="s")
if (is.null(s.region))
plot(x[,1:2],main="xy-locations",...)
@@ -13,7 +15,43 @@
title("xy-locations")
}
plot(x[,3],cumsum(x[,3]),type="l",xlab="t",ylab="",main="cumulative number",las=1,xlim=t.region)
- }
+ }
+ if (mark==TRUE)
+ {
+ l=dim(x)[1]
+ CEX=seq(mark.cexmin,mark.cexmax,length=l)
+ if (mark.col==0)
+ {
+ par(mfrow=c(1,1),pty="s")
+ if (is.null(s.region))
+ plot(x[,1:2],cex=CEX,...)
+ else
+ {
+ polymap(s.region,xlab="x",ylab="y")
+ points(x[,1:2],cex=CEX,...)
+ }
+ }
+ else
+ {
+ if (mark.col=="black" | mark.col==1)
+ COL=grey((l:1)/l)
+ if (mark.col=="red" | mark.col==2)
+ COL=rgb(l:0, 0, 0, max = l)
+ if (mark.col=="green" | mark.col==3)
+ COL=rgb(0, l:0, 0, max = l)
+ if (mark.col=="blue" | mark.col==4)
+ COL=rgb(0, 0, l:0, max = l)
+ par(mfrow=c(1,1),pty="s")
+ if (is.null(s.region))
+ plot(x[,1:2],col=COL,cex=CEX,...)
+ else
+ {
+ polymap(s.region,xlab="x",ylab="y")
+ points(x[,1:2],col=COL,cex=CEX,...)
+ }
+ }
+ }
+ }
}
getS3method("plot", "stpp", optional = FALSE)
Modified: pkg/R/rpcp.r
===================================================================
--- pkg/R/rpcp.r 2010-02-13 14:40:00 UTC (rev 33)
+++ pkg/R/rpcp.r 2010-08-31 09:46:16 UTC (rev 34)
@@ -1,18 +1,21 @@
-pcp.larger.region <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infecD=TRUE, maxradlarger, ...)
+itexp <- function(u, m, t) { -log(1-u*(1-exp(-t*m)))/m }
+rtexp <- function(n, m, t) { itexp(runif(n), m, t) }
+
+pcp.larger.region <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", dispersion, infecD=TRUE, larger.region, tronc=1, ...)
{
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(dispersion)) dispersion <- c(0.05,0.05)
+ if (length(dispersion)==1) dispersion=rep(dispersion,2)
+ dispersions <- dispersion[1]
+ dispersiont <- dispersion[2]
- if (missing(maxradlarger)) maxradlarger <- maxrad
- maxradlargers <- maxradlarger[1]
- maxradlargert <- maxradlarger[2]
+ if (missing(larger.region)) larger.region <- dispersion
+ larger.regions <- larger.region[1]
+ larger.regiont <- larger.region[2]
if (length(cluster)==1)
{
@@ -29,15 +32,15 @@
s.area <- areapl(s.region)
t.area <- t.region[2]-t.region[1]
- if (maxradlargers==0)
+ if (larger.regions==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))
+ s.larger <- rbind(s.region+larger.regions,s.region-larger.regions,cbind(s.region[,1]+larger.regions,s.region[,2]-larger.regions),cbind(s.region[,1]-larger.regions,s.region[,2]+larger.regions))
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)
+ t.larger <- c(t.region[1]-larger.regiont,t.region[2]+larger.regiont)
if (t.larger[1]<0) t.larger[1] <- 0
@@ -89,42 +92,49 @@
if (s.distr=="uniform")
{
- xp <- xpar+runif(nc,min=-maxrads,max=maxrads)
- yp <- ypar+runif(nc,min=-maxrads,max=maxrads)
+ xp <- xpar+runif(nc,min=-dispersions,max=dispersions)
+ yp <- ypar+runif(nc,min=-dispersions,max=dispersions)
}
if (s.distr=="normal")
{
- xp <- xpar+rnorm(nc,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(nc,mean=0,sd=maxrads/2)
+ xp <- xpar+rnorm(nc,mean=0,sd=dispersions/2)
+ yp <- ypar+rnorm(nc,mean=0,sd=dispersions/2)
}
if (s.distr=="exponential")
{
- xp <- xpar+rexp(nc,rate=1/maxrads)
- yp <- ypar+rexp(nc,rate=1/maxrads)
+ xp <- xpar+rexp(nc,rate=1/dispersions)
+ yp <- ypar+rexp(nc,rate=1/dispersions)
}
if (t.distr=="uniform")
{
if (infecD==TRUE)
- zp <- zpar+runif(nc,min=-maxradt,max=maxradt)
+ zp <- zpar+runif(nc,min=-dispersiont,max=dispersiont)
else
- zp <- zpar+runif(nc,min=0,max=maxradt)
+ zp <- zpar+runif(nc,min=0,max=dispersiont)
}
if (t.distr=="normal")
{
if (infecD==TRUE)
- zp <- zpar+abs(rnorm(nc,mean=0,sd=maxradt/2))
+ zp <- zpar+abs(rnorm(nc,mean=0,sd=dispersiont/2))
else
- zp <- zpar+rnorm(nc,mean=0,sd=maxradt/2)
+ zp <- zpar+rnorm(nc,mean=0,sd=dispersiont/2)
}
if (t.distr=="exponential")
{
if (infecD==TRUE)
- zp <- zpar+rexp(nc,rate=1/maxradt)
+ zp <- zpar+rexp(nc,rate=1/dispersiont)
else
- zp <- zpar+sample(c(-1,1),1)*rexp(nc,rate=1/maxradt)
+ zp <- zpar+sample(c(-1,1),1)*rexp(nc,rate=1/dispersiont)
}
+ if (t.distr=="trunc.exponential")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+rtexp(nc,m=1/dispersiont,t=tronc)
+ else
+ zp <- zpar+sample(c(-1,1),1)*rtexp(nc,m=1/dispersiont,t=tronc)
+ }
- 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))
+ mask <- ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/dispersions^2) + (((ypar-yp)^2)/dispersions^2)) < 1))
nchild2[ipar] <- sum(mask)
pattern.interm <- rbind(pattern.interm,cbind(xp[mask],yp[mask],zp[mask],rep(ipar,sum(mask))))
}
@@ -145,42 +155,49 @@
if (s.distr=="uniform")
{
- xp <- xpar+runif(1,min=-maxrads,max=maxrads)
- yp <- ypar+runif(1,min=-maxrads,max=maxrads)
+ xp <- xpar+runif(1,min=-dispersions,max=dispersions)
+ yp <- ypar+runif(1,min=-dispersions,max=dispersions)
}
if (s.distr=="normal")
{
- xp <- xpar+rnorm(1,mean=0,sd=maxrads/2)
- yp <- ypar+rnorm(1,mean=0,sd=maxrads/2)
+ xp <- xpar+rnorm(1,mean=0,sd=dispersions/2)
+ yp <- ypar+rnorm(1,mean=0,sd=dispersions/2)
}
if (s.distr=="exponential")
{
- xp <- xpar+rexp(1,rate=1/maxrads)
- yp <- ypar+rexp(1,rate=1/maxrads)
+ xp <- xpar+rexp(1,rate=1/dispersions)
+ yp <- ypar+rexp(1,rate=1/dispersions)
}
if (t.distr=="uniform")
{
if (infecD==TRUE)
- zp <- zpar+runif(1,min=-maxradt,max=maxradt)
+ zp <- zpar+runif(1,min=-dispersiont,max=dispersiont)
else
- zp <- zpar+runif(1,min=0,max=maxradt)
+ zp <- zpar+runif(1,min=0,max=dispersiont)
}
if (t.distr=="normal")
{
if (infecD==TRUE)
- zp <- zpar+abs(rnorm(1,mean=0,sd=maxradt/2))
+ zp <- zpar+abs(rnorm(1,mean=0,sd=dispersiont/2))
else
- zp <- zpar+rnorm(1,mean=0,sd=maxradt/2)
+ zp <- zpar+rnorm(1,mean=0,sd=dispersiont/2)
}
if (t.distr=="exponential")
{
if (infecD==TRUE)
- zp <- zpar+rexp(1,rate=1/maxradt)
+ zp <- zpar+rexp(1,rate=1/dispersiont)
else
- zp <- zpar+sample(c(-1,1),1)*rexp(1,rate=1/maxradt)
+ zp <- zpar+sample(c(-1,1),1)*rexp(1,rate=1/dispersiont)
}
+ if (t.distr=="trunc.exponential")
+ {
+ if (infecD==TRUE)
+ zp <- zpar+rtexp(1,m=1/dispersiont,t=tronc)
+ else
+ zp <- zpar+sample(c(-1,1),1)*rtexp(1,m=1/dispersiont,t=tronc)
+ }
- 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))
+ if ((inout(as.points(x=xp,y=yp),poly=s.region)==T) & (zp > t.region[1] & zp < t.region[2]) & (sqrt( (((xpar-xp)^2)/dispersions^2) + (((ypar-yp)^2)/dispersions^2)) < 1))
{
pattern.interm <- rbind(pattern.interm,cbind(xp,yp,zp,ipar))
ipt <- ipt+1
@@ -196,7 +213,7 @@
-rpcp <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", maxrad, infectious=TRUE, edge = "larger.region", ...)
+rpcp <- function(s.region, t.region, nparents=NULL, npoints=NULL, lambda=NULL, mc=NULL, nsim=1, cluster="uniform", dispersion, infectious=TRUE, edge = "larger.region", larger.region=larger.region, tronc=1,...)
{
#
# Simulate a space-time Poisson cluster process in a region D x T.
@@ -237,12 +254,12 @@
# first the spatial distribution of children and then
# the temporal distribution.
#
- # maxrad: vector of length 2 giving the maximum spatial and temporal
+ # dispersion: vector of length 2 giving the maximum spatial and temporal
# variation of the offspring.
- # maxrads = maximum distance between parent and child (radius of
+ # dispersions = maximum distance between parent and child (radius of
# a circle centred at the parent).
- # maxradt = maximum time separiting parent and child.
- # For a normal distribution of children, maxrad corresponds to
+ # dispersiont = maximum time separiting parent and child.
+ # For a normal distribution of children, dispersion 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.
@@ -266,9 +283,9 @@
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]
+ if (missing(dispersion)) dispersion <- c(0.05,0.05)
+ dispersions <- dispersion[1]
+ dispersiont <- dispersion[2]
if (length(cluster)==1)
{
@@ -291,10 +308,10 @@
while(ni<=nsim)
{
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=infectious, ...)$pts
+ pattern.interm <- pcp.larger.region(s.region=s.region, t.region=t.region, nparents=nparents, npoints=npoints, lambda=lambda, mc=mc, cluster=cluster, dispersion=dispersion, infecD=infectious, larger.region=larger.region,tronc=tronc, ...)$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=infectious, 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, dispersion=dispersion, infecD=infectious, larger.region=c(0,0), tronc=tronc, ...)$pts
if (nsim==1)
pattern <- as.3dpoints(pattern.interm)
Modified: pkg/R/rpp.r
===================================================================
--- pkg/R/rpp.r 2010-02-13 14:40:00 UTC (rev 33)
+++ pkg/R/rpp.r 2010-08-31 09:46:16 UTC (rev 34)
@@ -156,7 +156,7 @@
}
-ripp <- function(lambda, s.region, t.region, npoints=NULL, replace=T, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, Lambda=NULL, mut=NULL, ...)
+ripp <- function(lambda, s.region, t.region, npoints=NULL, replace=TRUE, discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, Lambda=NULL, ...)
{
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)
@@ -197,19 +197,15 @@
else
t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=(t.area/(nt-1)))
- if (is.null(Lambda) & is.null(mut))
+ if (is.null(Lambda))
{
- Lambda <- array(NaN,dim=c(nx,ny,nt))
- mut <- rep(0,nt)
- # maxlambda <- rep(0,nt)
+ Lambda <- array(NaN,dim=c(nx,ny,nt))
for(it in 1:nt)
{
L <- lambda(as.vector(s.grid$X),as.vector(s.grid$Y),t.grid$times[it],...)
M <- matrix(L,ncol=ny,nrow=nx,byrow=T)
M[!(s.grid$mask)] <- NaN
Lambda[,,it] <- M
- # maxlambda[it] <- max(Lambda[,,it],na.rm=T)
- mut[it] <- sum(Lambda[,,it],na.rm=T)
}
}
@@ -242,7 +238,7 @@
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=T))
+ samp <- sample(1:nt,npts,replace=replace)
times <- times.init[samp]
prob <- lambda(x,y,times,...)/lambdamax
u <- runif(npts)
@@ -270,8 +266,8 @@
if(dim(xy)[2]==1){wx <- xy[1]; wy <- xy[2]}
else{wx <- xy[,1]; wy <- xy[,2]}
if(replace==F)
- { wsamp <- sample(samp.remain,npoints-neffec,replace=replace,prob=mut[samp.remain]/max(mut[samp.remain],na.rm=T)) }
- else{ wsamp <- sample(1:nt,npoints-neffec,replace=replace,prob=mut/max(mut,na.rm=T)) }
+ { wsamp <- sample(samp.remain,npoints-neffec,replace=replace) }
+ else{ wsamp <- sample(1:nt,npoints-neffec,replace=replace) }
wtimes <- times.init[wsamp]
# lambdamax <- maxlambda[wsamp]
prob <- lambda(wx,wy,wtimes,...)/lambdamax
@@ -293,11 +289,11 @@
if (is.character(lambda))
{
- if (is.null(Lambda) & is.null(mut))
- stop("Lambda and mut must be specified")
+ if (is.null(Lambda))
+ stop("Lambda must be specified")
nx <- dim(Lambda)[1]
ny <- dim(Lambda)[2]
- nt <- length(mut)
+ nt <- dim(Lambda)[3]
s.grid <- make.grid(nx,ny,s.region)
s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
@@ -322,15 +318,15 @@
if (is.null(npoints))
{
- en <- sum(Lambda,na.rm=T)*s.grid$xinc*s.grid$yinc
+ en <- sum(Lambda,na.rm=T)*s.grid$xinc*s.grid$yinc*t.grid$tinc
npoints <- round(rpois(n=1,lambda=en),0)
}
if (is.null(lambdamax))
- lambdamax <- max(Lambda,na.rm=T)*max(mut)/npoints
+ lambdamax <- max(Lambda,na.rm=T)
# npts <- round(lambdamax/(s.area*t.area),0)
npts <- npoints
if (npts==0) stop("there is no data to thin")
-
+
if ((replace==F) && (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
if (discrete.time==T)
{
@@ -339,8 +335,8 @@
}
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=T))
+
+ samp <- sample(1:nt,npts,replace=replace)
times <- times.init[samp]
retain.eq.F <- F
@@ -349,14 +345,14 @@
xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
x <- xy[,1]
y <- xy[,2]
-
+
prob <- NULL
for(nx in 1:length(x))
{
nix <- iplace(X=s.grid$x,x=x[nx],xinc=s.grid$xinc)
niy <- iplace(X=s.grid$y,x=y[nx],xinc=s.grid$yinc)
nit <- iplace(X=t.grid$times,x=times[nx],xinc=t.grid$tinc)
- prob <- c(prob,(Lambda[nix,niy]*mut[nit]/npoints)/lambdamax)
+ prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
M <- which(is.na(prob))
@@ -376,6 +372,18 @@
}
# if (sum(retain==F)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
+
+ if (sum(retain==F)==length(retain))
+ {
+ lambdas <- matrix(0,nrow=nx,ncol=ny)
+ for(ix in 1:nx){for(iy in 1:ny){
+ lambdas[ix,iy] <- median(Lambda[ix,iy,],na.rm=T)}}
+ lambdamax <- max(lambdas,na.rm=T)
+ prob <- lambda(x,y,times,...)/lambdamax
+ retain <- u <= prob
+ if (sum(retain==F)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
+ }
+
x <- x[retain]
y <- y[retain]
samp <- samp[retain]
@@ -399,8 +407,8 @@
if(dim(xy)[2]==1){wx <- xy[1]; wy <- xy[2]}
else{wx <- xy[,1]; wy <- xy[,2]}
if(replace==F)
- { wsamp <- sample(samp.remain,npoints-neffec,replace=replace,prob=mut[samp.remain]/max(mut[samp.remain],na.rm=T)) }
- else{ wsamp <- sample(1:nt,npoints-neffec,replace=replace,prob=mut/max(mut,na.rm=T)) }
+ { wsamp <- sample(samp.remain,npoints-neffec,replace=replace) }
+ else{ wsamp <- sample(1:nt,npoints-neffec,replace=replace) }
wtimes <- times.init[wsamp]
# lambdamax <- maxlambda[wsamp]
@@ -410,7 +418,7 @@
nix <- iplace(X=s.grid$x,x=wx[nx],xinc=s.grid$xinc)
niy <- iplace(X=s.grid$y,x=wy[nx],xinc=s.grid$yinc)
nit <- iplace(X=t.grid$times,x=wtimes[nx],xinc=t.grid$tinc)
- prob <- c(prob,(Lambda[nix,niy]*mut[nit]/npoints)/lambdamax)
+ prob <- c(prob,(Lambda[nix,niy,nit])/lambdamax)
}
M <- which(is.na(prob))
if (length(M)!=0)
@@ -443,7 +451,7 @@
}
-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, ...)
+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, ...)
{
#
# Simulate a space-time Poisson process in a region D x T.
@@ -483,15 +491,14 @@
# 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.
+ # Lambda: array of spatial intensity if 'lambda' is a character.
#
- # mut: vector of temporal intensity if 'lambda' is a character.
#
#
# Value:
# 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.
+ # Lambda: an array of the intensity surface
#
# NB: the probability that an event occurs at a location s and a time t
# is:
@@ -567,21 +574,17 @@
t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=(t.area/(nt-1)))
Lambda <- array(NaN,dim=c(nx,ny,nt))
- mut <- rep(0,nt)
- # maxlambda <- rep(0,nt)
for(it in 1:nt)
{
L <- lambda(as.vector(s.grid$X),as.vector(s.grid$Y),t.grid$times[it],...)
M <- matrix(L,ncol=ny,nrow=nx,byrow=T)
M[!(s.grid$mask)] <- NaN
Lambda[,,it] <- M
- # maxlambda[it] <- max(Lambda[,,it],na.rm=T)
- mut[it] <- sum(Lambda[,,it],na.rm=T)
}
while(ni<=nsim)
{
- ipp <- ripp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time, nx=nx, ny=ny, nt=nt, lmax=lmax, Lambda=Lambda, mut=mut, ...)
+ ipp <- ripp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time, nx=nx, ny=ny, nt=nt, lmax=lmax, Lambda=Lambda, ...)
if (nsim==1)
{
@@ -601,7 +604,7 @@
{
while(ni<=nsim)
{
- ipp <- ripp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time, nx=nx, ny=ny, nt=nt, lmax=lmax, Lambda=Lambda, mut=mut, ...)
+ ipp <- ripp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time, nx=nx, ny=ny, nt=nt, lmax=lmax, Lambda=Lambda, ...)
if (nsim==1)
{
Modified: pkg/doc/doc.tex
===================================================================
--- pkg/doc/doc.tex 2010-02-13 14:40:00 UTC (rev 33)
+++ pkg/doc/doc.tex 2010-08-31 09:46:16 UTC (rev 34)
@@ -311,7 +311,7 @@
disease of farm livestock. The dataset \verb#fmd# contains a three-column matrix
of spatial locations and reported days (counted from 1st February) of FMD in north
Cumbria. Figure~\ref{fig:fmd} shows two static displays of the data from the county
-of north Cumbria.
+of north Cumbria.
The left-hand panel shows the locations of reported cases.
The right-hand panel shows the cumulative distribution of the times,
illustrating the characteristic S-shape of an epidemic process.
@@ -804,7 +804,7 @@
$g$ functions, "step" for $h$ functions and 0 for $\theta$ parameters.
The commands
\begin{verbatim}
-> inh1 = rinter(npoints=200, thetas=0, deltas=0.05, thetat=0, deltat=0.001,
+> inh1 = rinter(npoints=200, thetas=0, deltas=0.05, thetat=0, deltat=0.001,
inhibition=TRUE)
> stani(inh1$xyt)
\end{verbatim}
@@ -1316,283 +1316,3 @@
\end{document}
-
-\subsection{{\tt covst}}
-
-\subsubsection*{Description}
-
-Computes the covariances for pairs variables, given the separation distance
-and time of their locations in $S \times T$.
-
-\subsubsection*{Usage}
-
-\begin{verbatim}
- covst(dist, times, separable=TRUE, model, param=c(1,1,1,1,1,2), sigma2=1,
- scale=c(1,1), plot=TRUE, nlevels=10)
-\end{verbatim}
-
-\subsubsection*{Arguments}
-
-\begin{tabular}[c]{rl}
- \verb#dist#: & vector of distances h at which the covariance is computed, \\
- \verb#times#: & vector of times t at which the covariance is computed, \\
- \verb#separable#: & Logical. If \verb#TRUE#, the covariance function is
- separable. \\
- \verb#model#: & Vector of length 1 or 2 specifying the model(s) of
- covariance. \\
- & If \verb#separable=TRUE# and \verb#model# is of length 2,
- then the elements of \verb#model# define \\
- & the spatial and temporal covariances respectively.\\
- & If \verb#separable=TRUE# and \verb#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 \verb#separable=FALSE#, \verb#model# must be of length 1 and is either
- ``gneiting'' or \\
- & ``cesare'' . \\
- \verb#param#: & Vector of parameters of the covariance function (see
- Details) \\
- \verb#sigma2#: & Variance parameter.\\
- \verb#scale#: & Scale parameter. \\
- \verb#plot#: & Logical. If \verb#TRUE#, provide an image in grey scale
- and contour lines are added \\
- & to the plot. \\
- \verb#nlevels#: & Number of contour levels desired, if \verb#plot=TRUE#.
-\end{tabular}
-
-\subsubsection*{Details}
-
-We implemented stationary, isotropic spatio-temporal covariance functions.
-Covariance functions return the value of the covariance $\cov(h,t)$ which
-can be written as a product of a variance parameter $\sigma^2$ times a
-positive definite function $c(h,t)$: $\cov(h,t) = \sigma^2 c(h,t)$.
-
-\medskip
-
-\noindent
-{\em Separable covariance functions}
-\begin{equation*}
- c(h,t) = c_s(\| h \|) \, c_t(|t|) , h \in S, t \in T,
-\end{equation*}
-Models for purely spatial and purely temporal covariance functions:
-\begin{itemize}
-\item Mat{\'e}rn: $c(r) = \frac{(r/\beta)^\alpha}{2^{\alpha-1}\Gamma(\alpha)}
- {\cal K}_{\alpha}(r/\beta)$, $\alpha > 0$ and $\beta > 0$.
-
-$\cal{K}_{\alpha}$ is the modified Bessel function of second kind:
-\begin{equation*}
- {\mathcal{K}}_{\alpha}(x) = \frac{\pi}{2} \frac{I_{-\alpha}(x) -
- I_{\alpha}(x)}{\sin(\pi \alpha)}, \text{ with }
- 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}.
-\end{equation*}
-\item Exponential: $c(r) = \exp(-r/\beta)$, \ $\beta>0$,
-\item Stable: $c(r) = \exp(-(r/\beta)^\alpha)$, \ $\alpha \in [0,2]$,
- $\beta>0$,
-\item Cauchy: $c(r) = (1+(r/\beta)^2)^{-\alpha}$, \ $\alpha >0$,
-\item Wave: $c(r) = \frac{\beta}{r}\sin \left(\frac{r}{\beta} \right)$
- \ $r>0$, $\beta>0$ and $c(0)=1$,
-\end{itemize}
-Parameters $\alpha_1$ and $\alpha_2$ correspond to the parameters of
-the spatial and temporal covariance respectively and parameters
-$\beta_1$ and $\beta_2$ correspond to the spatial and temporal scales.
-
-\medskip
-
-\noindent
-{\em Non-separable covariance functions}
-
-\noindent
-The spatio-temporal covariance function can be:
-\begin{itemize}
-\item gneiting: $c(h,t) = \psi (t/\beta_2)^{-\alpha_6} \phi \left(\frac{h/
- \beta_1}{\psi (t/\beta_2)} \right),$ $\beta_1, \beta_2 >0$,
- \begin{itemize}
- \item If $\alpha_2=1$, $\phi(r)$ is the Stable model.
- \item if $\alpha_2=2$, $\phi(r)$ is the Cauchy model.
- \item if $\alpha_2=3$, $\phi(r)$ is the Mat{\'e}rn model.
- \item If $\alpha_5=1$, $\psi^2(r) = (r^{\alpha_3}+1)^{\alpha_4}$,
- \item If $\alpha_5=2$, $\psi^2(r) = (\alpha_4^{-1} r^{\alpha_3}
- +1)/(r^{\alpha_3}+1)$,
- \item If $\alpha_5=3$, $\psi^2(r) = -\log(r^{\alpha_3}+1/{\alpha_4})/
- \log {\alpha_4}$,
- \end{itemize}
- The parameter $\alpha_1$ is the respective parameter for the model of
- $\phi(.)$, $\alpha_3 \in (0,2]$, $\alpha_4 \in (0,1]$ and $\alpha_6 \geq 2$.
-\item cesare: $c(h,t) = \left( 1 + (h/\beta_1)^{\alpha_1} +
- (t/\beta_2)^{\alpha_2} \right)^{-\alpha_3}$,
- $\beta_1, \beta_2 >0$, $\alpha_1, \alpha_2 \in [1,2]$
- and $\alpha_3 \geq 3/2$.
-\end{itemize}
-
-
-\subsubsection*{Value}
-
-{\tiny $\bullet$} \verb#gs#: matrix of covariance values.
-
-\subsubsection*{References}
-
-\noindent
-Gneiting T. (2002). Nonseparable, stationary covariance functions
-for space-time data. {\em Journal of the American Statistical Association},
-{\bf 97}, 590--600.
-
-\subsubsection*{Author(s)}
-
-Edith Gabriel, {\tt edith.gabriel at univ-avignon.fr}.
-
-\subsubsection*{Example}
-
-\begin{verbatim}
- # separable
-
- M <- covst(dist=seq(0,1,length=100),times=seq(0,1,length=100),separable=TRUE,
- model=c("cauchy","exponential"),param=c(2,1,1,1,1,2))
-
- M <- covst(dist=seq(0,5,length=100),times=seq(0,20,length=100),separable=TRUE,
- model=c("matern","wave"),param=c(3,1,1,1,1,2),nlevels=5)
-
- # non separable
-
- M <- covst(dist=seq(0,1,length=100),times=seq(0,1,length=100),separable=FALSE,
- model="gneiting",param=c(1,2,1,1/2,2,2))
-
- M <- covst(dist=seq(0,1,length=100),times=seq(0,1,length=100),separable=FALSE,
- model="cesare",param=c(1.5,1.8,2,1/2,1,2))
-\end{verbatim}
-
-
-\clearpage
-
-\subsection{{\tt make.grid}}
-
-\subsubsection*{Description}
-
-Generate a rectangular grid of points in a polygon.
-
-\subsubsection*{Usage}
-
-\begin{verbatim}
- make.grid(nx, ny, poly)
-\end{verbatim}
-
-\subsubsection*{Arguments}
-
-\begin{tabular}[c]{rl}
- \verb#nx, ny# & number of pixels in each row and each column of the
- rectangular grid,\\
- \verb#poly# & two-column matrix containing the $x, y$-coordinates of the
- vertices of the polygon \\
- & boundary. If missing, \verb#poly# is the unit square.
-\end{tabular}
-
-%\subsubsection*{Details}
-%
-%Let us denote ${\cal P}$ the polygon and $(x_p,y_p)$ the points defining
-%${\cal P}$. The {\tt nx}$\times${\tt ny} grid is built according to the
-%algorithm:
-%
-%1. Compute the mesh size of the grid:
-%$$x_{inc}=(\max(x_p)-\min(x_p))/{\tt nx} \ \ ; \ \
-%y_{inc}=(\max(y_p)-\min(y_p))/{\tt ny}.$$
-%
-%2. Define the grid nodes $(x ,y)$ by:
-%
-%\hspace{1cm} a) $x_1=\min(x_p)+x_{inc}/2 ; y_1=\min(y_p)+y_inc/2$
-%
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/stpp -r 34
More information about the Stpp-commits
mailing list