[Stpp-commits] r56 - in pkg: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 15 15:40:25 CEST 2012
Author: gabriele
Date: 2012-08-15 15:40:25 +0200 (Wed, 15 Aug 2012)
New Revision: 56
Added:
pkg/R/PCFhat.R
pkg/R/plotPCF.R
pkg/man/PCFhat.Rd
pkg/man/plotPCF.Rd
pkg/man/stpp.Rd
pkg/src/pcffunction.f
Modified:
pkg/DESCRIPTION
pkg/R/STIKhat.r
pkg/R/plotK.R
pkg/R/rinfec.r
pkg/R/rlgcp.r
pkg/R/rpp.r
pkg/man/STIKhat.Rd
pkg/man/plotK.Rd
pkg/man/rlgcp.Rd
pkg/man/rpcp.Rd
pkg/man/rpp.Rd
pkg/man/sim.stpp.Rd
pkg/src/stikfunction.f
Log:
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-06-04 12:29:36 UTC (rev 55)
+++ pkg/DESCRIPTION 2012-08-15 13:40:25 UTC (rev 56)
@@ -1,11 +1,11 @@
Package: stpp
Type: Package
Title: Space-Time Point Pattern simulation, visualisation and analysis
-Version: 0.3
-Date: 2012-04-21
+Version: 1.0
+Date: 2012-08-08
Author: Edith Gabriel, Peter J Diggle, stan function by Barry Rowlingson
Maintainer: Edith Gabriel <edith.gabriel at univ-avignon.fr>
-Depends: R (>= 2.10), splancs
+Depends: R (>= 2.10), splancs, KernSmooth
Suggests: rpanel, rgl
Description: A package for analysing, simulating and displaying space-time point patterns
License: GPL-3
Added: pkg/R/PCFhat.R
===================================================================
--- pkg/R/PCFhat.R (rev 0)
+++ pkg/R/PCFhat.R 2012-08-15 13:40:25 UTC (rev 56)
@@ -0,0 +1,101 @@
+PCFhat <- function(xyt, s.region, t.region, dist, times, lambda, ks="box", hs, kt="box", ht, correction = TRUE)
+{
+require(splancs)
+require(KernSmooth)
+
+ if (missing(s.region)) s.region <- sbox(xyt[,1:2],xfrac=0.01,yfrac=0.01)
+ if (missing(t.region)) t.region <- range(xyt[,3],na.rm=T)
+
+ pts <- xyt[,1:2]
+ xytimes <- xyt[,3]
+ ptsx <- pts[, 1]
+ ptsy <- pts[, 2]
+ ptst <- xytimes
+ npt <- length(ptsx)
+ ndist <- length(dist)
+ dist <- sort(dist)
+ ntimes <- length(times)
+ times <- sort(times)
+ bsupt <- max(t.region)
+ binft <- min(t.region)
+
+ area <- areapl(s.region)*(bsupt-binft)
+
+ np <- length(s.region[, 1])
+ polyx <- c(s.region[, 1], s.region[1, 1])
+ polyy <- c(s.region[, 2], s.region[1, 2])
+ pcfhat <- array(0, dim = c(ndist,ntimes))
+
+ frac=1
+ if (missing(hs))
+ {
+ d=dist(pts)
+ if (ks=="gaussian") hs = dpik(d,kernel="normal",range.x=c(min(d),max(d)/frac))
+ else hs = dpik(d,kernel=ks,range.x=c(min(d),max(d)/frac))
+ }
+ if (missing(ht))
+ {
+ d=dist(ptst)
+ if (kt=="gaussian") ht = dpik(d,kernel="normal",range.x=c(min(d),max(d)/frac))
+ else ht = dpik(d,kernel=kt,range.x=c(min(d),max(d)/frac))
+ }
+
+ kernel=c(ks=ks,hs=hs,kt=kt,ht=ht)
+
+ if (ks=="box") ks=1
+ else if (ks=="epanech") ks=2
+ else if (ks=="gaussian") ks=3
+ else if (ks=="biweight") ks=4
+
+ if (kt=="box") kt=1
+ else if (kt=="epanech") kt=2
+ else if (kt=="gaussian") kt=3
+ else if (kt=="biweight") kt=4
+
+ if(missing(lambda))
+ {
+ misl <- 1
+ lambda <- rep(1,npt)
+ }
+ else misl <- 0
+ if (length(lambda)==1) lambda <- rep(lambda,npt)
+ if (correction==TRUE){edg <- 1} else edg <- 0
+ storage.mode(pcfhat) <- "double"
+
+ if (area>10) lambda <- lambda*area
+
+ nev <- rep(0,ntimes)
+ klist <- .Fortran("pcffunction", as.double(ptsx),
+ as.double(ptsy), as.double(ptst),
+ as.integer(npt), as.double(polyx),
+ as.double(polyy), as.integer(np),
+ as.double(dist), as.integer(ndist),
+ as.double(times), as.integer(ntimes),
+ as.double(bsupt), as.double(binft),
+ as.double(lambda), as.integer(ks), as.integer(kt),
+ as.integer(edg), as.double(hs),
+ as.double(ht), (pcfhat))
+ pcfhat <- klist[[20]]
+
+ if (misl==1)
+ {
+ if (area>10)
+ pcfhat <- ((area^3)/(npt*(npt-1)))*pcfhat
+ else
+ pcfhat <- (area/(npt*(npt-1)))*pcfhat
+ }
+ else
+ {
+ if (area>10)
+ pcfhat <- pcfhat*area
+ else
+ pcfhat <- pcfhat/area
+ }
+
+ pcfhat <- pcfhat/(4*pi*dist)
+
+ if(dist[1]==0) pcfhat[1,]=NA
+ if(times[1]==0) pcfhat[,1]=NA
+
+ invisible(return(list(pcf=pcfhat,dist=dist,times=times,kernel=kernel)))
+}
Modified: pkg/R/STIKhat.r
===================================================================
--- pkg/R/STIKhat.r 2012-06-04 12:29:36 UTC (rev 55)
+++ pkg/R/STIKhat.r 2012-08-15 13:40:25 UTC (rev 56)
@@ -38,6 +38,7 @@
if (area>10) lambda <- lambda*area
+
nev <- rep(0,ntimes)
klist <- .Fortran("stikfunction", as.double(ptsx),
as.double(ptsy), as.double(ptst),
Modified: pkg/R/plotK.R
===================================================================
--- pkg/R/plotK.R 2012-06-04 12:29:36 UTC (rev 55)
+++ pkg/R/plotK.R 2012-08-15 13:40:25 UTC (rev 56)
@@ -1,5 +1,19 @@
plotK <- function(K,n=15,L=TRUE,persp=FALSE,legend=TRUE,...)
{
+devl=dev.list()
+if(is.null(devl))
+{
+old.par <- par(no.readonly = TRUE)
+on.exit(par(old.par))
+}
+else
+{
+dev.off(devl[length(devl)])
+X11()
+old.par <- par(no.readonly = TRUE)
+on.exit(par(old.par))
+}
+
if (isTRUE(L))
{
k <- K$Khat-K$Ktheo
@@ -23,7 +37,7 @@
if (isTRUE(persp))
{
mask <- matrix(0,ncol=length(K$times),nrow=length(K$dist))
- for(i in 1:length(K$dist)){ for(j in 1:length(K$times)){mask[i,j] <- COL[.iplace(x=k[i,j],X=M,xinc=diff(range(M))/(n-1))]}}
+ for(i in 1:length(K$dist)){ for(j in 1:length(K$times)){mask[i,j] <- COL[findInterval(x=k[i,j],vec=M)]}}
COL <- mask[1:(length(K$dist)-1),1:(length(K$times)-1)]
if(isTRUE(legend))
@@ -33,8 +47,8 @@
persp(x=K$dist, y=K$times, z=k, xlab="u",ylab="v", zlab="",expand=1, col=COL, ...)
title(titl,cex.main=2)
par(fig=c(0.825,1,0,1))
- mini <- .iplace(x=min(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
- maxi <- .iplace(x=max(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
+ mini <- findInterval(x=min(k,na.rm=TRUE),vec=M)
+ maxi <- findInterval(x=max(k,na.rm=TRUE),vec=M)
legend("right",fill=colo(n)[maxi:mini],legend=M[maxi:mini],horiz=F,bty="n")
}
else
@@ -60,8 +74,8 @@
axis(2,at=at[length(at)],labels="v",cex.axis=2)
title(titl,cex.main=2)
par(fig=c(0,1,0.1,1))
- mini <- .iplace(x=min(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
- maxi <- .iplace(x=max(k,na.rm=T),X=M,xinc=diff(range(k))/(n-1))
+ mini <- findInterval(x=min(k,na.rm=TRUE),vec=M)
+ maxi <- findInterval(x=max(k,na.rm=TRUE),vec=M)
legend("right",fill=colo(n)[maxi:mini],legend=M[maxi:mini],horiz=F,bty="n")
}
else
@@ -78,5 +92,6 @@
title(titl,cex.main=2)
}
}
+par(old.par)
}
Added: pkg/R/plotPCF.R
===================================================================
--- pkg/R/plotPCF.R (rev 0)
+++ pkg/R/plotPCF.R 2012-08-15 13:40:25 UTC (rev 56)
@@ -0,0 +1,88 @@
+plotPCF <- function(PCF,n=15,persp=FALSE,legend=TRUE,...)
+{
+devl=dev.list()
+if(is.null(devl))
+{
+old.par <- par(no.readonly = TRUE)
+on.exit(par(old.par))
+}
+else
+{
+dev.off(devl[length(devl)])
+X11()
+old.par <- par(no.readonly = TRUE)
+on.exit(par(old.par))
+}
+
+ k=PCF$pcf
+ K=PCF
+
+ titl <- expression(hat(g)* group("(",list(u,v),")") )
+
+
+ colo <- colorRampPalette(c("red", "white", "blue"))
+ M <- max(abs(range(k)))
+ M <- pretty(c(-M,M),n=n)
+ n <- length(M)
+ COL <- colo(n)
+ if (isTRUE(persp))
+ {
+ mask <- matrix(0,ncol=length(K$times),nrow=length(K$dist))
+ for(i in 1:length(K$dist)){ for(j in 1:length(K$times)){mask[i,j] <- COL[findInterval(x=k[i,j],vec=M)]}}
+ COL <- mask[1:(length(K$dist)-1),1:(length(K$times)-1)]
+
+ if(isTRUE(legend))
+ {
+ par(cex.lab=2,cex.axis=1.5,font=2,lwd=1,mar=c(0,0,3,0))
+ par(fig=c(0,0.825,0,1))
+ persp(x=K$dist, y=K$times, z=k, xlab="u",ylab="v", zlab="",expand=1, col=COL, ...)
+ title(titl,cex.main=2)
+ par(fig=c(0.825,1,0,1))
+ mini <- findInterval(x=min(k,na.rm=TRUE),vec=M)
+ maxi <- findInterval(x=max(k,na.rm=TRUE),vec=M)
+ legend("right",fill=colo(n)[maxi:mini],legend=M[maxi:mini],horiz=F,bty="n")
+ }
+ else
+ {
+ par(cex.lab=2,cex.axis=1.5,font=2,lwd=1)
+ persp(x=K$dist, y=K$times, z=k, xlab="u",ylab="v", zlab="", expand=1, col=COL, ...)
+ title(titl,cex.main=2)
+ }
+ }
+ else
+ {
+ if(isTRUE(legend))
+ {
+ par(cex.lab=1.5,cex.axis=1.5,font=2,lwd=1,plt=c(0,1,0,1),mar=c(0.5,0.5,2.5,.5),las=1)
+ par(fig=c(0.1,0.825,0.1,1))
+ contour(K$dist, K$times, k, labcex=1.5,levels=M,drawlabels=F,col=colo(n),zlim=range(M),axes=F)
+ box(lwd=2)
+ at <- axTicks(1)
+ axis(1,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(1,at=at[length(at)],labels="u",cex.axis=2)
+ at <- axTicks(2)
+ axis(2,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(2,at=at[length(at)],labels="v",cex.axis=2)
+ title(titl,cex.main=2)
+ par(fig=c(0,1,0.1,1))
+ mini <- findInterval(x=min(k,na.rm=TRUE),vec=M)
+ maxi <- findInterval(x=max(k,na.rm=TRUE),vec=M)
+ legend("right",fill=colo(n)[maxi:mini],legend=M[maxi:mini],horiz=F,bty="n")
+ }
+ else
+ {
+ par(cex.lab=2,cex.axis=1.5,font=2,lwd=2,las=1)
+ contour(K$dist, K$times, k, labcex=1.5,levels=M,drawlabels=T,col=colo(n),zlim=range(M),axes=F)
+ box(lwd=2)
+ at <- axTicks(1)
+ axis(1,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(1,at=at[length(at)],labels="u",cex.axis=2)
+ at <- axTicks(2)
+ axis(2,at=at[1:(length(at)-1)],labels=at[1:(length(at)-1)])
+ axis(2,at=at[length(at)],labels="v",cex.axis=2)
+ title(titl,cex.main=2)
+ }
+ }
+par(old.par)
+}
+
Modified: pkg/R/rinfec.r
===================================================================
--- pkg/R/rinfec.r 2012-06-04 12:29:36 UTC (rev 55)
+++ pkg/R/rinfec.r 2012-08-15 13:40:25 UTC (rev 56)
@@ -139,7 +139,7 @@
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)
+ ITK <- findInterval(vec=t.grid$times,x=ti)
continue <- FALSE
while(continue==FALSE)
@@ -173,7 +173,7 @@
}
else
{
- itk <- .iplace(t.grid$times,tk,t.grid$tinc)
+ itk <- findInterval(vec=t.grid$times,x=tk)
if ((mu[itk]==0) || ((h=="gaussian") && (mu[itk]<ug)))
{
@@ -224,8 +224,8 @@
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)
+ nix <- findInterval(vec=s.grid$x,x=xp)
+ niy <- findInterval(vec=s.grid$y,x=yp)
Up <- Lambda[nix,niy]/lmax
retain <- up <= Up
if ((retain==FALSE) || is.na(retain))
Modified: pkg/R/rlgcp.r
===================================================================
--- pkg/R/rlgcp.r 2012-06-04 12:29:36 UTC (rev 55)
+++ pkg/R/rlgcp.r 2012-08-15 13:40:25 UTC (rev 56)
@@ -351,9 +351,9 @@
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)
+ nix <- findInterval(vec=s.grid$x,x=x[ix])
+ niy <- findInterval(vec=s.grid$y,x=y[ix])
+ nit <- findInterval(vec=t.grid$times,x=times[ix])
prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
@@ -404,9 +404,9 @@
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)
+ nix <- findInterval(vec=s.grid$x,x=wx[ix])
+ niy <- findInterval(vec=s.grid$y,x=wy[ix])
+ nit <- findInterval(vec=t.grid$times,x=wtimes[ix])
prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
M <- which(is.na(prob))
Modified: pkg/R/rpp.r
===================================================================
--- pkg/R/rpp.r 2012-06-04 12:29:36 UTC (rev 55)
+++ pkg/R/rpp.r 2012-08-15 13:40:25 UTC (rev 56)
@@ -1,543 +1,534 @@
-.make.grid <- function(nx,ny,poly)
- {
- if (missing(poly)) poly <- matrix(c(0,0,1,0,1,1,0,1),4,2,T)
-
- if ((nx < 2) || (ny < 2)) stop("the grid must be at least of size 2x2")
-
- xrang <- range(poly[, 1], na.rm = TRUE)
- yrang <- range(poly[, 2], na.rm = TRUE)
- xmin <- xrang[1]
- xmax <- xrang[2]
- ymin <- yrang[1]
- ymax <- yrang[2]
-
- xinc <- (xmax-xmin)/nx
- yinc <- (ymax-ymin)/ny
-
- xc <- xmin-xinc/2
- yc <- ymin-yinc/2
- xgrid <- rep(0,nx)
- ygrid <- rep(0,ny)
- xgrid[1] <- xc + xinc
- ygrid[1] <- yc + yinc
-
- for (i in 2:nx)
- {
- xgrid[i] <- xgrid[i-1]+xinc
- }
- for (i in 2:ny)
- {
- ygrid[i] <- ygrid[i-1]+yinc
- }
-
- yy <- matrix(xgrid,nx,ny)
- xx <- t(yy)
- yy <- matrix(ygrid,nx,ny)
-
- X <- as.vector(xx)
- Y <- as.vector(yy)
-
- poly <- rbind(poly,poly[1,])
- pts <- inpip(pts=cbind(X,Y),poly)
-
- X[pts] <- TRUE
- X[X!=TRUE] <- FALSE
- mask <- matrix(X,ncol=ny,nrow=nx,byrow=TRUE)
-
- invisible(return(list(x=xgrid,y=ygrid,X=xx,Y=yy,pts=pts,xinc=xinc,yinc=yinc,mask=matrix(as.logical(mask),nx,ny))))
- }
-
-.rhpp <- function(lambda, s.region, t.region, npoints=NULL, replace=TRUE, discrete.time=FALSE)
- {
- 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)
-
- s.area <- areapl(s.region)
- t.region <- sort(t.region)
- t.area <- t.region[2]-t.region[1]
-
- if (missing(lambda) & !(is.null(npoints)))
- {
- if (t.area==0) lambda <- npoints/(s.area)
- else lambda <- npoints/(s.area * t.area)
- }
-
- pattern <- list()
- index.t <- list()
- if (is.numeric(lambda))
- {
- if (is.null(npoints)==TRUE)
- {
- if (t.area==0)
- { npoints <- round(rpois(n=1,lambda=lambda * s.area),0) }
- else
- { npoints <- round(rpois(n=1,lambda=lambda * s.area * t.area),0) }
- }
- xy <- matrix(csr(poly=s.region,npoints=npoints),ncol=2)
- x <- xy[,1]
- y <- xy[,2]
- npoints <- length(x)
-
- if (discrete.time==TRUE)
- {
- vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
- if ((length(vect)<npoints) & (replace==FALSE))
- stop("when replace=FALSE and discrete.time=TRUE, the length of seq(t.region[1],t.region[2],by=1) must be greater than the number of points")
- names(vect) <- 1:length(vect)
- M <- sample(vect,npoints,replace=replace)
- times <- M
- names(times) <- NULL
- samp <- as.numeric(names(M))
- }
- else
- {
- times <- runif(npoints,min=t.region[1],max=t.region[2])
- samp <- sample(1:npoints,npoints,replace=replace)
- times <- times[samp]
- }
- times <- sort(times)
- index.times <- sort(samp)
- pattern.interm <- list(x=x,y=y,t=times,index.t=index.times)
- pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)
- index.t <- pattern.interm$index.t
- }
- else
- stop("lambda must be numeric")
-
- invisible(return(list(pts=pattern,index.t=index.t)))
- }
-
-
-.iplace <- function(X,x,xinc)
- {
- n <- length(X)
- i <- 0
- repeat
- {
- i <- i+1
- if ((x >= X[i]-xinc) & (x < X[i]+xinc))
- break
- }
- ip <- i
- return(ip)
- }
-
-
-.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)
-
- s.area <- areapl(s.region)
- t.region <- sort(t.region)
- t.area <- t.region[2]-t.region[1]
-
- if (missing(lambda) & !(is.null(npoints)))
- {
- if (t.area==0) lambda <- npoints/(s.area)
- else lambda <- npoints/(s.area * t.area)
- }
-
- lambdamax <- lmax
- pattern <- list()
- index.t <- list()
-
- if (is.function(lambda))
- {
- 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)))
-
- if (is.null(Lambda))
- {
- 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=TRUE)
- M[!(s.grid$mask)] <- NaN
- Lambda[,,it] <- M
- }
- }
-
- if (is.null(npoints)==TRUE)
- {
- if (t.area==0)
- { en <- sum(Lambda,na.rm=TRUE)*s.grid$xinc*s.grid$yinc }
- else
- {
- en <- sum(Lambda,na.rm=TRUE)*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=TRUE)
- npts <- round(lambdamax/(s.area*t.area),0)
- if (npts==0) stop("there is no data to thin")
-
- xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
- x <- xy[,1]
- y <- xy[,2]
-
- 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)
- times <- times.init[samp]
- prob <- lambda(x,y,times,...)/lambdamax
- u <- runif(npts)
- retain <- u <= prob
- if (sum(retain==FALSE)==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=TRUE)}}
- lambdamax <- max(lambdas,na.rm=TRUE)
- 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]
- samp.remain <- (1:nt)[-samp]
- times <- times[retain]
-
- neffec <- length(x)
- 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(replace==FALSE)
- { 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
- u <- runif(npoints-neffec)
- 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 <- list(x=x,y=y,t=times,index.t=index.times)
- pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)
- index.t <- pattern.interm$index.t
- }
-
- if (is.character(lambda))
- {
- if (is.null(Lambda))
- stop("Lambda must be specified")
- nx <- dim(Lambda)[1]
- ny <- dim(Lambda)[2]
- nt <- dim(Lambda)[3]
-
- 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)))
-
- if (is.null(npoints))
- {
- en <- sum(Lambda,na.rm=TRUE)*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=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)
- 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(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,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==F)==length(retain)) retain.eq.F <- FALSE
- else retain.eq.F <- TRUE
- }
-# if (sum(retain==FALSE)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
-
-
- if (sum(retain==FALSE)==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=TRUE)}}
- lambdamax <- max(lambdas,na.rm=TRUE)
- prob <- lambda(x,y,times,...)/lambdamax
- retain <- u <= prob
- if (sum(retain==FALSE)==length(retain)) stop ("no point was retained at the first iteration, please check your parameters")
- }
-
- 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(replace==FALSE)
- { 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 <- NULL
- for(nx in 1:length(wx))
- {
- 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,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 <- list(x=x,y=y,t=times,index.t=index.times)
- pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)
- index.t <- pattern.interm$index.t
- }
-
- invisible(return(list(pts=pattern,index.t=index.t)))
-}
-
-
-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, ...)
-{
-
- 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)
-
- s.area <- areapl(s.region)
- t.region <- sort(t.region)
- t.area <- t.region[2]-t.region[1]
-
- if (missing(lambda) & !(is.null(npoints)))
- {
- if (t.area==0) lambda <- npoints/(s.area)
- else lambda <- npoints/(s.area * t.area)
- }
-
- lambdamax <- lmax
- pattern <- list()
- index.t <- list()
- ni <- 1
-
- #
- # Homogeneous Poisson Process
- #
-
- if (is.numeric(lambda))
- {
- while(ni<=nsim)
- {
- hpp <- .rhpp(lambda=lambda, s.region=s.region, t.region=t.region, npoints=npoints, replace=replace, discrete.time=discrete.time)
- if (nsim==1)
- {
- pattern <- as.3dpoints(hpp$pts)
- index.t <- hpp$index.t
- }
- else
- {
- pattern[[ni]] <- as.3dpoints(hpp$pts)
- index.t[[ni]] <- hpp$index.t
- }
- ni <- ni+1
- }
- Lambda <- NULL
- }
-
- #
- # Inhomogeneous Poisson Process
- #
-
- if (is.function(lambda))
- {
- 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)))
-
- 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=TRUE)
- M[!(s.grid$mask)] <- NaN
- Lambda[,,it] <- M
- }
-
- 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, ...)
-
- if (nsim==1)
- {
- pattern <- as.3dpoints(ipp$pts)
- index.t <- ipp$index.t
- }
- else
- {
- pattern[[ni]] <- as.3dpoints(ipp$pts)
- index.t[[ni]] <- ipp$index.t
- }
- ni <- ni+1
- }
- }
-
- if (is.character(lambda))
- {
- 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, ...)
-
- if (nsim==1)
- {
- pattern <- as.3dpoints(ipp$pts)
- index.t <- ipp$index.t
- }
- else
- {
- pattern[[ni]] <- as.3dpoints(ipp$pts)
- index.t[[ni]] <- ipp$index.t
- }
- ni <- ni+1
- }
- }
-
- invisible(return(list(xyt=pattern,index.t=index.t,s.region=s.region,t.region=t.region,lambda=lambda,Lambda=Lambda)))
-}
-
-
-
+.make.grid <- function(nx,ny,poly)
+ {
+ if (missing(poly)) poly <- matrix(c(0,0,1,0,1,1,0,1),4,2,T)
+
+ if ((nx < 2) || (ny < 2)) stop("the grid must be at least of size 2x2")
+
+ xrang <- range(poly[, 1], na.rm = TRUE)
+ yrang <- range(poly[, 2], na.rm = TRUE)
+ xmin <- xrang[1]
+ xmax <- xrang[2]
+ ymin <- yrang[1]
+ ymax <- yrang[2]
+
+ xinc <- (xmax-xmin)/nx
+ yinc <- (ymax-ymin)/ny
+
+ xc <- xmin-xinc/2
+ yc <- ymin-yinc/2
+ xgrid <- rep(0,nx)
+ ygrid <- rep(0,ny)
+ xgrid[1] <- xc + xinc
+ ygrid[1] <- yc + yinc
+
+ for (i in 2:nx)
+ {
+ xgrid[i] <- xgrid[i-1]+xinc
+ }
+ for (i in 2:ny)
+ {
+ ygrid[i] <- ygrid[i-1]+yinc
+ }
+
+ yy <- matrix(xgrid,nx,ny)
+ xx <- t(yy)
+ yy <- matrix(ygrid,nx,ny)
+
+ X <- as.vector(xx)
+ Y <- as.vector(yy)
+
+ poly <- rbind(poly,poly[1,])
+ pts <- inpip(pts=cbind(X,Y),poly)
+
+ X[pts] <- TRUE
+ X[X!=TRUE] <- FALSE
+ mask <- matrix(X,ncol=ny,nrow=nx,byrow=TRUE)
+
+ invisible(return(list(x=xgrid,y=ygrid,X=xx,Y=yy,pts=pts,xinc=xinc,yinc=yinc,mask=matrix(as.logical(mask),nx,ny))))
+ }
+
+.rhpp <- function(lambda, s.region, t.region, npoints=NULL, replace=TRUE, discrete.time=FALSE)
+ {
+ 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)
+
+ s.area <- areapl(s.region)
+ t.region <- sort(t.region)
+ t.area <- t.region[2]-t.region[1]
+
+ if (missing(lambda) & !(is.null(npoints)))
+ {
+ if (t.area==0) lambda <- npoints/(s.area)
+ else lambda <- npoints/(s.area * t.area)
+ }
+
+ pattern <- list()
+ index.t <- list()
+ if (is.numeric(lambda))
+ {
+ if (is.null(npoints)==TRUE)
+ {
+ if (t.area==0)
+ { npoints <- round(rpois(n=1,lambda=lambda * s.area),0) }
+ else
+ { npoints <- round(rpois(n=1,lambda=lambda * s.area * t.area),0) }
+ }
+ xy <- matrix(csr(poly=s.region,npoints=npoints),ncol=2)
+ x <- xy[,1]
+ y <- xy[,2]
+ npoints <- length(x)
+
+ if (discrete.time==TRUE)
+ {
+ vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
+ if ((length(vect)<npoints) & (replace==FALSE))
+ stop("when replace=FALSE and discrete.time=TRUE, the length of seq(t.region[1],t.region[2],by=1) must be greater than the number of points")
+ names(vect) <- 1:length(vect)
+ M <- sample(vect,npoints,replace=replace)
+ times <- M
+ names(times) <- NULL
+ samp <- as.numeric(names(M))
+ }
+ else
+ {
+ times <- runif(npoints,min=t.region[1],max=t.region[2])
+ samp <- sample(1:npoints,npoints,replace=replace)
+ times <- times[samp]
+ }
+ times <- sort(times)
+ index.times <- sort(samp)
+ pattern.interm <- list(x=x,y=y,t=times,index.t=index.times)
+ pattern <- cbind(x=pattern.interm$x,y=pattern.interm$y,t=pattern.interm$t)
+ index.t <- pattern.interm$index.t
+ }
+ else
+ stop("lambda must be numeric")
+
+ invisible(return(list(pts=pattern,index.t=index.t)))
+ }
+
+
+
+.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)
+
+ s.area <- areapl(s.region)
+ t.region <- sort(t.region)
+ t.area <- t.region[2]-t.region[1]
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/stpp -r 56
More information about the Stpp-commits
mailing list