| RFparameters {RandomFields} | R Documentation |
RFparameters sets and returns control parameters for the simulation
of random fields
RFparameters(..., no.readonly=FALSE)
... |
arguments in |
no.readonly |
If |
The possible parameters are
General options
PracticalRangelogical or .
If not FALSE the range of the primitive
covariance functions is
adjusted so that cov(1) is zero for models with finite range.
The value of
cov(1) is about 0.05 (for scale=1)
for models without range.
See CovarianceFct or type
PrintModelList()
for the list of
primitive models.
FALSE : the practical range ajustment is not used.
TRUE : PracticalRange is applicable only if
the value is known exactly, or, at least, can be approximated by
a closed formula.
2 : if the practical range is not known exactly it
is approximated numerically.
Default: FALSE [init].
PrintLevelIf PrintLevel<=0
there is not any output on the screen. The
higher the number the more tracing information is given.
Default: 1 [init, do].
1 : error messages
2 : messages about partial failures of the algorithm
>2 : additional informations
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.
General options for simulating
pchCharacter or empty string.
The character is printed after each
performed simulation if more than one simulation is performed at
once. If 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].
StoringLogical.
If 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 performed with DoSimulateRF for the same
model parameters are performed faster.
See alse CE.several, TBMCE.several and
local.several for related parameters.
Default: FALSE [do].
skipcheckslogical.
If 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.onlyLogical or NA. Used for the automatic choice of methods. See also RFMethods.
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].
exactnesslogical or NA.
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].
everyinteger.
if greater than zero, then every everyth iteration is
printed if simulated by TBM or random coin method.
Default: 0 [do].
anisological.
If TRUE missing anisotropy parameter
or scale parameter in model $ are
intepreted as identity matrix. If FALSE
missing anisotropy or scale parameter is intepreted
as scale parameter of value 1.0. The latter coding tends to faster
simulation, but also to more error messages. The parameter should
be changed only by advanced users.
Default: TRUE [init].
Options for simulating with the standard circulant embedding method
CE.forceLogical. Circulant embedding does not work if a
certain circulant
matrix has negative eigenvalues. Sometimes it is convenient
to replace all the negative eigenvalues by zero
(CE.force=TRUE) after CE.trials number of trials.
Default: FALSE [init].
CE.mminScalar or vector, integer if positive.
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 |\code{CE.mmin}|; here CE.mmin must be smaller
than -1; the
value -1 takes over the minimal starting size.
Note: in any cases, the initial size might be increased according
to CE.useprimes.
Default: 0 [init].
CE.strategy0 : if the circulant
matrix has negative eigenvalues then the
size in each direction is doubled;
1 : the size is enhanced only in
one direction, namely that one where the covariance function has the
largest value at the end point of the grid — note that
the default value of 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.maxmemmaximal total size of the circulant matrix.
The total amount of memory needed for the internal calculations
is about 16 (=2 * sizeof(double))
times as large as 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.tolImIf the modulus of the imaginary part is less than
CE.tolIm then the eigenvalue is considered as real.
Default: 1E-3 [init].
CE.tolReEigenvalues between CE.tolRe and 0 are considered as
0 and set 0. Default: -1E-7 [init].
CE.trialsA larger circulant matrix is likely to make more eigenvalues
non-negative. If at least one of the thresholds 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.severallogical.
If FALSE only half the memory is need, but
only a single independent realisation can created.
Default: TRUE [init].
CE.useprimesLogical. If 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.dependentLogical. If 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].
CE.method0 or 1. Decomposition of the covariance matrix. This parameter is only relevant if multivariate random fields are simulated.
Default (CE.method=0):
If .Random.seed is fixed, Cholesky decomposition allows for
fixing the random field of the first component whilst
parameters for the second field are changed.
If CE.method=0 and approximate circulant embedding is
allowed, i.e. CE.force=TRUE,
SVD is tried for the last failed attempt of Cholesky
decomposition. (SVD, in contrast to Cholesky decomposition, allows
for approximate circulant embedding.)
If CE.method=1, Cholesky
decomposition will not be attempted, but singular value
decomposition used instead in all attempts. (SVD is slower, but
more precise.)
Default: 0 [init].
Options for simulating by simple matrix decomposition
direct.bestvariablesinteger.
When searching for an appropriate simuation method
the matrix decomposition method (see ‘direct method’ below)
is preferred if the number of variables is less than or equal to
direct.bestvariables.
Default: 800 [init]
direct.maxvariablesIf the number of variables to generate is
greater than 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.methodDecomposition of the covariance matrix.
If direct.method=1, Cholesky
decomposition will not be attempted, but singular value
decomposition
used instead if direct.svdtolerance positive.
In case of a multivariate random field, direct.method = 2
or 3 orders the covariance such that first all components are
considered for the first variable, then all components for the
second one, and so on. If direct.method = 0 or 1
it starts with the first component of all locations, then the
second components follow, etc.
Default: 0 [init].
direct.svdtoleranceIf SVD decomposition is used for calculating the square root of
the covariance matrix then the absolute componentwise difference between
the covariance matrix and square of the square root must be less
than direct.svdtolerance. No check is performed if
direct.svdtolerance is negative.
Default: 1e-12 [init].
Options for simulating hyperplane tessellations
hyper.superposinteger.
number of superposed hyperplane tessellations.
Default: 300 [do].
hyper.maxlinesinteger.
Maximum number of allowed lines.
Default: 1000 [init].
hyper.mar.distrinteger. code for the marginal distribution used in the simulation:
0uniform distribution
1Frechet distribution with form parameter
hyper.mar.param
2Bernoulli distribution (Binomial with n=1) with
parameter hyper.mar.param
The parameter should not be changed yet.
Default: 0 [do].
hyper.mar.paramParameter used for the marginal
distribution. The parameter should not be changed yet.
Default: 0 [do].
Options for simulating with the local ce
methods (cutoff, intrinsic)
local.forcesee CE.force above.
Default: FALSE [init].
local.mminsee CE.mmin above.
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.maxmemsee CE.maxmem above.
Default: 20000000 [init].
local.tolImsee CE.tolIm above.
Default: 1E-7 [init].
local.tolResee CE.tolRe above.
Default: -1E-9 [init].
local.severalsee CE.several above.
Default: 1 [init].
local.useprimessee CE.useprimes above.
Default: TRUE [init].
local.dependentsee CE.dependent above.
Default: FALSE [init].
Options for simulating a Gaussian Markov random fields
markov.neighbours2 or 3.
Number of neighbours in each direction.
Default: 2 [init].
markov.precisionSee also the GMRF manual.
Default: 1.0 [init].
markov.cycliclogical.
simulation on a torus; should not be changed
except the user knows what he/she does. (The covariance function
will be (slightly?!) changed then without further notice!)
Default: FALSE [init].
markov.maxmeminteger.
maximum numer of points of the grid.
Default: 250000 [init].
Options for simulating by random coins
mpp.locationsinteger. Way of defining the point process.
-1 : automatic choice
0 : grid
1 : Poisson point process
2 : uniform distribution
3 : single point follows d dimensional Gaussian distribution
Default: -1 [init].
mpp.intensityreal.
Number of superposed
realisations (to approximate the normal distribution; total number
for all (additive) components with same anisotropy);
if mpp.intensity<=0.0 then only a single
value is simulated (for checking).
Default: 100 [do].
mpp.plusIn order avoid edge effects, the simulation area is enlarged by
a constant r so that all marks have their
(supposed) support in the ball with radius r centred at
the origin; see also mpp.approxzero.
If mpp.plus>0 the true radius r is replaced by
mpp.radius.
Default: 0.0 [init].
mpp.relRadius is added to the effectiveRadius.
The latter gives the approximate radius for the support of the
grain.
Default: 2.5 [init].
mpp.approxzeroFunctions that
do not have
compact support are set to zero outside the ball outside for which the
function has absolute values less than MPP.approxzero.
Default: 0.001 [init].
mpp.samplingdistmpp.samplingdist gives the grid distance
to numerically approximate the volume of the grain.
[-mpp.samplingr,mpp.sampling]^d gives the
domain over which the approximation is calculated.
Default: 0.01 [init].
mpp.samplingrsee mpp.samplingdist
Default: 5 [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.tolpoints at a distance less than or equal to 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].
nugget.methlogical. If TRUE any method given the user for the
simulation of the nugget effect is replaced by the method
‘nugget’ whenever appropriate (zonal nugget or a method
different form circulant embedding, direct and sequential).
Options for using the sequential method
The method works only for spatio-temporal settings (and
grids).
sequ.maxinteger.
maximum number of allowed variables on which
the conditional distribution is conditionned.
Default: 5000 [init].
sequ.backinteger.
number of previous instances on which
the algorithm should condition.
If less than one then the number of previous instances
equals sequ.max / (number of spatial points).
Default: 5 [init].
sequ.initialFirst, N=(number of spatial points) * sequ.back
number of points are simulated. Then, sequentially,
all spatial points for the next time instance
are simulated at once, based on the previous sequ.back
instances. The distribution of the first N points
is the correct distribution, but
differs, in general, from the distribution of the sequentially
simulated variables. We prefer here to have the same distribution
all over (although only approximatively the correct one),
hence do some initial sequential steps first.
If sequ.initial is non-negative, then sequ.initial
first steps are performed.
If sequ.initial is negative, then
sequ.back - sequ.initial
initial steps are performed. The latter ensures that
none of the very first N variables are returned.
Default: -10 [init].
Options for special simulation methods
Special methods exist for the following covariance functions
this covariance function is simulated by a certain moving average of a spatially independent, but temporally dependent random field Y.
spec.MA.rgives the radius beyond which the bivariate
standard normal density is considered as being zero
Default: 2.5 [init].
spec.MA.distthe random field Y is
approximated by a grid; the grid length is given by
spec.MA.dist.
Default: 0.1 [init].
Options for simulating with a turning bands method
Currently, there are 3 variants of the turning bands method
implemented:
spectralThe spectral turning bands method is implemented for 2 (and 1) dimensions only.
TBM2It is based on the two dimensional turning bands operator and is applicable for 1 and 2 dimensions. As an additional dimension the time dimension can be added.
TBM3It is based on the three dimensional turning bands operator and is applicable for 1,2,3 dimensions. As an additional dimension the time dimension can be added.
The following parameters are used.
spectral.gridLogical.
The angle of the lines is random if
spectral.grid=FALSE,
and k*pi/spectral.lines
for k in 1:spectral.lines,
otherwise. This parameter is only considered
if the spectral measure, not the density is used.
Default: TRUE [do].
spectral.ergodicIn case of an additive model and spectral.ergodic=FALSE,
the additive component are chosen proportional to their
variance. In total spectral.lines are simulated. If
spectral.ergodic=TRUE, the components are simulated
separately and then added.
Default: FALSE [do].
spectral.linesNumber of lines used (in total for all additive components of the
covariance function). Default: 500 [do].
spectral.metroLogical. If TRUE then preference is given
to the (slower) metropolis sampling algorithm of the spectral
density.
Default: FALSE [init].
spectral.nmetrointeger. Considered if the Metropolis
algorithm is used. It gives the number of metropolis steps
before returning the next draw from the spectral density.
If spectral.nmetro
is not positive thenRandomFields tries to find a good
choice for
spectral.nmetro by itself.
Default: 0 [init].
spectral.sigmareal. Considered if the Metropolis
algorithm is used. It gives the standard deviation of the
multivariate normal distribution of the proposing
distribution.
If spectral.sigma
is not positive thenRandomFields tries to find a good
choice for
spectral.sigma itself.
Default: 0 [init].
TBM.methodcharacter.
The preferred method to simulate on the line for TBM2 and
TBM3;
If ‘Nothing’ then automatic choice.
Default: "Nothing" [init].
TBM.centerScalar or vector.
If not 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.pointsinteger. If greater than 0,
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.everyIf TBM2.every>0 then every
TBM2.everyth iteration is announced.
Default: 0 [do].
TBM2.linesNumber of lines used.
Default: 60 [do].
TBM2.linesimufactorTBM2.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.linesimustepIf TBM2.linesimustep is positive the grid on the line has lag
TBM2.linesimustep.
See also TBM2.linesimufactor.
Default: 0.0 [init].
TBM2.layersLogical or integer. If 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.everyIf TBM3.every>0 then every
TBM3.everyth iteration is announced.
Default: 0 [do].
TBM3.linesNumber of lines used.
Default: 500 [do].
TBM3.linesimufactorSee TBM2.linesimufactor for
the meaning.
Default: 2.0 [init].
TBM3.linesimustepSee
TBM2.linesimustep for the meaning.
Default: 0.0 [init].
TBM3.layersSee
TBM2.layers for the meaning.
Default: FALSE [init].
TBMCE.forcesee TBM.method and CE.force
Default: FALSE [init].
TBMCE.mminsee TBM.method and
CE.mmin. Default: 0 [init].
TBMCE.strategysee TBM.method and
CE.strategy. Default: 0 [init].
TBMCE.maxmemsee TBM.method and
CE.maxmem. Default: 10000000 [init].
TBMCE.tolImsee TBM.method and
CE.tolIm. Default: 1E-3 [init].
TBMCE.tolResee TBM.method and
CE.tolRe. Default: -1E-7 [init].
TBMCE.trialssee TBM.method and
CE.trials. Default: 3 [init].
TBMCE.useprimessee TBM.method and
CE.useprimes. Default: TRUE [init].
TBMCE.dependentsee TBM.method and CE.dependent.
Default: FALSE [init].
Options specific to simulating max-stable random fields
maxstable.maxGaussMax-stable random fields.
The simulation of the max-stable process based on random fields uses
a stopping rule that necessarily needs a finite upper endpoint
of the marginal distribution of the random field.
In the case of extremal Gaussian random fields,
see 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/~schlather
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, i.e. grid, non-grid) 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)
z1 <- GaussRF(x1, y1, model=model, param=param, grid=TRUE,
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)
tbm.points <- length(GetRegisterInfo(meth="TBM")$S$line)
z2 <- GaussRF(x2, x2, model=model, param=param, grid=TRUE, register=1,
method=meth, TBM.center=0, TBM.points=tbm.points)
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("stable", alpha=1.0),
list(model="gencauchy", alpha=1.0, beta=2.0)
)
for (exactness in c(NA, FALSE, TRUE)) {
readline(paste("\n\nexactness: `", exactness, "'; press return"))
z <- GaussRF(x, x, grid=TRUE, gridtriple=FALSE,
model=model, exactness=exactness,
stationary.only=NA, Print=1, n=1,
TBM2.linesimustep=1, Storing=TRUE)
Print(GetRegisterInfo()$meth$name)
}
#############################################################
## 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 for n independentent repetitions ##
## and 2*n locally dependent dependent repetitions . ##
## To get the precision, the procedure is repeated m times.##
#############################################################
# In the example below, local.dependent speeds up the simulation
# by about factor 16 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("$", var=1.1, aniso=matrix(nc=2, c(2,1,0.5,1)),
list(model="exp"))
(CovarianceFct(matrix(c(1, -1), ncol=2), model=model)) ## true value
m <- if (interactive()) 100 else 2
n <- if (interactive()) 100 else 10
# using local.dependent=FALSE (which is the default)
c1 <- numeric(m)
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="")
c1[i] <- cov(z[1,length(y), ], z[length(x), 1, ])
}) # many times slower than with local.dependent=TRUE
mean(c1)
sd(c1)
# using local.dependent=TRUE...
c2 <- numeric(m)
unix.time(
for (i in 1:m) {
cat("", i)
z <- GaussRF(x, y, model=model, grid=TRUE, method="cu", n=2 * n,
local.dependent=TRUE, pch="")
c2[i] <- cov(z[1,length(y),], z[length(x), 1 , ])
})
mean(c2)
sd(c2) # the sd is samller (using more locally dependent realisations)
## but it is (much) faster! For n=n2 instead of n=2 * n, the
## value of sd(c2) would be larger due to the local dependencies
## in the realisations.