fitvario(x, y=NULL, z=NULL, T=NULL, data, model, param,
lower=NULL, upper=NULL, sill=NA, ...)fitvario.default(x, y=NULL, z=NULL, T=NULL, data, model, param,
lower=NULL, upper=NULL, sill=NA, trend,
use.naturalscaling=TRUE, PrintLevel=RFparameters()$Print,
optim.control=NULL, bins=20, nphi=1, ntheta=1, ntime=20,
distance.factor=0.5,
upperbound.scale.factor=10, lowerbound.scale.factor=20,
lowerbound.scale.LS.factor=5,
upperbound.var.factor=10, lowerbound.var.factor=100,
lowerbound.sill=1E-10, scale.max.relative.factor=1000,
minbounddistance=0.001, minboundreldist=0.02,
approximate.functioncalls=50, refine.onborder=TRUE,
pch=RFparameters()$pch,
var.name="X",
time.name="T", transform=NULL, standard.style=NULL,
lsq.methods=c("self", "plain", "sqrt.nr", "sd.inv", "internal"),
mle.methods=c("ml", "reml"),
cross.methods=NULL, users.guess=NULL, users.beta=NULL, table.format=FALSE)
coord;
If also a time component is given, the in the data the indices for
the spatial components run the fastest.modelparam=c(mean, variance, nugget, scale,...);
the parameters must be given
in this order. Further parameters are to be added in case of a
parametrised class of covariance functions, slower and param are vectors and length(lower) <
length(param)
then lower must match the number of additional parameters
$\kappa$.
NA the sill is kept fix. Only used if the
standard format for the covariance model is given. See Details.missing(param);
linear formula : uses X1, X2,... and T as internal
parameters for the coordinates; all parameters are estimated
matrix : must have thefitvario.default and listed in the
following.TRUE then internally, rescaled
covariance functions will be used for which
cov(1)$\approx$0.05.
use.naturalscaling has the advantage thatntimes is a vector, then the first component are the
maximum time distance (in units of the grid length T[3]) and the
second component gives the step size (in units of the grid lenrefine.onborder=TRUE and if the result of
any maximum likelihood method or cross validation method
is on a borderline, then the optimisation is redone
in a modified way (which takes about double extra time)pch!="" then one or two
additional steps in the MLE methods are
marked by transform allows for the definition of a parameter as a
function of other estimated or fixed parameters.
All the parameters are supposed to be in a vector called NULL. This variable should only be
set by the advanced user. If NULL then standard.style will be
TRUE if the covariance model allows for a param (except that no NA's should
be contained) or model.table.format=TRUE then a matrix is returned where the
first rows contain the estimated parameters, followed by
the target values of all methods for the given set of parameters;
the last rows give the lower and upper bounds used in the estimations.
The columns correspond to the various estimation methods for the parameters. If table.format=FALSE then
the function returns a list with the following elements
fitvario in RandomFields V1.0 to
geoR whose homepage is at
First, by means of various control parameters, see below, the algorithm
first tries to estimate the bounds for the parameters to be estimated,
if the bounds for the parameters are not given.
Independently whether users.guess is given,
the algorithm guesses initial values for the parameters.
The automatic guess and the user's guess will be called
primitive methods in the following.
Second, the variogram model is fitted by various least squares
methods (according to the value of lsq.methods) using
the best parameter set among the primitive methods as initial value
if the effective number of parameters is greater than 1.
[Remarks: (i)
Third, the model is fitted by various maximum likelihood methods
(according to the value of mle.methods) using
the best parameter set among the primitive methods and the lsq methods
as initial value (if the effective number of parameters is greater
than 1).
Comments on specific parameters:
lower
The lower bounds are technical bounds that should not
restrict really restrict the domaine of the value.
However, if these values are too small the optimisation
algorithm will frequently run into local minima or get
stuck close the border of the parameter domain.
It is advised to limit seriously the domain of the additional
parameters of the covariance model and/or the total number of
parameters to be estimated, iflowermay containNAs, then these values
are generated by the If a nested model is given, the bounds may again be supplied for all parameters or only for the additional form parameters of the model. The bounds given apply uniformely to all submodels of the nested model.
If themodelis given in list format, thenloweris a list, where components may be missing
orNA. These are generated by the algorithm, then.
IflowerisNULLall lower bounds are generated
automatically.
upper.kappa
Seelower.kappa.sill
Additionally to estimatingnuggetandvarianceseparately, they may also be estimated together under the
condition thatnugget+variance=sill.
For the latter a finite value forsillhas to be supplied,
andnuggetandvarianceare set toNA.sillis only used for the standard model.use.naturalscaling
logical. IfTRUEthen internally, rescaled
covariance functions will be used for which
cov(1)$\approx$0.05. However this parameter
does not influence
the output offitvario: the parameter vector
returned byfitvariorefersalwaysto the standard covariance model as given inPracticalRangeinuse.naturalscaling==TRUE:scaleand the shape parameter of a parameterised
covariance model can be estimated better if they are estimated
simultaneously.upperbound.scale.factorandlowerbound.scale.factor,
etc. might be more realistic.use.naturalscaling==TRUE:fitvario.TRUE.PrintLevel
0 : no message
1 : error messages
2 : warnings
3 : minimum debugging information
5 : extended debugging information, including graphics
Default:0.trace.optim
see control parametertraceof0.bins
vector of explicit boundaries for the bins or the
number of bins for the empirical variogram (used in the
LSQ target function, which is described at the beginning
of the Details).
Note that for anisotropic models, the value ofbinsmight
be enlarged.
Default:20.distance.factor
right boundary of the last
bin is calculated asdistance.factor* (maximum distance
between all pairs of points). Only used ifbinsis a scalar.
Default:0.5.upperbound.scale.factor
The upper bound for the scale is determined
asupperbound.scale.factor* (maximum distance
between all pairs of points).
Default:10.lowerbound.scale.factor
The lower bound for the scale is determined
as$\hbox{(minimum distance
between different pairs of points)} /
\code{lowerbound.scale.factor}$.
Default:20.
lowerbound.scale.LS.factor
For the LSQ target function a different lower bound
for the scale is used. It is determined as$\hbox{(minimum distance
between different pairs of points)} /
\code{lowerbound.scale.LS.factor}$.
Default:5.
upperbound.var.factor
The upper bound for the variance and the nugget is determined
as$$\code{upperbound.var.factor} *
{\hbox{var}}(\code{data}).$$Default:10.lowerbound.var.factor
The lower bound for the variance and the nugget is determined
as$$\hbox{var}(\code{data}) /
\code{lowerbound.var.factor}.$$If a standard model definition is given and
either the nugget or the variance is fixed,
the parameter to be estimated
must also be greater thanlowerbound.sill.
If a non-standard model definition is given
thenlowerbound.var.factoris only used
for the first model; the other lower bounds for the
variance are zero.
Default:100.lowerbound.sill
Seelowerbound.var.factor.
Default:1E-10.scale.max.relative.factor
If the initial scale value for the ML estimation
obtained by the LSQ target function is
less than$\hbox{(minimum distance
between different pairs of points)} /
\code{scale.max.relative.factor}$a warning is given that probably a nugget effect
is present.
Note: ifscale.max.relative.factoris greater
thanlowerbound.scale.LS.factorthen
no warning is given as
the scale has the lower bound (minimum distance
between different pairs of points) /lowerbound.scale.LS.factor.
Default:1000.minbounddistance
If any value of the parameter vector
returned from the ML estimation
is closer thanminbounddistanceto any of the bounds or if any value
has a relative distance smaller thanminboundreldist, then it is assumed that
the MLE algorithm has dropped into a local minimum,
and it will be continued with evaluating the
ML target function on a grid, cf. the beginning paragraphs
of the Details.
Default:0.001.minboundreldist
Seeminbounddistance.
Default:0.02.approximate.functioncalls
In case the parameter vector is too close to the given
bounds, the ML target function is evaluated on a grid
to get a new initial value for the ML estimation.
The number of points of the grid is approximatelyapproximate.functioncalls.
Default:50.lsq.methods
Variogram fit by least squares methods; first, a preliminary
trend is estimated
by a simple regression; second, the variogram is fitted; third,
the trend is fitted using the estimated covariance structure."self"weighted lsq. Weights are the values of the
fitted variogram to the power of -2"plain"model fitted by least squares; trends are never taken
into account"sqrt.nr"weighted lsq. Weight is the square root
of the number
of points in the bin"sd.inv"weighted lsq. Weight is the inverse the
standard deviation of the
variogram cloud within the binmle.methods
Model fit by various maximum likelihood methods
(according to the value ofmle.methods) using
the best parameter set among the primitive methods and the lsq methods
as initial value (if the effective number of parameters is greater
than 1). If the best parameter vector of the MLE found so far is too close
to some given bounds, see the specific parameters above, it is
assumed thatrefine.onborder=TRUEthe MLE target function is calculated on a grid, the
best parameter vector is taken, and the optimisation is restarted with
this parameter vector."ml"maximum likelihood;
since ML and REML give the same result if there are not any
covariates, ML is performed in that case, independently whether
it is given or not."reml"restricted maximum likelihoodCressie, N.A.C. (1993)Statistics for Spatial Data.New York: Wiley.
% \item Cross validation methods % Gneiting, T., M. Schlather, B. Huwe (2005) % \emph{} % In preparation.
RandomFields,model <-"gencauchy"
param <- c(0, 1, 0, 1, 1, 2)
estparam <- c(0, NA, 0, NA, NA, 2) ## NA means: "to be estimated"
## sequence in `estparam' is
## mean, variance, nugget, scale, (+ further model parameters)
## So, mean, variance, and scale will be estimated here.
## Nugget is fixed and equals zero.
points <- 100
x <- runif(points,0,3)
y <- runif(points,0,3) ## 100 random points in square [0, 3]^2
d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=10)
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
lower=c(0.1, 0.1), upper=c(1.9, 5), cross.me=NULL))
#########################################################
## The next two estimations give about the same result.
## For the first the sill is fixed to 1.5. For the second the sill
## is reached if the estimated variance is smaller than 1.5
estparam <- c(0, NA, NA, NA, NA, NA)
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
sill=1.5, cross.me=NULL))
estparam <- c(0, NA, NaN, NA, NA, NA)
parampositions(model=model, param=estparam)
f <- function(param) {
param[5] <- max(0, 1.5 - param[1])
return(param)
}
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
transform=f, cross.me=NULL))
#########################################################
## esimation of coupled parameters (kappa1 = kappa2, here)
estmodel <- list(list(model="gencauchy", var=NA, scale=NA,
kappa=list(NA, function(m) m[[1]]$kappa[1])),
"+",
list(model="nugget", var=NA))
system.time(str(fitvario(x=cbind(x,y), data=d, model=estmodel,
cross.me=NULL), vec.len=6))
## much faster implementation, but trickier coding
estparam <- c(0, NA, NA, NA, NA, NaN)
parampositions(model=model, param=estparam)
f <- function(param) {param[3] <- param[2]; param}
system.time(str(fitvario(x=cbind(x,y), data=d, model=model,
param=estparam, transform=f, cross.me=NULL), vec.len=6))
## about as fast as previous coding
## about the same precision as previous coding, or better
## But, a warning is given, since the user may programme
## strange things in this setup, and the programme cannot check it.
system.time(str(fitvario(x=cbind(x,y), data=d, model=model,
param=estparam, transform=f, cross.me=NULL,
standard.style=TRUE), vec.len=6))
#########################################################
## estimation in a anisotropic framework
x <- y <- (1:6)/4
model <- list(list(m="exp", v=1.5, aniso=c(4,2,-2,1)))
z <- GaussRF(x=x, y=y, grid=TRUE, model=model, n=10)
estmodel <- list(list(m="exp", v=NA, aniso=c(NA, NA,-2,1)))
system.time(str(fitvario(as.matrix(expand.grid(x, y)), data=z,
model=estmodel, nphi=20)))Run the code above in your browser using DataLab