RFparameters {RandomFields} | R Documentation |
RFparameters
sets and returns control parameters for the simulation
of random fields
RFparameters(..., no.readonly=FALSE)
... |
arguments in tag = value form, or a list of tagged
values. |
no.readonly |
If RFparameters is called without
parameter then all parameters are returned in a list. If
no.readonly=TRUE then only rewritable parameters are returned.
|
The possible parameters are
General options
PracticalRange
FALSE
the range of the covariance functions is
adjusted so that cov(1) is about 0.05 (for scale=1
).
FALSE
: the practical range ajustment is not used.
TRUE
: PracticalRange
is applicable only if
the value is known exactly.
2
: PracticalRange
is applicable if the value is
known pretty well
3
: PracticalRange
is applicable if the value is
roughly known
11
: if the practical range is not known exactly it
is approximated numerically.
12
: if the practical range is not known pretty well it
is approximated numerically.
13
: if the practical range is even not known
approximately it is approximated numerically.
Note that values beyond FALSE
, TRUE
, and 11
,
are only used for specialists' purposes.
Default: FALSE
[init].
PrintLevel
PrintLevel
<=0
there is not any output on the screen. The
higher the number the more tracing information is given.
Default: 1 [init, do].
Note that PrintLevel
is also used in other packages
as a default, for example in SoPhy
(risk.index
and
create.roots
). The changing of
PrintLevel
here may cause some unexpected effects in these
functions. See the documentation there.
pch
pch='!'
then an absolute
counter is shown instead of the character.
If pch='%'
then a
counter of percentages is shown instead of the character.
Note that also '^H's are printed in
the last two cases,
which may have undesirable interactions with some few other R
functions, e.g. Sweave
.
Default: '*'
[do].
Storing
FALSE
then the intermediate results are
destroyed after the simulation of the random field(s)
or if an error had occured.
On the other hand, if Storing=TRUE
, then
several simulations for the same
model parameters are performed faster, see the examples.
Note that in subsequent calls of GaussRF
intermediate changes of RFparameters with flag "[init]"
do not have any influence on the algorithm, if
Storing=TRUE
.
See alse CE.several
, TBMCE.several
and
local.several
for related parameters.
Default: FALSE
[do].
skipchecks
TRUE
, the check whether the given parameter values
and the dimension are within the allowed range is skipped.
Do not change the value of this variable except you really
know what you do.
Default: FALSE
[init].
stationary.only
TRUE
: the simulation of non-stationary random
fields is refused. In particular, the intrinsic
embedding method is excluded and
the simulation of Brownian motion is rejected.
FALSE
: intrinsic embedding is always allowed,
actually it's the first one considered in the automatic
selection algorithm.
NA
: the simulation of the Brownian motion allowed,
but intrinsic embedding is not used for stationary random fields.
Default: NA
[init].
exactness
TRUE
: add.MPP, hyperplanes and
all turning bands methods are excluded.
If the circulant embedding method is considered as badly
behaved, then the matrix decomposition methods are preferred.
FALSE
: if the circulant embedding method is
considered as badly behaved or the number of points to be
simulated is large, the turning bands methods are
rather preferred.
NA
: approximative non-exact methods are excluded,
i.e. TBM2 if the Abel transform of the covariance function
cannot be given explicitely.
Default: NA
[init].
Options for simulating with the standard circulant embedding method
CE.force
CE.force=TRUE
) after CE.trials
number of trials.
Default: FALSE
[init].
CE.mmin
CE.mmin
determines the initial size of the circulant
matrix. If CE.mmin=0
the minimal starting size is
determined automatically according to the
dimensions of the grid.
If CE.mmin>0
then the absolute starting size is given.
If CE.mmin<0
then the automatically determined
matrix size is multiplied
by |CE.mmin
|; here CE.mmin
must be smaller than -1; the
value -1 takes over the minimal starting size.CE.useprimes
.
Default: 0
[init].CE.strategy
CE.trials
is probably too small
in that case.
In some cases CE.strategy=0
works better, in other cases
CE.strategy=1
. Just try.
Default: 0
[init].
CE.maxmem
CE.maxmem
if RFparameters
()$Storing=FALSE
and 32 (=4 * sizeof(double)) time as large if Storing=TRUE
.
Note that CE.maxmem
can be used to control the automatic
choice of the simulation algorithm. Namely, in case of huge
circulant matrices, other simulation
methods (TBM) are faster and might be preferred be the user.
Default: 4096^2 = 16777216
[init].
CE.tolIm
CE.tolIm
then the eigenvalue is considered as real.
Default: 1E-3
[init].CE.tolRe
CE.tolRe
and 0 are considered as
0 and set 0. Default: -1E-7
[init].CE.trials
CE.tolRe
and
CE.tolIm
are missed then the matrix size is doubled
according to CE.strategy
,
and the matrix is checked again. This procedure is repeated
up to CE.trials-1
times. If there are still negative
eigenvalues, the simulation method fails if CE.force=FALSE
.
Default: 3
[init].
CE.several
FALSE
only half the memory is need, but
only a single independent realisation can created.
Default: TRUE
[init].
CE.useprimes
FALSE
the columns of the circulant matrix
have length 2^k for some k. Otherwise the algorithm
tries to find a nicely factorizable number close to the size of the
given matrix. Default: TRUE
[init].CE.dependent
FALSE
then independent random fields are created. If TRUE
then at least 4 non-overlapping rectangles are taken out of the
the expanded grid defined by the circulant matrix.
These simulations are dependent.
See below for an example.
See CE.trials
for some
more information on the circulant matrix.
Default: FALSE
[init].
Options for simulating with the local ce
methods (cutoff, intrinsic)
local.force
FALSE
[init].
local.mmin
Difference: if local.mmin=0
the automatic determination
of the initial size of the circulant matrix takes into
account an expansion factor. This expansion factor is
intended to make the circulant matrix positive definite
and is either theoretically or numerically known, or guessed.
If the usual strategy of circulant embedding
(doubling the grid sizes) should be taken over then
local.mmin
must
be set to -1
.
Default: 0
[init].
local.maxmem
CE.maxmem
above.
Default: 20000000
[init].local.tolIm
CE.tolIm
above.
Default: 1E-7
[init].local.tolRe
CE.tolRe
above.
Default: -1E-9
[init].local.several
CE.several
above.
Default: 1
[init].local.useprimes
CE.useprimes
above.
Default: TRUE
[init].local.dependent
CE.dependent
above.
Default: FALSE
[init].
Options for simulating by simple matrix decomposition
direct.bestvariables
direct.bestvariables
Default: 800
[init]
direct.maxvariables
direct.maxvariables
, then any matrix decomposition
method is rejected. It is important that this option is set
conveniently to avoid great losses of time during the automatic
search of a simulation method (method=NULL
in
GaussRF
).
Default: 4096
[init]direct.method
direct.method=1
, Cholesky
decomposition will not be attempted, but singular value
decomposition
used instead.
Default: 0
[init].direct.svdtolerance
direct.svdtolerance
. No check is performed if
direct.svdtolerance
is negative.
Default: 1e-12
[init].
Options for simulating nugget effects
Simulating a nugget effect seems trivial. It gets complicated
and best methods (including direct
and circulant
embedding
!) fail if zonal anisotropies are considered,
where sets of points have to be identified that belong to the
same subspace of eigenvalue 0 of the anisotropy matrix.
nugget.tol
nugget.tol
are considered as being identical. This strategy applies to
the simulation method and the covariance function itself.
Hence, the covariance function is only positive definite
if nugget.tol=0.0
. However, if the anisotropy matrix
does not have full rank and nugget.tol=0.0
then,
the simulations are likely to be odd.
The value of nugget.tol
should be of order 1e-15.
Default: 0.0
[init].
Options for simulating with a turning bands method
Currently, there are 3 variants of the turning bands method
implemented:
spectral
TBM2
TBM3
The following parameters are used.
spectral.grid
spectral.grid=FALSE
,
and k*pi/spectral.lines
for k in 1:spectral.lines
,
otherwise. Default: TRUE
[do].spectral.lines
500
[do].TBM.method
TBM2
and
TBM3
;
currently either 'circulant embedding'
or 'direct'
.
If 'direct'
then the method is overwritten if the number of points on the grid
is larger than direct.maxvariables
.
If the circulant embedding method is used, then the TBMCE
parameters below determine the behaviour of the circulant
embedding algorithm.
Default: "circulant embedding"
[init].
TBM.center
NA
, the TBM.center
is used as the center of
the turning bands for TBM2
and TBM3
.
Otherwise the center is determined
automatically such that the line length is minimal.
See also TBM.points
and the examples below.
Default: NA
[init].
TBM.points
TBM.points
gives the number of points simulated on the TBM
line, hence
must be greater than the minimal number of points given by
the size of the simulated field and the two paramters
TBMx.linesimufactor
and TBMx.linesimustep
.
If TBM.points
is not positive the number of points is
determined automatically.
The use of TBM.center
and TBM.points
is highlighted
in an example below.
Default: 0
[init].
TBM2.every
TBM2.every>0
then every
TBM2.every
th iteration is announced.
Default: 0
[do].TBM2.lines
60
[do].TBM2.linesimufactor
TBM2.linesimufactor
or
TBM2.linesimustep
must be non-negative; if
TBM2.linesimustep
is positive then TBM2.linesimufactor
is ignored.
If both
parameters are naught then TBM.points
is used (and must be
positive).
The grid on the line is TBM2.linesimufactor
-times
finer than the smallest distance.
See also TBM2.linesimustep
.
Default: 2.0
[init].TBM2.linesimustep
TBM2.linesimustep
is positive the grid on the line has lag
TBM2.linesimustep
.
See also TBM2.linesimufactor
.
Default: 0.0
[init].TBM2.num
TRUE
then the covariance function on the line
is approximated numerically.
If FALSE
only those models are allowed
that have an analytic representation on
the line.
Default: TRUE
[init].TBM2.layers
TRUE
then the turning layers are used whenever
a time component is given.
If FALSE
the turning layers are used only when the
traditional TBM is not applicable.
If negative then turning layers may never be used.
If greater than 1 then only turning layers may be used.
Default: FALSE
[init].TBM3.every
TBM3.every>0
then every
TBM3.every
th iteration is announced.
Default: 0
[do].TBM3.lines
500
[do].TBM3.linesimufactor
TBM2.linesimufactor
for
the meaning.
Default: 2.0
[init].TBM3.layers
TBM2.layers
for the meaning.
Default: FALSE
[init].TBM3.linesimus
TBM2.linesimustep
for the meaning.
Default: 0.0
[init].TBMCE.force
TBM.method
and CE.force
Default: FALSE
[init].TBMCE.mmin
TBM.method
and
CE.mmin
. Default: 0
[init].TBMCE.strategy
TBM.method
and
CE.strategy
. Default: 0
[init].TBMCE.maxmem
TBM.method
and
CE.maxmem
. Default: 10000000
[init].TBMCE.tolIm
TBM.method
and
CE.tolIm
. Default: 1E-3
[init].TBMCE.tolRe
TBM.method
and
CE.tolRe
. Default: -1E-7
[init].TBMCE.trials
TBM.method
and
CE.trials
. Default: 3
[init].TBMCE.useprimes
TBM.method
and
CE.useprimes
. Default: TRUE
[init].TBMCE.dependent
TBM.method
and CE.dependent
.
Default: FALSE
[init].
Options for simulating with Poisson point processes
add.MPP.realisations
100
[do].MPP.approxzero
MPP.approxzero
.
Default: 0.001
[init].MPP.radius
MPP.approxzero
.
If MPP.radius>0
the true radius r is replaced by
MPP.radius
.
Default: 0.0
[init].
Options for simulating hyperplane tessellations
hyper.superpos
300
[do].
hyper.maxlines
1000
[init].
hyper.mar.distr
0
1
hyper.mar.param
2
hyper.mar.param
The parameter should not be changed yet.
Default: 0
[do].
hyper.mar.param
0
[do].
Options specific to simulating max-stable random fields
maxstable.maxGauss
MaxStableRF
, the upper endpoint is
approximated by maxstable.maxGauss
.
Default: 3.0
[init].
General comments on the options
The following refers to the simulation of Gaussian random fields
(InitGaussRF
, GaussRF
), but most
parts also apply
for the simulation of max-stable random fields
(InitMaxStableRF
, MaxStableRF
).
Some of the global parameters determine the basic settings of a
simulation, e.g. direct.method
(which chooses a square
root of a positive definite matrix). The values of
such parameters are read by
InitGaussRF
and stored in an internal register.
Changing
such a parameter between calling InitGaussRF
and calling
DoSimulateRF
or between subsequent calls of
GaussRF
will not have any effect. These parameters have
the flag "[init]".
Parameters like TBM2.lines
(which determines the number of
i.i.d. processes to be simulated on the line)
are only relevant when generating
random numbers. These parameters are read by DoSimulateRF
(or by the second part of GaussRF
), and
are marked by "[do]".
Storing
has an influence on both, InitGaussRF
and
DoSimulateRF
. InitGaussRF
may reserve
more memory if Storing=TRUE
. DoSimulateRF
will
free the register
if Storing=FALSE
, whatever the value of Storing
was
when InitGaussRF
was called.
The distinction between [init] and [do] is also relevant if
GaussRF
is used and called a second time
with the same parameters for the random field and if
RFparameters()$Storing=TRUE
.
Then GaussRF
realises that the second call has the
same random field parameters, and
takes over the stored intermediate results (that have been calculated
with the RFparameters()
at that time). To prevent the use of
stored intermediate results or to take into account intermediate
changes of RFparameters
set RFparameters(Storing=FALSE)
or use
DeleteRegister()
between calls of GaussRF
.
A programme that checks whether the parameters are well
adapted to a specific simulation problem is given as an example of
EmpiricalVariogram()
.
For further details on the implemented methods, see RFMethods.
If any parameter has been given
RFparameters
returns an invisible list of
the given parameters in full name.
Otherwise the complete list of parameters is returned. Further the
values of the following internal readonly variables are returned:
* covmaxchar
: max. name length for variogram/covariance
models
* covnr
: number of currently implemented
variogram/covariance models
* distrmaxchar
: max. name length for a distribution
* distrnr
: number of currently implemented
distributions
* maxdim
: maximum number of dimensions for a
random field
* maxmodels
: maximum number of elementary models in
definition of a complex covariance model
* methodmaxchar
: max. name length for methods
* methodnr
: number of currently implemented simulation methods
Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/institute
Schlather, M. (1999) An introduction to positive definite functions and to unconditional simulation of random fields. Technical report ST 99-10, Dept. of Maths and Statistics, Lancaster University.
GaussRF
,
GetPracticalRange
,
MaxStableRF
,
RandomFields
,
and RFMethods
.
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) do.call(getOption("device"), list(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)) do.call(getOption("device"), list(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("\n\nexactness: ", 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("\n", 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("\n", 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.0087