fields: mKrig.MLE – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

mKrig.MLE

Maximizes likelihood for the process marginal variance (rho) and nugget standard deviation (sigma) parameters (e.g. lambda) over a many covariance models or covariance parameter values.

Description

These functions are designed to explore the likelihood surface for different covariance parameters with the option of maximizing over sigma and rho. They are depreciated and my be omitted in later versions of fields with their roles being replaced by other functions. See details below.

Usage

mKrig.MLE(x, y, weights = rep(1, nrow(x)), cov.fun="stationary.cov",
cov.args = NULL, 
          Z = NULL, par.grid = NULL, lambda = NULL, lambda.profile = TRUE, 
          verbose = FALSE, relative.tolerance = 1e-04, ...)

mKrig.MLE.joint(x, y, weights = rep(1, nrow(x)), 
                lambda.guess = 1, cov.params.guess=NULL, 
                cov.fun="stationary.cov", cov.args=NULL, 
                Z = NULL, optim.args=NULL, find.trA.MLE = FALSE, 
                ..., verbose = FALSE)

fastTps.MLE(x, y, weights = rep(1, nrow(x)), Z = NULL, ...,
                 par.grid=NULL, theta, lambda = NULL, lambda.profile = TRUE,
                 verbose = FALSE, relative.tolerance = 1e-04)

Arguments

`cov.args`	Additional arguments that would also be included in calls to the covariance function to specify the fixed part of the covariance model.
`cov.fun`	The name, a text string, of the covariance function.
`cov.params.guess`	A list of initial guesses for covariance parameters over which the user wishes to perform likelihood maximization. The list contains the names of the parameters as well as the values.
`find.trA.MLE`	If TRUE will estimate the effective degrees of freedom using a simple Monte Carlo method throughout joint likelihood maximization. Either way, the trace of the mKrig object with the best log-likelihood is calculated depending on `find.trA`. Computing the trace will add to the computational burden by approximately NtrA solutions of the linear system but the cholesky decomposition is reused.
`lambda`	If `lambda.profile=FALSE` the values of lambda to evaluate the likelihood if "TRUE" the starting values for the optimization. If lambda is NA then the optimum value from previous search is used as the starting value. If lambda is NA and it is the first value the starting value defaults to 1.0.
`lambda.guess`	The initial guess for lambda in the joint log-likelihood maximization process.
`lambda.profile`	If `TRUE` maximize likelihood over lambda.
`optim.args`	Additional arguments that would also be included in calls to the optim function in joint likelihood maximization. If `NULL`, this will be set to use the "BFGS-" optimization method. See `optim` for more details. The default value is: `optim.args = list(method = "BFGS", control=list(fnscale = -1, ndeps = rep(log(1.1), length(cov.params.guess)+1), reltol=1e-04, maxit=10))` Note that the first parameter is lambda and the others are the covariance parameters in the order they are given in `cov.params.guess`. Also note that the optimization is performed on a log-scale, and this should be taken into consideration when passing arguments to `optim`.
`par.grid`	A list or data frame with components being parameters for different covariance models. A typical component is "theta" comprising a vector of scale parameters to try. If par.grid is "NULL" then the covariance model is fixed at values that are given in ....
`relative.tolerance`	Relative tolerance used to declare convergence when maximizing likelihood over lambda.
`theta`	Range parameter for compact Wendland covariance. (seefastTps)
`verbose`	If `TRUE` print out interesting intermediate results.
`weights`	Precision ( 1/variance) of each observation
`x`	Matrix of unique spatial locations (or in print or surface the returned mKrig object.)
`y`	Vector or matrix of observations at spatial locations, missing values are not allowed! Or in mKrig.coef a new vector of observations. If y is a matrix the columns are assumed to be independent observations vectors generated from the same covariance and measurment error model.
`Z`	Linear covariates to be included in fixed part of the model that are distinct from the default low order polynomial in `x`
`...`	Additional arguments that would also be included in a call to `mKrig` to specify the covariance model and fixed model covariables.

Details

The "mKrig" prefixed functions are depreciated and are replaced in functionality by mKrigMLEJoint and mKrigMLEGrid.

The observational model follows the same as that described in the Krig function and thus the two primary covariance parameters for a stationary model are the nugget standard deviation (sigma) and the marginal variance of the process (rho). It is useful to reparametrize as rho and\ lambda= sigma^2/rho. The likelihood can be maximized analytically over rho and the parameters in the fixed part of the model the estimate of rho can be substituted back into the likelihood to give a expression that is just a function of lambda and the remaining covariance parameters. It is this expression that is then maximized numerically over lambda when lambda.profile = TRUE.

Note that fastTps.MLE is a convenient variant of this more general version to use directly with fastTps, and mKrig.MLE.joint is similar to mKrig.MLE, except it uses the optim function to optimize over the specified covariance parameters and lambda jointly rather than optimizing on a grid. Unlike mKrig.MLE, it returns an mKrig object.

Value

mKrig.MLE returns a list with the components:

`summary`	A matrix giving the results for evaluating the likelihood for each covariance model.
`par.grid`	The par.grid argument used.
`cov.args.MLE`	The list of covariance arguments (except for lambda) that have the largest likelihood over the list covariance models. To fit the surface at the largest likelihood among those tried `do.call( "mKrig", c(obj$mKrig.args, obj$cov.args.MLE,list(lambda=obj$lambda.opt)) )` where `obj` is the list returned by this function.
`call`	The calling arguments to this function.

mKrig.MLE.joint returns an mKrig object with the best computed log-likelihood computed in the maximization process with the addition of the summary table for the mKrig object log-likelihood and:

lnLike.eval

A table containing information on all likelihood evaluations performed in the maximization process.

Author(s)

Douglas W. Nychka, John Paige

References

https://github.com/NCAR/Fields

Examples

# some synthetic data
  N<- 100
  set.seed(123)
  x<- matrix(runif(2*N), N,2)
  theta<- .2
  Sigma<-  Matern( rdist(x,x)/theta , smoothness=1.0)
  Sigma.5<- chol( Sigma)
  sigma<- .1
  M<-5 #  Five (5) independent spatial data sets
  F.true<- t( Sigma.5)%*% matrix( rnorm(N*M), N,M)
  Y<-  F.true +  sigma* matrix( rnorm(N*M), N,M)
# find MLE for lambda with range and smoothness fixed in Matern for first
# data set
  obj<- mKrig.MLE( x,Y[,1], Covariance="Matern", theta=.2, smoothness=1.0)
  obj$summary # take a look
  fit<- mKrig( x,Y[,1], Covariance="Matern", theta=.2,
                                   smoothness=1.0, lambda= obj$lambda.best) 
#
# search over the range parameter and use all 5 replications for combined
# likelihood
## Not run: 
  par.grid<- list( theta= seq(.1,.25,,6))
# default starting value for lambda is .02 subsequent ones use previous optimum.
  obj<- mKrig.MLE( x,Y, Covariance="Matern",lambda=c(.02,rep(NA,4)),
                                  smoothness=1.0, par.grid=par.grid)

## End(Not run)

#perform joint likelihood maximization over lambda and theta. 
#optim can get a bad answer with poor initial guesses.
set.seed(123)
obj<- mKrig.MLE.joint(x,Y[,1], 
                      cov.args=list(Covariance="Matern", smoothness=1.0), 
                      cov.params.guess=list(theta=.2), lambda.guess=.1)

#look at lnLik evaluations
obj$lnLik.eval

## Not run: 
#perform joint likelihood maximization over lambda, theta, and smoothness.  
#optim can get a bad answer with poor initial guesses.
set.seed(123)
obj<- mKrig.MLE.joint(x,Y[,1], 
                      cov.args=list(Covariance="Matern"), 
                      cov.params.guess=list(theta=.2, smoothness=1), lambda.guess=.1)

#look at lnLik evaluations
obj$lnLik.eval

#generate surface plot of results of joint likelihood maximization
#NOTE: mKrig.MLE.joint returns mKrig object while mKrig.MLE doesn't, 
#so this won't work for mKrig.MLE.
surface(obj)

## End(Not run)

fields

Tools for Spatial Data

v11.6

GPL (>= 2)

Authors

Douglas Nychka [aut, cre], Reinhard Furrer [aut], John Paige [aut], Stephan Sain [aut], Florian Gerber [aut], Matthew Iverson [aut], University Corporation for Atmospheric Research [cph]

Initial release

2020-10-06