fields (version 8.3-6)

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

This function is designed to explore the likelihood surface for different covariance parameters with the option of maximizing over sigma and rho.

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)

Value

mKrig.MLE returns a list with the components:summaryA matrix giving the results for evaluating the likelihood for each covariance model.par.gridThe par.grid argument used.cov.args.MLEThe 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.

callThe 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.evalA table containing information on all likelihood evaluations performed in the maximization process.

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

References

http://cran.r-project.org/web/packages/fields/fields.pdf http://www.image.ucar.edu/~nychka/Fields/

See Also

mKrig Krig stationary.cov optim

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
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)

#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

#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)