#############################################################
## ##
## Examples using the symmetric stable model, also called ##
## "powered exponential model" ##
## ##
#############################################################
PrintModelList() ## the complete list of implemented models
model <- "stable"
mean <- 0
variance <- 4
nugget <- 1
scale <- 10
alpha <- 1 ## see help("CovarianceFct") for additional
## parameters of the covariance functions
step <- 1 ## nicer, but also time consuming if step <- 0.1
x <- seq(0, 20, step)
y <- seq(0, 20, step)
f <- GaussRF(x=x, y=y, model=model, grid=TRUE,
param=c(mean, variance, nugget, scale, alpha))
image(x, y, f)
#############################################################
## ... using gridtriple
step <- 1 ## nicer, but also time consuming if step <- 0.1
x <- c(0, 20, step) ## note: vectors of three values, not a
y <- c(0, 20, step) ## sequence
f <- GaussRF(grid=TRUE, gridtriple=TRUE,
x=x ,y=y, model=model,
param=c(mean, variance, nugget, scale, alpha))
image(seq(x[1],x[2],x[3]), seq(y[1],y[2],y[3]), f)
#############################################################
## arbitrary points
x <- runif(100, max=20)
y <- runif(100, max=20)
z <- runif(100, max=20) # 100 points in 3 dimensional space
(f <- GaussRF(grid=FALSE,
x=x, y=y, z=z, model=model,
param=c(mean, variance, nugget, scale, alpha)))
#############################################################
## usage of a specific method
## -- the complete list can be obtained by PrintMethodList()
x <- runif(100, max=20) # arbitrary points
y <- runif(100, max=20)
(f <- GaussRF(method="dir", # direct matrix decomposition
x=x, y=y, model=model, grid=FALSE,
param=c(mean, variance, nugget, scale, alpha)))
#############################################################
## simulating several random fields at once
step <- 1 ## nicer, but also time consuming if step <- 0.1
x <- seq(0, 20, step) # grid
y <- seq(0, 20, step)
f <- GaussRF(n=3, # three simulations at once
x=x, y=y, model=model, grid=TRUE,
param=c(mean, variance, nugget, scale, alpha))
image(x, y, f[,,1])
image(x, y, f[,,2])
image(x, y, f[,,3])
#############################################################
## ##
## Examples using the extended definition forms ##
## ##
## ##
#############################################################
## note that the output seems plausible but not checked!!!!
## tbm may also be used for multiplicate models (if they have
## *exactly* the same anisotropy parameters)
x <- (0:100)/10
m <- matrix(c(1,2,3,4),ncol=2)/5
z <- GaussRF(x=x, y=x, grid=TRUE,
model=list(
list(m="power",v=1,k=2,a=m),
"*", list(m="sph", v=1, a=m)
),
me="TBM3", reg=0,n=1)
print(c(mean(as.double(z)),var(as.double(z))))
image(z,zlim=c(-3,3))
## non-separable space-time model applied for two space dimensions
## note that tbm method does not work nicely, but at least
## in some special cases.
x <- y <- (1:32)/2 ## grid definition, but as a sequence
T <- c(1,32,1)*10 ## note necessarily gridtriple definition
aniso <- diag(c(0.5,8,1))
k <- c(1,phi=1,1,0.5,psi=1,dim=2)
model <- list(list(m="nsst", v=1, k=k, a=aniso))
z <- GaussRF(x=x, y=y, T=T, grid=TRUE, model=model)
rl <- function() if (interactive()) readline("Press return")
for (i in 1:dim(z)[3]) { image(z[,,i]); rl();}
for (i in 1:dim(z)[2]) { image(z[,i,]); rl();}
for (i in 1:dim(z)[1]) { image(z[i,,]); rl();}
#############################################################
## ##
## Example of a 2d random field based on ##
## covariance functions valid in 1d only ##
## ##
#############################################################
x <- seq(0, 10, 1/10)
model <- list(list(model="fractgauss", var=1, kappa=0.5,
aniso=c(1, 0, 0, 0)),
"*",
list(model="fractgauss", var=1, kappa=0.5,
aniso=c(0, 0, 0, 1)))
z <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE, model=model)
image(x, x, z)
#############################################################
## ##
## Brownian motion ##
## (using Stein's method) ##
## ##
#############################################################
# 2d
step <- 0.3 ## nicer, but also time consuming if step <- 0.1
x <- seq(0, 10, step)
kappa <- 1 # in [0,2)
z <- GaussRF(x=x, y=x, grid=TRUE, model="fractalB",
param=c(0,1,0,1,kappa))
image(z,zlim=c(-3,3))
# 3d
x <- seq(0, 3, step)
kappa <- 1 # in [0,2)
z <- GaussRF(x=x, y=x, z=x, grid=TRUE, model="fractalB",
param=c(0,1,0,1,kappa))
rl <- function() if (interactive()) readline("Press return")
for (i in 1:dim(z)[1]) { image(z[i,,]); rl();}
#############################################################
## This example shows the benefits from stored, ##
## intermediate results: in case of the circulant ##
## embedding method, the speed is doubled in the second ##
## simulation. ##
#############################################################
DeleteAllRegisters()
RFparameters(Storing=TRUE, PrintLevel=1)
y <- x <- seq(0, 50, 0.2)
(p <- c(runif(3), runif(1)+1))
ut <- system.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
method="circ", param=p))
image(x, y, f)
hist(f)
c( mean(as.vector(f)), var(as.vector(f)) )
cat("system time (first call)", format(ut,dig=3),"")
# second call with the *same* parameters is much faster:
ut <- system.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
method="circ",param=p))
image(x, y, f)
hist(f)
c( mean(as.vector(f)), var(as.vector(f)) )
cat("system time (second call)", format(ut,dig=3),"")
#############################################################
## ##
## Example how the cutoff method can be set ##
## explicitly using hypermodels ##
## ##
#############################################################
## NOTE: this feature is still in an experimental stage
## which has not been yet tested intensively;
## further: parameters and algorithms may change in
## future.
## simuation of the stable model using the cutoff method
#RFparameters(Print=8, Storing=FALSE)
x <- seq(0, 1, 1/24)
scale <- 1.0
model1 <- list( list(model="stable", var=1, scale=scale, kappa=1.0) )
rs <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
z1 <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE,
model=model1, n=1, meth="cutoff", Storing=TRUE)
size <- GetRegisterInfo()$method[[1]]$mem$new$method[[1]]$mem$size
#str(GetRegisterInfo(), vec=15)
## simulation of the same random field using the circulant
## embedding method and defining the hypermodel explicitely
model2 <- list(list(model="cutoff", var=1, kappa=c(1, sqrt(2), 1),
scale=scale),
"(",
list(model="stable", var=1, scale=scale, kappa=1.0)
)
assign(".Random.seed", rs, envir=.GlobalEnv)
z2 <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE,
model=model2, n=1, meth="circulant", CE.mmin=size)print(range(z1-z2)) ## essentially no difference between the fields!
#############################################################
## The cutoff method simulates on a torus and a (small) ##
## rectangle is taken as the required simulation. ##
## ##
## The following code shows a whole such torus. ##
## The main part of the code sets local.dependent=TRUE and ##
## local.mmin to multiples of the basic rectangle lengths ##
#############################################################
# definition of the realisation
x <- seq(0, 2, len=20)
y <- seq(0, 1, len=40)
grid.size <- c(length(x), length(y))
model <- list(list(model="exp", var=1.1, aniso=c(2,1,0.5,1)))
# determination of the (minimal) size of the torus
DeleteRegister()
InitGaussRF(x, y, model=model, grid=TRUE, method="cu")
ce.info <- GetRegisterInfo()$method[[1]]$mem$new$method[[1]]$mem
blocks <- ceiling(ce.info$size / grid.size)
size <- blocks * grid.size
# simulation and plot of the torus
DeleteRegister()
z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=prod(blocks),
local.dependent=TRUE, local.mmin=size)
hei <- 8
do.call(getOption("device"),
list(hei=hei, wid=hei / blocks[2] / diff(range(y)) *
blocks[1] * diff(range(x))))
close.screen(close.screen())
split.screen(rev(blocks))
k <- 0
for (j in 1:blocks[2]) {
for (i in 1:blocks[1]) {
k <- k + 1
screen(k)
par(mar=rep(1, 4) * 0.02)
image(z[,,(blocks[2]-j) * blocks[1] + i], zlim=c(-3, 3), axes=FALSE)
}
}Run the code above in your browser using DataLab