RandomFields (version 2.0.71)

fitvario: LSQ and Maximum Likelihood Estimation of Random Field Parameters

Description

The function estimates arbitrary parameters of a random field specification with various methods.

Usage

fitvario(x, y=NULL, z=NULL, T=NULL, data, model, param,
         lower=NULL, upper=NULL, sill=NA, grid=!missing(gridtriple),
          gridtriple, ...)

fitvario.default(x, y=NULL, z=NULL, T=NULL, data, model, param, grid=!missing(gridtriple), gridtriple=FALSE, trend = NULL, BC.lambda, ## if missing then no BoxCox-Trafo BC.lambdaLB=-10, BC.lambdaUB=10, lower=NULL, upper=NULL, sill=NA, use.naturalscaling=FALSE, PrintLevel, optim.control=NULL, bins=20, nphi=1, ntheta=1, ntime=20, distance.factor=0.5, upperbound.scale.factor=3, lowerbound.scale.factor=3, 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, minmixedvar=1/1000, maxmixedvar=1000, pch=RFparameters()$pch, transform=NULL, standard.style=NULL, var.name="X", time.name="T", lsq.methods=c("self", "plain", "sqrt.nr", "sd.inv", "internal"), mle.methods=c("ml"), cross.methods=NULL, users.guess=NULL, only.users = FALSE, Distances=NULL, truedim, solvesigma = NA, # if NA then use algorithm -- ToDo allowdistanceZero = FALSE, na.rm = TRUE)

Arguments

x
($n\times2$)-matrix of coordinates, or vector of x-coordinates. All locations must be given explicitely and cannot be passed via a grid definition as in GaussRF
y
vector of y coordinates
z
vector of z coordinates
T
vector of T coordinates; these coordinates are given in triple notation, see GaussRF
data
vector or matrix of values measured at coord; If a matrix is given then the columns are interpreted as independent realisations. If also a time component is given, then in the data the indices for the spatial components run the fa
model
string or list; covariance model, see CovarianceFct and Covariance, or type
param
vector or matrix or NULL. If vector then param=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, s
lower
list or vector. Lower bounds for the parameters. If param is a vector, lower has to be a vector as well and its length must equal the number of parameters to be estimated. The order of param has to be mai
upper
list or vector. Upper bounds for the parameters. See also lower.
sill
If not NA the sill is kept fix. Only used if the standard format for the covariance model is given. See Details.
grid
boolean. Weather coordinates give a grid
gridtriple
boolean. Format, see GaussRF
BC.lambda
a vector of at most two numerical components (just one component corresponds to two identical ones) which are the parameters of the box-cox-transformation: $\frac{x^\lambda_1-1}\lambda_1+\lambda_2$ If the model is univariate, the firs
BC.lambdaLB
lower bound for the first box-cox-parameter
BC.lambdaUB
upper bound for the first box-cox-parameter
trend
If a univariate model is used, the following trend types are possible: number: the constant mean (not to be estimated any more) NA: there is a constant mean to be estimated formula : uses X1, X2,... and T as interna
...
arguments as given in fitvario.default and listed in the following.
use.naturalscaling
logical. Only used if model is given in standard (simple) way. If TRUE then internally, rescaled covariance functions will be used for which cov(1)$\approx$0.05. use.naturalscaling has the advantage that
PrintLevel
level to which messages are shown. See Details.
optim.control
control list for optim, which uses 'L-BFGS-B'. However 'parscale' may not be given.
bins
number of bins of the empirical variogram. See Details.
nphi
scalar or vector of 2 components. If it is a vector then the first component gives the first angle of the xy plane and the second one gives the number of directions on the half circle. If scalar then the first angle is assumed to be zero.
ntheta
scalar or vector of 2 components. If it is a vector then the first component gives the first angle in the third direction and the second one gives the number of directions on the half circle. If scalar then the first angle is assumed to be
ntime
scalar or vector of 2 components. if ntimes 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 len
distance.factor
relative right bound for the bins. See Details.
upperbound.scale.factor
relative upper bound for scale in LSQ and MLE. See Details.
lowerbound.scale.factor
relative lower bound for scale in MLE. See Details.
lowerbound.scale.LS.factor
relative lower bound for scale in LSQ. See Details.
upperbound.var.factor
relative upper bound for variance and nugget. See Details.
lowerbound.var.factor
relative lower bound for variance. See Details.
lowerbound.sill
absolute lower bound for variance and nugget. See Details.
scale.max.relative.factor
relative lower bound for scale below which an additional nugget effect is detected. See Details.
minbounddistance
absolute distance to the bounds below which a part of the algorithm is considered as having failed. See Details.
minboundreldist
relative distance to the bounds below which a part of the algorithm is considered as having failed. See Details.
approximate.functioncalls
approximate evaluations of the ML target function on a grid. See Details.
refine.onborder
logical. If refine.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)
minmixedvar
lower bound for variance in a mixed model; so, the covariance model for mixed model part might be calibrated appropriately
maxmixedvar
upper bound for variance in a mixed model; so, the covariance model for mixed model part might be calibrated appropriately
pch
character shown before evaluating any method; if pch!="" then one or two additional steps in the MLE methods are marked by + and #. Default: "*".
var.name
basic name for the coordinates in the formula of the trend. Default: X
time.name
basic name for the time component in the formula of the trend. Default: X
transform
function. Essentially, transform allows for the definition of a parameter as a function of other estimated parameters. All the parameters are supposed to be in a vector called param where the positions are giv
standard.style
logical or 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 standard definiti
lsq.methods
variants of the least squares fit of the variogram. See Details.
mle.methods
variants of the maximum likelihood fit of the covariance function. See Details.
cross.methods
Not implemented yet.
users.guess
User's guess of the parameters. All the parameters must be given using the same rules as for either param (except that no NA's should be contained) or model.
only.users
boolean. If true then only users.guess is used as a starting point for the fitting algorithms
Distances
alternatively to coordinates x, y, and z the distances themselves can be given. Then truedim must be indicated.
truedim
see Distances
solvesigma
Boolean -- experimental stage! If a mixed effect part is present where the variance has to be estimated, then this variance parameter is solved iteratively within the profile likelihood function, if solvesigma=TRUE.This makes
allowdistanceZero
boolean. If true, then multiple observations are allowed within a single data set. In this case, the coordinates are slightly scattered, so that the points have some tiny distances.
na.rm
boolean -- experimental stage. Only the data may have missing values. If na.rm=TRUE then lines of (repeated) data are deleted if at least one missing value appears. If na.rm=FALSE then the repetitions are treat

Value

  • The function returns a list with the following elements
  • evlist returned by EmpiricalVariogram
  • tableMatrix. 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.
  • lowerboundslower bounds
  • lowerboundsupper bounds
  • transformtransformation function
  • varioobsolete
  • selflist containing
    • model
    {the variogram or covariance model} residuals{NULL} ml.value{the likelihood value for the model}

item

plain, sqrt.nr, sd.inv, internal, ml, reml

code

NULL

Acknowledgement

Thanks to Paulo Ribeiro for hints and comparing the perliminary versions of fitvario in RandomFields V1.0 to likfit of the package geoR whose homepage is at http://www.est.ufpr.br/geoR/.

Details

The optimisations are performed using optimize if one parameter has to be estimated only and optim, otherwise.

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) best with respect to the target value of the respective lsq method; (ii) the effective number of parameters in the optimisation algorithm can be smaller than the number of estimated parameters, since in some cases, some parameters can be calculated explicitely; relevant for the choice between optimize and optim is the effective number of parameters; (iii) optim needs]

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:

  • BC.lambdaIf you want to estimateBC.lambdayou should assert that all data values are positive; otherwise errors will probably occur because of the box-cox-transformation. The second parameter of the box-cox-transformation cannot be estimated since it corresponds to the mean. So the mean should be estimated instead.
  • trendAmong the formes mentioned above it is possible to use just one matrix for the trend instead of a list of identical ones.
  • lower The lower bounds are technical bounds that should not 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, ifmanyparameters of the covariance model are estimated. If the model is given in standard form, the user may supply the lower bounds for the whole parameter vector, or only for the additional form parameters of the model. The lower bound for the mean will be ignored.lowermay 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 inCovarianceFct. (In contrast toPracticalRangeinRFparameters.) Advantages ifuse.naturalscaling=TRUE:
    • scaleand the shape parameter of a parameterised covariance model can be estimated better if they are estimated simultaneously.
    • The estimated bounds calculated by means ofupperbound.scale.factorandlowerbound.scale.factor, etc. might be more realistic.
    • in case of anisotropic models, the inverse of the elements of the anisotropy matrix should be in the above bounds.
    Disadvantages ifuse.naturalscaling=TRUE:
    • For some covariance models with additional parameters, the rescaling factor has to be determined numerically. Then, more time is needed to performfitvario.
    Default: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 parametertraceofoptim. Default:0.
  • 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 bin
  • mle.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 thatoptimran into a local minimum because of a bad starting value. In this case and ifrefine.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 likelihood
    % \item \code{cross.methods}\cr % Here, the parameters are chosen so that % the sum of the distances between the measured and estimated % values in the cross validation are minimised. % The various methods differ in the distance definition. % % The best parameter set among the primitive, lsq and mle methods % are used as initial parameter value in the optimisation % (if the effective number of parameters is greater % than 1). If the best parameter vector of the cross validation % algorithm % found so far is too close % to some given bounds, see the specific parameters above, it is % assumed that % \command{\link[stats]{optim}} ran into a local minimum because of a bad starting % value. % In this case and if \code{refine.onborder=TRUE} % the target function is calculated on a grid, the % best parameter vector is taken, and the optimisation is restarted with % this parameter vector. % % \itemize{ % \item \code{"cross.sq"} squared distance % \item \code{"cross.abs"} absolute value % \item \code{"cross.ign"} \sQuote{ignorance} distance, see % Gneiting et al. (2005) % \item \code{"cross.crps"} \sQuote{xxxxxxx} distance, see % Gneiting et al. (2005) % }
  • References

    • Least squares and mle methods

      Cressie, 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.

    • Related software Ribeiro, P. and Diggle, P. (2001) Software for geostatistical analysis using R and S-PLUS: geoR and geoS, version 0.6.15.http://www.maths.lancs.ac.uk/~ribeiro/geoR.html.
    • REML (rml) LaMotte, L.R. (2007) A direct derivation of the REML likelihood functionStatistical Papers48, 321-327.

    See Also

    Covariance, CovarianceFct, GetPracticalRange, parampositions RandomFields, weather.

    Examples

    Run this code
    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) ## 300 random points in square [0, 3]^2
     ## simulate data according to the model:
     d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=1000) #1000
     ## fit the data:
    
    Print(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
        lower=c(0.1, 0.1, 0.1), upper=c(1.9, 5, 2)))
    
    
    #########################################################
    ## 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) 
    Print(v <- fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
        sill=1, use.nat=FALSE)) ## gencauchy works better with use.nat=FALSE
    
    estmodel <-  list("+",
                      list("$", var=NA, scale=NA,
                           list("gencauchy", alpha=NA, beta=NA)
                           ),
                      list("$", var=NA, list("nugget"))
                     )
    parampositions(model=estmodel, dim=2)
    f <- function(variab) c(variab, max(0, 1.0 - variab[1]))
    Print(v2 <- fitvario(x=cbind(x,y), data=d, model=estmodel, 
                      lower = c(TRUE, TRUE, TRUE, TRUE, FALSE),
                      transform=f, use.nat=FALSE))
    
    #########################################################
    ## estimation of coupled parameters (alpha = beta, here)
    # source("RandomFields/tests/source.R")
    f <- function(param) param[c(1:3,3,4)]
    Print(fitvario(x=cbind(x,y), data=d, model=estmodel,
                 lower=c(TRUE, TRUE, TRUE, FALSE, TRUE),
                 transform=f))
    
    
    
    #########################################################
    ## estimation in a anisotropic framework
    x <- y <- (1:6)/4
    model <- list("$", aniso=matrix(nc=2, c(4,2,-2,1)), var=1.5,
                  list("exp"))
    z <- GaussRF(x=x, y=y, grid=TRUE, model=model, n=10)
    estmodel <-list("$", aniso=matrix(nc=2, c(NA,NA,-2,1)), var=NA,
                    list("exp"))
    Print(fitvario(as.matrix(expand.grid(x, y)), data=z,
                 model=estmodel, nphi=20))
    #########################################################
    ## estimation with trend (formula)
    model <- list("$", var=1, scale=2, list("gauss"))
    estmodel <- list("$", var=NA, scale=NA, list("gauss"))
    x <- seq(-pi,pi,pi/2)
    n <- 5
    data <- GaussRF(x, x, gridtri=FALSE, model=model,
           trend=function(X1,X2) sin(X1) + 2*cos(X2),n=n)
    Print(v <- fitvario(x, x, data=data, gridtrip=FALSE,
                      model=estmodel,
                      trend=~sin(X1)+cos(X1)+sin(X2)+cos(X2)))
    
    
    
    ########################################################
    ## estimation of anisotropy matrix with two identical ##
    ## diagonal elements                                  ##
    x <- c(0, 5, 0.4)
    model <- list("$", var=1, scale=1, list("exponential"))
    z <- GaussRF(x, x, x, model=model, gridtriple=TRUE, n=10, Print=2)
    
    est.model <- list("+",
                     list("$", var=NA, aniso=diag(c(NA, NA, NA)), list("exponen")),
                     list("$", var=NA, list("nugget")))
    parampositions(est.model, dim=3)
    trafo <- function(variab) {variab[c(1:2, 2:4)]}
    lower <- c(TRUE, TRUE, FALSE, TRUE, TRUE) # which parameter to be estimated
    fitlog <- fitvario(x, x, x, gridtriple=TRUE, data=z, model=est.model,
                       transform=trafo, lower=lower)
    str(fitlog$ml)

    Run the code above in your browser using DataCamp Workspace