Nothing
.set.cov <- function(separable,model,param,sigma2)
{
mods <- 0
modt <- 0
mod <- 0
models <- c("exponential","cauchy","stable","wave","gneiting","cesare","matern","none")
for(i in 1:length(model))
{
M <- which(models==model[i])
if (length(M)==0) stop("the model is not implemented")
}
for (i in 1:length(unique(model)))
{
if (((isTRUE(separable)) & ((model[i]==models[5]) | (model[i]==models[6]))) | ((!(isTRUE(separable))) & ((model[i]==models[1]) | (model[i]==models[2]) | (model[i]==models[3]) | (model[i]==models[4]) | (model[i]==models[7])))) stop("'stcov' does not match with 'model'")
}
if (isTRUE(separable))
{
if ((length(model)!=1) & (length(model)!=2))
stop("for separable covariance functions, 'model' must be of length 1 or 2")
if (length(model)==1)
{
if (model=="none")
{
mods <- 0
modt <- 0
}
if (model=="exponential")
{
mods <- 1
modt <- 1
}
if (model=="stable")
{
mods <- 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("Stable model parameter must lie in [0,2]")
}
if (model=="cauchy")
{
mods <- 3
if (param[1]<=0) stop("Cauchy model parameter must be strictly positive")
modt <- 3
if (param[2]<=0) stop("Cauchy model parameter must be strictly positive")
}
if (model=="wave")
{
mods <- 4
modt <- 4
}
if (model=="matern")
{
mods <- 7
if (param[2]<=0 | param[1]<=0) stop("Matern model parameters must be strictly positive")
modt <- 7
if (param[3]<=0 | param[4]<=0) stop("Matern model parameters must be strictly positive")
}
}
if (length(model)==2)
{
if (model[1]=="none")
mods <- 0
if (model[2]=="none")
modt <- 0
if (model[1]=="exponential")
mods <- 1
if (model[2]=="exponential")
modt <- 1
if (model[1]=="stable")
{
mods <- 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("Stable model parameter must lie in [0,2]")
}
if (model[1]=="cauchy")
{
mods <- 3
if (param[1]<=0) stop("Cauchy model parmaeter must be strictly positive")
}
if (model[2]=="cauchy")
{
modt <- 3
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[3]<=0 | param[1]<=0) stop("Matern model parameters must be strictly positive")
}
if (model[2]=="matern")
{
modt <- 7
if (param[2]<=0 | param[4]<=0) stop("Matern model parameters must be strictly positive")
}
}
}
if (!(isTRUE(separable)))
{
if (length(model)!=1)
stop("for non-separable covariance functions, 'model' must be of length 1")
if (model=="gneiting")
{
mod <- 5
if (param[6]<2) stop("for Gneiting's covariance function, the sixth parameter must be greater than 2")
if ((param[3]<=0) | (param[3]>2)) stop("for Gneiting's covariance function, the third parameter must lie in (0,2]")
if ((param[4]<=0) | (param[4]>1)) stop("for Gneiting's covariance function, the fourth parameter must lie in (0,1]")
if ((param[5]!=1) & (param[5]!=2) & (param[5]!=3)) stop("for Gneiting's covariance function, the fifth parameter must be 1, 2 or 3")
if ((param[2]!=1) & (param[2]!=2) & (param[2]!=3)) stop("for Gneiting's covariance function, the second parameter must be 1, 2 or 3")
if ((param[2]==1) & ((param[1]<0) | (param[1]>2))) stop("for Gneiting's covariance function, if the second parameter equals 1, the first parameter must lie in [0,2]")
if ((param[2]==2) & (param[1]<=0)) stop("for Gneiting's covariance function, if the second parameter equals 2, the first parameter must be strictly positive")
}
if (model=="cesare")
{
mod <- 6
if (((param[1]>2) | (param[1]<1)) | ((param[2]>2) | (param[2]<1))) stop("for De Cesare's model, the first and second parameters must lie in [1,2]")
if (param[3]<3/2) stop("for De Cesare's model, the third parameter must be greater than 3/2")
}
}
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)
{
nt <- length(times)
np <- length(dist)
model <- .set.cov(separable,model,param,sigma2)
gs <- array(0, dim = c(np,nt))
storage.mode(gs) <- "double"
gs <- .Fortran("covst",
(gs),
as.double(dist),
as.integer(np),
as.double(times),
as.integer(nt),
as.integer(model),
as.double(param),
as.double(sigma2),
as.double(scale))[[1]]
if (plot==TRUE)
{
image(dist,times,gs,col=grey((1000:1)/1000),xlab="h",ylab="t",cex.axis=1.5,cex.lab=2,font=2)
contour(dist,times,gs,add=T,col=4,labcex=1.5,nlevels=nlevels)
}
return(gs)
}
.gauss3D <- function(nx=100,ny=100,nt=100,xlim,ylim,tlim,separable=TRUE,model="exponential",param=c(1,1,1,1,1,2),scale=c(1,1),var.grf=1,mean.grf=0,exact=TRUE)
{
N <- c(nx,ny,nt)
mod <- .set.cov(separable,model,param,var.grf)
g <- floor(log(2*(N-1))/log(2))+1
M <- 2^g
count <- 1
changG <- TRUE
while(changG==TRUE)
{
L <- rep(-9999,M[1]*M[2]*M[3])
storage.mode(L) <- "double"
res <- .Fortran("circ",
(L),
as.integer(M),
as.integer(N),
as.double(xlim),
as.double(ylim),
as.double(tlim),
as.integer(mod),
as.double(param),
as.double(var.grf),
as.double(scale))[[1]]
L <- array(res,dim=M)
FTL <- fft(L)
if (isTRUE(exact))
{
if (min(Re(FTL))<0)
{
g <- g+1
M <- 2^g
changG <- TRUE
count <- count+1
}
else
changG <- FALSE
}
else
{
FTL[Re(FTL)<0] <- 0
changG <- FALSE
}
}
print(count)
X <- array(rnorm(M[1]*M[2]*M[3],0,1),dim=M)
X <- fft(X,inverse=TRUE)
A <- sqrt(FTL)*X
G <- (Re(fft(A,inverse=FALSE))/(M[1]*M[2]*M[3]))[1:N[1],1:N[2],1:N[3]]
G <- G+mean.grf
## Remarks
# --------
#
# The right way is:
# X1 <- array(rnorm(M[1]*M[2]*M[3],0,1),dim=M)
# X2 <- array(rnorm(M[1]*M[2]*M[3],0,1),dim=M)
# X <- fft(X1+1i*X2,inverse=TRUE)
# A <- sqrt(FTL)*X
# G <- (Re(fft(A,inverse=FALSE))/(M[1]*M[2]*M[3]))[1:N[1],1:N[2],1:N[3]]
# but it doesn't change the results and it is slightly faster.
#
# Taking the first elements of FTL and then computing G (changing M by N)
# is faster than taking the first elements of G, but it does not provide
# the same results
return(G)
}
rlgcp <- function(s.region, t.region, replace=TRUE, npoints=NULL, nsim=1, nx=100, ny=100, nt=100,separable=TRUE,model="exponential",param=c(1,1,1,1,1,2),scale=c(1,1),var.grf=1,mean.grf=0,lmax=NULL,discrete.time=FALSE,exact=FALSE)
{
if (missing(s.region)) s.region <- matrix(c(0,0,1,1,0,1,1,0),ncol=2)
if (missing(t.region)) t.region <- c(0,1)
t.region <- sort(t.region)
s.area <- areapl(s.region)
t.area <- t.region[2]-t.region[1]
tau <- c(start=t.region[1],end=t.region[2],step=(t.region[2]-t.region[1])/(nt-1))
bpoly <- bbox(s.region)
lambdamax <- lmax
pattern <- list()
index.t <- list()
Lambdafin <- list()
ni <- 1
s.grid <- .make.grid(nx,ny,s.region)
s.grid$mask <- matrix(as.logical(s.grid$mask),nx,ny)
if (discrete.time==TRUE)
{
vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
if (nt>length(vect))
{
nt <- length(vect)
warning("nt used is less than the one given in argument")
t.grid <- list(times=vect,tinc=1)
}
else
{
vect <- round(seq(floor(t.region[1]),ceiling(t.region[2]),length=nt))
t.grid <- list(times=vect,tinc=round(t.area/(nt-1)))
}
}
else
t.grid <- list(times=seq(t.region[1],t.region[2],length=nt),tinc=(t.area/(nt-1)))
while(ni<=nsim)
{
S <- .gauss3D(nx=nx,ny=ny,nt=nt,xlim=range(s.region[,1]),ylim=range(s.region[,2]),tlim=range(t.region),separable=separable,model=model,param=param,scale=scale,var.grf=var.grf,mean.grf=mean.grf,exact=exact)
Lambda <- exp(S)
mut <- rep(0,nt)
for (it in 1:nt)
{
Lambda[,,it][s.grid$mask==FALSE] <- NaN
mut[it] <- sum(Lambda[,,it],na.rm=TRUE)
}
if (is.null(npoints))
{
en <- sum(Lambda,na.rm=TRUE)*s.grid$xinc*s.grid$yinc*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")
if ((replace==FALSE) & (nt < max(npts,npoints))) stop("when replace=FALSE, nt must be greater than the number of points used for thinning")
if (discrete.time==TRUE)
{
vect <- seq(floor(t.region[1]),ceiling(t.region[2]),by=1)
times.init <- sample(vect,nt,replace=replace)
}
else
times.init <- runif(nt,min=t.region[1],max=t.region[2])
samp <- sample(1:nt,npts,replace=replace,prob=mut/max(mut,na.rm=TRUE))
times <- times.init[samp]
retain.eq.F <- FALSE
while(retain.eq.F==FALSE)
{
xy <- matrix(csr(poly=s.region,npoints=npts),ncol=2)
x <- xy[,1]
y <- xy[,2]
prob <- NULL
for(ix in 1:length(x))
{
nix <- 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])
if (nix==0 | niy==0 | nit==0)
prob=c(prob,NA)
else
prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
M <- which(is.na(prob))
if (length(M)!=0)
{
x <- x[-M]
y <- y[-M]
times <- times[-M]
prob <- prob[-M]
npts <- length(x)
}
u <- runif(npts)
retain <- u <= prob
if (sum(retain==FALSE)==length(retain)) retain.eq.F <- FALSE
else retain.eq.F <- TRUE
}
x <- x[retain]
y <- y[retain]
samp <- samp[retain]
samp.remain <- (1:nt)[-samp]
times <- times[retain]
neffec <- length(x)
if (neffec > npoints)
{
retain <- 1:npoints
x <- x[retain]
y <- y[retain]
samp <- samp[retain]
samp.remain <- (1:nt)[-samp]
times <- times[retain]
}
while(neffec < npoints)
{
xy <- as.matrix(csr(poly=s.region,npoints=npoints-neffec))
if (dim(xy)[2]==1) {wx <- xy[1]; wy <- xy[2]}
else {wx <- xy[,1]; wy <- xy[,2]}
if (isTRUE(replace))
wsamp <- sample(1:nt,npoints-neffec,replace=replace,prob=mut/max(mut,na.rm=TRUE))
else
wsamp <- sample(samp.remain,npoints-neffec,replace=replace,prob=mut[samp.remain]/max(mut[samp.remain],na.rm=TRUE))
wtimes <- times.init[wsamp]
prob <- NULL
for(ix in 1:length(wx))
{
nix <- 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])
if (nix==0 | niy==0 | nit==0)
prob=c(prob,NA)
else
prob <- c(prob,Lambda[nix,niy,nit]/lambdamax)
}
M <- which(is.na(prob))
if (length(M)!=0)
{
wx <- wx[-M]
wy <- wy[-M]
wtimes <- wtimes[-M]
prob <- prob[-M]
}
if (neffec > 0)
{
u <- runif(length(prob))
retain <- u <= prob
x <- c(x,wx[retain])
y <- c(y,wy[retain])
times <- c(times,wtimes[retain])
samp <- c(samp,wsamp[retain])
samp.remain <- (1:nt)[-samp]
neffec <- length(x)
}
}
times <- sort(times)
index.times <- sort(samp)
pattern.interm <- cbind(x=x,y=y,t=times)
if (nsim==1)
{
pattern <- as.3dpoints(pattern.interm)
index.t <- index.times
Lambdafin <- Lambda
}
else
{
pattern[[ni]] <- as.3dpoints(pattern.interm)
index.t[[ni]] <- index.times
Lambdafin[[ni]] <- Lambda
}
ni <- ni+1
}
invisible(return(list(xyt=pattern,s.region=s.region,t.region=t.region,Lambda=Lambdafin,index.t=index.t)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.