RFparameters(Storing=TRUE)
str(RFparameters())
############################################################
## ##
## use of TBM.points and TBM.center ##
## ##
############################################################
## The following example shows that the same realisation
## can be obtained on different grid geometries (or point
## configurations) using TBM3 (or TBM2)
x1 <- seq(-150,150,1)
y1 <- seq(-15, 15, 1)
x2 <- seq(-50, 50, 1)
model <- "exponential"
param <- c(0, 1, 0, 10)
meth <- "TBM3"
###### simulation of a random field on long thing stripe
runif(1)
rs <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
DeleteAllRegisters()
z1 <- GaussRF(x1, y1, model=model, param=param, grid=TRUE, register=1,
method=meth, TBM.center=0, Storing=TRUE)get(getOption("device"))(height=1.55, width=12)
par(mar=c(2.2, 2.2, 0.1, 0.1))
image(x1, y1, z1, col=rainbow(100))
polygon(range(x2)[c(1,2,2,1)], range(y1)[c(1,1,2,2)], border="red", lwd=3)
###### definition of a random field on a square of shorter diagonal
assign(".Random.seed", rs, envir=.GlobalEnv)
z2 <- GaussRF(x2, x2, model=model, param=param, grid=TRUE, register=2,
method=meth, TBM.center=0,
TBM.points=length(GetRegisterInfo(1)$method[[1]]$mem$l))get(getOption("device"))(height=4.3, width=4.3)
par(mar=c(2.2, 2.2, 0.1, 0.1))
image(x2, x2, z2, zlim=range(z1), col=rainbow(100))
polygon(range(x2)[c(1,2,2,1)], range(y1)[c(1,1,2,2)], border="red", lwd=3)
############################################################
## ##
## use of exactness ##
## ##
############################################################
x <- seq(0, 1, 1/30)
model <- list(list(model="stable", var=1, scale=1, kappa=1.0),
"+",
list(model="gencauchy", var=1, scale=1, kappa=c(1, 2)),
)
for (exactness in c(NA, FALSE, TRUE)) {
readline(paste("exactness: ", exactness, "; press return"))
DeleteRegister()
z <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE,
model=model, exactness=exactness,
stationary.only=NA, Print=4, n=1,
TBM2.linesimustep=1, Storing=TRUE)
str(lapply(GetRegisterInfo()$method,
function(x) x[c("name", "covnr")]))
}
#############################################################
## The following gives a tiny example on the advantage of ##
## local.dependent=TRUE (and CE.dependent=TRUE) if in a ##
## study most of the time is spent with simulating the ##
## Gaussian random fields. Here, the covariance at a pair ##
## of points is estimated. ##
#############################################################
# In the example below, local.dependent speeds up the simulation
# by about factor 27 at the price of an increased variance of
# factor 1.5
x <- seq(0, 1, len=10)
y <- seq(0, 1, len=10)
grid.size <- c(length(x), length(y))
model <- list(list(model="exp", var=1.1, aniso=c(2,1,0.5,1)))
CovarianceFct(matrix(c(1, -1), ncol=2), model=model) ## true value
RFparameters(Storing=TRUE)
m <- if (interactive()) 1000 else 2
# determine number of non-overlapping realisations on the torus
DeleteRegister()
InitGaussRF(x, y, model=model, grid=TRUE, method="cu")
blocks <- GetRegisterInfo()$method[[1]]$mem$new$method[[1]]$mem$simupos
(n <- prod(blocks) * 1) ## n any multiple of prod(blocks) to avoid
## dependencies between the m estimated covariance if
## if local.dep=TRUE; or put RFparameters(Storing=FALSE),
## but this is slower
# using local.dependent=TRUE...
c1 <- numeric(m)
DeleteRegister()
unix.time(
for (i in 1:m) {
cat("", i)
z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=n,
local.dependent=TRUE, pch="")
c1[i] <- cov(z[1,length(y),], z[length(x), 1 , ])
}) # about 35 0.3 35 0 0
var(c1) # about 0.013
# using local.dependent=FALSE...
c2 <- numeric(m)
DeleteRegister()
unix.time(
for (i in 1:m) {
cat("", i)
z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=n,
local.dependent=FALSE, pch="")
c2[i] <- cov(z[1,length(y),], z[length(x), 1 , ])
}) # about 950 3 950 0 0
var(c2) # about 0.0087Run the code above in your browser using DataLab