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
1 2 3 4 5 6 7 8 9 10 11 12 13  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 = 1e04, ...)
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 = 1e04)

Arguments
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. 
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. 
weights 
Precision ( 1/variance) of each observation 
Z 
Linear covariates to be included in fixed part of the
model that are distinct from the default low order
polynomial in 
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 .... 
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. 
lambda 
If 
lambda.guess 
The initial guess for lambda in the joint loglikelihood maximization process 
lambda.profile 
If 
optim.args 
Additional arguments that would also be included in calls
to the optim function in joint likelihood maximization. If

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
loglikelihood is calculated depending on 
... 
Additional arguments that would also be included in a call to

verbose 
If 
relative.tolerance 
Tolerance used to declare convergence when maximizing likelihood over lambda. 
theta 
Range parameter for compact Wendland covariance. (see fastTps) 
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.
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

call 
The calling arguments to this function. 
mKrig.MLE.joint
returns an mKrig object with the best
computed loglikelihood computed in the maximization process
with the addition of the summary table for the mKrig object
loglikelihood and:
lnLike.eval 
A table containing information on all likelihood evaluations performed in the maximization process. 
Author(s)
Douglas W. Nychka, John Paige
References
http://cran.rproject.org/web/packages/fields/fields.pdf http://www.image.ucar.edu/~nychka/Fields/
See Also
mKrig
Krig
stationary.cov
optim
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55  # 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)
