mKrigMLE: 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 function are designed to explore the likelihood surface for different covariance parameters with the option of maximizing over sigma and rho. They used the mKrig base are designed for computational efficiency.

Usage

mKrigMLEGrid(x, y, weights = rep(1, nrow(x)), Z = NULL, mKrig.args = NULL, 
    cov.fun = "stationary.cov", cov.args = NULL, na.rm = TRUE, 
    par.grid = NULL, lambda = NULL, lambda.profile = TRUE,
    relative.tolerance = 1e-04, 
    REML = FALSE, verbose = FALSE)
mKrigMLEJoint(x, y, weights = rep(1, nrow(x)), Z = NULL, mKrig.args
                 = NULL, na.rm = TRUE, cov.fun = "stationary.cov",
                 cov.args = NULL, lambda.start = 0.5, cov.params.start
                 = NULL, optim.args = NULL, abstol = 1e-04,
                 parTransform = NULL, REML = FALSE, verbose = FALSE)

fastTpsMLE(x, y, weights = rep(1, nrow(x)), Z = NULL, ...,
                 par.grid=NULL, theta, lambda = NULL, lambda.profile = TRUE,
                 verbose = FALSE, relative.tolerance = 1e-04)
mKrigJointTemp.fn(parameters, mKrig.args, cov.args, parTransform,
           parNames, REML = FALSE, capture.env)

Arguments

abstol

Absolute convergence tolerance used in optim.

capture.env

For the ML obective function the frame to save the results of the evaluation. This should be the environment of the function calling optim.

cov.fun

The name, a text string, of the covariance function.

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.params.start

A list of initial starts 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.

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.start

The initial guess for lambda in the joint log-likelihood maximization process

lambda.profile

If TRUE maximize likelihood over lambda.

mKrig.args

A list of additional parameters to supply to the base mKrig function that are distinct from the covariance model. For example mKrig.args= list( m=1 ) will set the polynomial to be just a constant term (degree = m -1 = 0).

na.rm

Remove NAs from data.

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.start)+1), abstol=1e-04, maxit=20)) Note that the first parameter is lambda and the others are the covariance parameters in the order they are given in cov.params.start. Also note that the optimization is performed on a transformed scale (based on the function parTransform ), and this should be taken into consideration when passing arguments to optim.

parameters

The parameter values for evaluate the likelihood.

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 ….

parNames

Names of the parameters to optimize over.

parTransform

A function that maps the parameters to a scale for optimization or effects the inverse map from the transformed scale into the original values. See below for more details.

relative.tolerance

Tolerance used to declare convergence when maximizing likelihood over lambda.

REML

Currently using REML is not implemented.

theta

Range parameter for compact Wendland covariance. (see fastTps)

verbose

If TRUE print out interesting intermediate results.

weights

Precision ( 1/variance) of each observation

Matrix of unique spatial locations (or in print or surface the returned mKrig object.)

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.

Linear covariates to be included in fixed part of the model that are distinct from the default low order polynomial in x

…

Other arguments to pass to the mKrig function.

Value

mKrigMLEGrid 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. NOTE: 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.

mKrigMLEJoint returns a list with components:

summary

A vector giving the MLEs and the log likelihood at the maximum

lnLike.eval

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

optimResults

The list returned from the optim function.

par.MLE

The maximum likelihood estimates.

parTransform

The transformation of the parameters used in the optimziation.

Details

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 fastTpsMLE is a convenient variant of this more general version to use directly with fastTps, and mKrigMLEJoint is similar to mKrigMLEGrid, except it uses the optim function to optimize over the specified covariance parameters and lambda jointly rather than optimizing on a grid. Unlike mKrigMLEJoint, it returns an mKrig object.

For mKrigMLEJoint the default transformation of the parameters is set up for a log/exp transformation:

 parTransform <- function(ptemp, inv = FALSE) {
            if (!inv) {
                log(ptemp)
            }
            else {
                exp(ptemp)
            }
        }

References

https://github.com/NCAR/Fields

Examples

Run this code

# NOT RUN {
#perform joint likelihood maximization over lambda and theta. 
#optim can get a bad answer with poor initial starts.
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
obj<- mKrigMLEJoint(x,y, 
                      cov.args=list(Covariance="Matern", smoothness=1.0), 
                      cov.params.start=list(theta=.2), lambda.start=.1)
#                      
# check  lnLikeihood evaluations that were culled from optim
# these are in obj$lnLike.eval
# funny ranges are set to avoid  very low likelihood values

quilt.plot( log10(cbind(obj$lnLike.eval[,1:2])), obj$lnLike.eval[,5],
             xlim=c(-1.2,-.40), ylim=c( -1,1), zlim=c( -625, -610))
             points( log10(obj$pars.MLE[1]), log10(obj$pars.MLE[2]),
	     pch=16, col="grey" )

# some synthetic data with replicates
  N<- 50
  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
  #  250 independent spatial data sets but a common covariance function 
  #    -- there is little overhead in
  #        MLE across independent realizations and a good test of code validity.
  M<-250
  #F.true<- t( Sigma.5)%*% matrix( rnorm(N*M), N,M)  
  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 grid of ranges 
# and smoothness fixed in Matern                     
 par.grid<- list( theta= seq( .1,.35,,8))
  obj1b<- mKrigMLEGrid( x,Y,
     cov.args = list(Covariance="Matern", smoothness=1.0), 
        par.grid = par.grid
                    )
  obj$summary # take a look
# profile over theta
  plot( par.grid$theta, obj1b$summary[,"lnProfileLike.FULL"],
    type="b", log="x")
  
  
# }
# NOT RUN {
# m=0 is a simple switch to indicate _no_ fixed spatial drift
# (the default and highly recommended  is linear drift, m=2). 
# this results in MLEs that are less biased -- in fact it nails it !
  obj1a<- mKrigMLEJoint(x,Y, 
                    cov.args=list(Covariance="Matern", smoothness=1.0), 
                    cov.params.start=list(theta=.5), lambda.start=.5,
                     mKrig.args= list( m=0))
 
 test.for.zero( obj1a$summary["sigmaMLE"], sigma, tol=.0075)
 test.for.zero( obj1a$summary["theta"], theta, tol=.05)
# }
# NOT RUN {
# }
# NOT RUN {
#perform joint likelihood maximization over lambda, theta, and smoothness.  
#note: finding smoothness is not robust optim 
#      can get a bad answer with poor initial guesses.
obj2<- mKrigMLEJoint(x,Y, 
                      cov.args=list(Covariance="Matern"), 
                      cov.params.start=list(theta=.18, smoothness=1.1),
                       lambda.start=.08)

#look at lnLikelihood  evaluations
obj2$summary
#compare to REML
obj3<- mKrigMLEJoint(x,Y, 
                      cov.args=list(Covariance="Matern"), 
                      cov.params.start=list(theta=.18, smoothness=1.1),
                        lambda.start=.08
                       , REML=TRUE)
# }
# NOT RUN {
#look at lnLikelihood  evaluations
obj3$summary
# check convergence of MLE to true fit with no fixed part
# 
obj4<- mKrigMLEJoint(x,Y, 
                      mKrig.args= list( m=0),
                      cov.args=list(Covariance="Matern", smoothness=1), 
                      cov.params.start=list(theta=.2),
                       lambda.start=.1, REML=TRUE)
#look at lnLikelihood  evaluations
obj4$summary
# nails it!
# }

Run the code above in your browser using DataCamp Workspace