R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Simple function to search over covariance parameters for...

LKrig.MLE

R Documentation

Simple function to search over covariance parameters for Lattice Krig.

Description

Given a list of different covariance parameters for the Lattice Krig covariance model this function computes the likelihood or a profiled version (over lambda) and approximates a generalized cross-validation function at each of the parameter settings. This is an experimental function that has been productively used with a Latin hypercube design package to efficiently search through the LatticeKrig covariance parameter space.

Usage

LKrigFindLambda(x, y, ..., LKinfo,
                            use.cholesky=NULL,                          
                            lambda.profile=TRUE,
                            lowerBoundLogLambda=-16,tol = 0.005,
                            verbose=FALSE)

LKrig.MLE( x,y,..., LKinfo, use.cholesky = NULL,
                    par.grid=NULL,
                    lambda.profile=TRUE,
                    verbose=FALSE,
                    lowerBoundLogLambda = -16,
                    nTasks = 1, taskID = 1,
                    tol = 0.005)
LKrig.make.par.grid(par.grid=NULL, LKinfo = NULL)

Arguments

`x`	The spatial locations.
`y`	The observations.
`par.grid`	A list with components llambda, alpha, a.wght giving the different sets of parameters to evaluate. If M is the number of parameter setting to evaluate llambda is a vector length M and alpha and a.wght are matrices with M rows and nlevel columns. Thus, the kth trial has parameters par.grid$llambda[k], par.grid$alpha[k,] and par.grid$a.wght[k,]. Currently this function does not support passing a non-stationary spatial parameterization for alpha. The LKinfo object details the other parts of the covariance specification (e.g. number of levels, grid sizes) that do not change. Note that par.grid assumes ln lambda not lambda. See details below for some other features of the par.grid arguments.
`lambda.profile`	A logical value controlling whether the likelihood is maximized over lambda. For LKrigFindLambda if TRUE the likelihood is maximized over lambda at the covariance values in LKinfo and if FALSE the likelihood is just evaluated at LKinfo including the lambda value in this list. For LKrig.MLE if TRUE for each set of parameters in par.grid the value of lambda is found that maximizes the likelihood. In this case the llambda value is the starting value for the optimizer. If llammbda[k] is NA then the lambda value found from the k-1 maximization is used as a starting value for the k step. (In the source code this is `llambda.opt`.) Of course this only makes sense if the other parameters are ordered so that the results for k-1 make sense as a lambda starting value for k. If FALSE the likelihood is evaluated for the covariance parameters at the kth positions in the par.grid list including lambda.
`LKinfo`	An LKinfo object that specifies the LatticeKrig covariance. Usually this is obtained by a call to `LKrigSetup` or as the component LKinfo from the LKrig object. The search sequentially replaces the alpha and a.wght arguments in this list by the values in par.grid but leaves everything else the same. If par.grid is not passed the parameter values in LKinfo are used to evaluate the likelihood. This option is most useful if one has fixed values of alpha and a.wght and the goal is to maximize the likelihood over lambda.
`lowerBoundLogLambda`	Lower limit for lambda in searching for MLE.
`nTasks`	If using Rmpi the number of slaves available.
`tol`	Tolerance on log likelihood use to determine convergence.
`taskID`	If using Rmpi the slave id.
`verbose`	If TRUE prints out intermediate results.
`use.cholesky`	If not NULL then this object is used as the symbolic cholesky decomposition of the covariance matrix for computing the likelihood.
`...`	Any arguments to be passed to LKrig. E.g. `x`, `y`, `Z` a covariate or `weights`.

Details

LKrigFindLambda: Uses a simple one dimensional optimizer optimize. To maximize the log likelihood for log lambda over the range: llambda.start + [-8,5]. This function is used to determine lambda in LatticeKrig.

LKrig.MLE: This is a simple wrapper function to accomplish repeated calls to the LKrig function to evaluate the profile likelihood and/or to optimize the likelihood over the lambda parameters. The main point is that maximization over the lambda parameter (or equivalently for sigma and rho) is the most important and should be done before considering variation of other parameters. If lambda is specified then one has closed form expressions for sigma, rho that can then be substituted back into the log full likelihood. This operation that is the default throughout LatticeKrig (and fields) can concentrate the likelihood on a reduced set of parameters. The further refinement when lambda.profile==TRUE is to maximize the concentrated likelihood over lambda and report this result. This will be a profile likelihood over the remaining parameters of the covariance.

The covariance/model parameters are alpha, a.wght, and log lambda and are separate matrix or vector components of the par.grid list. The cleanest version of this function would just require the par.grid list, however, to be easier to use there are several options to give partial information and let the function itself create the master parameter list. For example, just a search over lambda should be easy and not require creating par.grid outside the function. To follow this option one can just give an LKinfo object. The value for the lambda component in this object will be the starting value with the default starting value being lambda =1.0.

In the second example below most of the coding is getting the grid of parameters to search in the right form. It is useful to normalize the alpha parameters to sum to one so that the marginal variance of the process is only parameterized by rho. To make this easy to implement there is the option to specify the alpha parameters in the form of a mixture model so that the components are positive and add to one (the gamma variable below). If a component gamma is passed as a component of par.grid then this is assumed to be in the mixture model form and the alpha weights are computed from this. Note that gamma will be a matrix with (nlevel - 1) columns while alpha has nlevel columns.

For those readers that use which.max these functions are natural extensions and are handy for looking at interpolated surfaces of the likelihood function.

which.max.matrix: Finds the maximum value in a matrix and returns the row/column index.

which.max.image Finds the maximum value in an image matrix and returns the index and the corresponding grid values.

LKrig.make.par.grid: This is usually used as an internal function that converts the list of parameters in par.grid and the LKinfo object into an more complex data structure used by LKrig.MLE. Its returned value is a "list of lists" to make the search over different parameters combinations simple.

Value

LKrigFindLambda

`summary`	Giving information on the optimization over lambda.
`LKinfo`	Covariance information object.
`llambda.start, lambda.MLE`	Initial and final values for lambda.
`lnLike.eval`	Matrix with all values of log likelihood that were evaluated
`call`	Calling arguments.
`Mc`	Cholesky decomposition.

LKrig.MLE

`summary`	A matrix with columns: effective degrees of freedom, ln Profile likelihood, Generalized cross-validation function, MLE sigma, MLE rho, full likelihood and number of parameter evaluations. The rows correspond to the different parameters in the rows of the `par.grid` components.
`par.grid`	List of parameters used in search. Some parameters might be filled in from the initial par.grid list passed and also from `LKinfo`.
`LKinfo`	LKinfo list that was either passed or created.
`index.MLE`	Index for case that has largest Likelihood value.
`index.GCV`	Index for case that has largest GCV value.
`LKinfo.MLE`	LKinfo list at the parameters with largest profile likelihood.
`lambda.MLE`	Value of lambda from grid with largest profile likelihood.
`call`	Calling sequence for this function.

which.max.matrix Returns a 2 column matrix with row and column index of maximum.

which.max.image For an object in image format returns components x,y,z giving the location and the maximum value for the image. Also included is the component ind that is the row and column indices for the maximum in the image matrix.

LKrig.make.par.grid Returns a list with components, alpha, a.wght. Each component is a list where each component of the list is a separate set of parameters. This more general format is useful for the nonstationary case when the parameters alpha might be a list of nlevel matrices.

Author(s)

Douglas Nychka

Examples

# 
# fitting summer precip for  sub region of North America
  data(NorthAmericanRainfall)
# rename for less typing
  x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
# total precip in 1/10 mm for JJA 
 y<- log10(NorthAmericanRainfall$precip)
# cut down the size of this data set so examples run quickly
# examples also work with  the full data set. Also try NC= 100 for a
# nontrivial model.
  ind<- x[,1] > -90 & x[,2] < 35 #
  x<- x[ind,]
  y<- y[ind]

# This is a single level smoother
  LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=4.05, alpha=1.0)
  lambdaFit<- LKrigFindLambda( x,y,LKinfo=LKinfo)
  lambdaFit$summary

  NG<-3
 #NOTE: make this larger ( e.g. 15) for a better (real!) grid search on log lambda
 # values
  par.grid<- list( a.wght= rep( 4.05,NG),alpha= rep(1, NG),
                      llambda=  seq(-8,-2,,NG))
  LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=5, alpha=1.0)
  lambda.search.results<-LKrig.MLE( x,y,LKinfo=LKinfo,
                                    par.grid=par.grid,
                                    lambda.profile=FALSE)
  lambda.search.results$summary
# profile likelihood
  plot( lambda.search.results$summary[,1:2], 
         xlab="effective degrees freedom",
         ylab="ln profile likelihood")
# fit at largest likelihood value:
    lambda.MLE.fit<- LKrig( x,y,
                    LKinfo=lambda.search.results$LKinfo.MLE)
## Not run:                     
# optimizing  Profile likelihood over lambda using optim
# consider 3 values for a.wght (range parameter)
# in this case the log lambdas passed are the starting values for optim.
  NG<-3
  par.grid<- list( a.wght= c( 4.05,4.1,5),alpha= rep(1, NG),
                      llambda= c(-5,NA,NA))
# NOTE: NAs in llambda mean use the previous MLE for llambda as the
# current starting value. 
  LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=5, alpha=1.0) 
  lambda.search.results<-LKrig.MLE(
                              x,y,LKinfo=LKinfo, par.grid=par.grid,
                              lambda.profile=TRUE, verbose=TRUE)
  print(lambda.search.results$summary)
# note first result a.wght = 4.05 is the optimized result for the grid
# search given above.

## End(Not run)
########################################################################    
# search over two multi-resolution levels varying the  levels of alpha's
########################################################################
## Not run: 
# NOTE: search ranges found largely by trial and error to make this
# example work also the grid is quite coarse ( and NC is small) to
# be quick as a help file example
  data(NorthAmericanRainfall)
# rename for less typing
  x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
# total precip in 1/10 mm for JJA 
 y<- log10(NorthAmericanRainfall$precip)
# cut down the size of this data set so examples run quickly
# examples also work with  the full data set. Also try NC= 100 for a
# nontrivial model.
  ind<- x[,1] > -90 & x[,2] < 35 #
  x<- x[ind,]
  y<- y[ind]
  
  Ndes<- 10  # NOTE: this is set very small just to make example run fast
  set.seed(124)
  par.grid<- list()
# create grid of alphas to sum to 1 use a mixture model parametrization
#  alpha1 = (1/(1 + exp(gamma1)) ,
# alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
# 
  par.grid$gamma<- cbind(runif( Ndes, -3,2), runif( Ndes, -3,2))
  par.grid$a.wght<- matrix( 4.5, nrow=Ndes, ncol=3)
# log lambda grid search values
  par.grid$llambda<- runif( Ndes,-5,-3)  
   LKinfo1<- LKrigSetup( x, NC=5, nlevel=3, a.wght=5, alpha=c(1.0,.5,.25))
# NOTE: a.wght in call is not used. Also a better search is to profile over
#  llambda

 alpha.search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
                                    lambda.profile=FALSE)

########################################################################
# Viewing the search results
########################################################################

# this scatterplot is good for a quick look because  effective degrees
# of freedom is a useful summary of fit. 
  plot( alpha.search.results$summary[,1:2], 
         xlab="effective degrees freedom",
         ylab="ln profile likelihood")
#

## End(Not run)

## Not run: 
# a two level model search 
# with profiling over lambda.
#  This takes a few minutes
  Ndes<- 40 
  nlevel<-2 
  par.grid<- list()
## create grid of alphas to sum to 1 use a mixture model parametrization:
#    alpha1 = (1/(1 + exp(gamma1)) ,
#   alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
  set.seed(123)
  par.grid$gamma<- runif( Ndes,-3,4)
## values for range (a.wght)
  a.wght<-  4 + 1/exp(seq( 0,4,,Ndes))
  par.grid$a.wght<- cbind( a.wght, a.wght)
# log lambda grid search values (these are the starting values)
  par.grid$llambda<- rep(-4, Ndes)

  LKinfo1<- LKrigSetup( x, NC=15, nlevel=nlevel, 
                          a.wght=5, alpha=rep( NA,2) ) 
##
## the search over the parameter list in par.grid  maximizing over lambda 
  search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
                                 lambda.profile=TRUE)
# plotting results
set.panel(1,2)
 plot( search.results$summary[,1:2], 
         xlab="effective degrees freedom",
         ylab="ln profile likelihood")
 xtemp<- matrix(NA, ncol=2, nrow=Ndes)
 for( k in 1:Ndes){
   xtemp[k,] <- c( (search.results$par.grid$alpha[[k]])[1],
                  (search.results$par.grid$a.wght[[k]])[1] )
}
 quilt.plot( xtemp,search.results$summary[,2])
#  fit using Tps
 tps.out<- Tps(  xtemp,search.results$summary[,2], lambda=0)
 contour( predictSurface(tps.out), lwd=3,add=TRUE)
 set.panel()

## End(Not run)
## Not run: 
# searching over nu 
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
good<- !is.na(y)
y<- y[good]
x<- x[good,]
par.grid<- expand.grid( nu=seq(.5,1.5,,4), a.wght=c(4.01,4.1, 4.2,4.5,5))
par.grid$llambda<- rep( NA, length(par.grid$nu))
LKinfo<- LKrigSetup(x,  nlevel=3,NC=5)
out<- LKrig.MLE( x,y, LKinfo=LKinfo, par.grid=par.grid, verbose=TRUE)

## End(Not run)
## Not run: 
# an MLE fit taking advantage of replicated fields
set.seed(222)
N<- 500
M<-100
sigma<- .1
x<- matrix( runif(N*2), N,2)
LKinfo<- LKrigSetup( x, NC=5, nlevel=4,nu=1.0, a.wght=4.2 )
# the replicate fields
y<- LKrig.sim(x,LKinfo=LKinfo, M=M ) + sigma*matrix( rnorm(N*M), N,M)
# with correct lambda
obj<- LKrig( x,y, LKinfo=LKinfo, lambda=.1)
print( obj$sigma.MLE)
print( obj$rho.MLE)

## End(Not run)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(LatticeKrig)
Loading required package: spam
Loading required package: grid
Spam version 1.3-0 (2015-10-24) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: 'spam'

The following objects are masked from 'package:base':

    backsolve, forwardsolve

Loading required package: fields
Loading required package: maps

 # maps v3.1: updated 'world': all lakes moved to separate new #
 # 'lakes' database. Type '?world' or 'news(package="maps")'.  #


> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/LatticeKrig/LKrig.MLE.Rd_%03d_medium.png", width=480, height=480)
> ### Name: LKrig.MLE
> ### Title: Simple function to search over covariance parameters for Lattice
> ###   Krig.
> ### Aliases: LKrig.MLE LKrigFindLambda LKrig.make.par.grid
> ### Keywords: spatial
> 
> ### ** Examples
> 
> # 
> # fitting summer precip for  sub region of North America
>   data(NorthAmericanRainfall)
> # rename for less typing
>   x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
> # total precip in 1/10 mm for JJA 
>  y<- log10(NorthAmericanRainfall$precip)
> # cut down the size of this data set so examples run quickly
> # examples also work with  the full data set. Also try NC= 100 for a
> # nontrivial model.
>   ind<- x[,1] > -90 & x[,2] < 35 #
>   x<- x[ind,]
>   y<- y[ind]
> 
> # This is a single level smoother
>   LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=4.05, alpha=1.0)
>   lambdaFit<- LKrigFindLambda( x,y,LKinfo=LKinfo)
>   lambdaFit$summary
        EffDf    lnProfLike           GCV     sigma.MLE       rho.MLE 
 3.461127e+01  1.994752e+02  8.941907e-04  2.501174e-02  8.082863e-02 
  llambda.opt        lnLike  counts value          grad 
-4.861396e+00            NA  1.300000e+01            NA 
> 
>   NG<-3
>  #NOTE: make this larger ( e.g. 15) for a better (real!) grid search on log lambda
>  # values
>   par.grid<- list( a.wght= rep( 4.05,NG),alpha= rep(1, NG),
+                       llambda=  seq(-8,-2,,NG))
>   LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=5, alpha=1.0)
>   lambda.search.results<-LKrig.MLE( x,y,LKinfo=LKinfo,
+                                     par.grid=par.grid,
+                                     lambda.profile=FALSE)
>   lambda.search.results$summary
        EffDf lnProfLike          GCV  sigma.MLE    rho.MLE llambda.MLE lnLike
[1,] 56.57676   168.5842 0.0009847125 0.01725275 0.88730410          -8     NA
[2,] 38.16865   199.4183 0.0009475525 0.02448106 0.08894732          -5     NA
[3,] 17.87395   183.9233 0.0015738388 0.03958272 0.01157711          -2     NA
     counts value grad
[1,]           NA   NA
[2,]           NA   NA
[3,]           NA   NA
> # profile likelihood
>   plot( lambda.search.results$summary[,1:2], 
+          xlab="effective degrees freedom",
+          ylab="ln profile likelihood")
> # fit at largest likelihood value:
>     lambda.MLE.fit<- LKrig( x,y,
+                     LKinfo=lambda.search.results$LKinfo.MLE)
> ## Not run: 
> ##D                     
> ##D # optimizing  Profile likelihood over lambda using optim
> ##D # consider 3 values for a.wght (range parameter)
> ##D # in this case the log lambdas passed are the starting values for optim.
> ##D   NG<-3
> ##D   par.grid<- list( a.wght= c( 4.05,4.1,5),alpha= rep(1, NG),
> ##D                       llambda= c(-5,NA,NA))
> ##D # NOTE: NAs in llambda mean use the previous MLE for llambda as the
> ##D # current starting value. 
> ##D   LKinfo<- LKrigSetup(x,NC=12,nlevel=1, a.wght=5, alpha=1.0) 
> ##D   lambda.search.results<-LKrig.MLE(
> ##D                               x,y,LKinfo=LKinfo, par.grid=par.grid,
> ##D                               lambda.profile=TRUE, verbose=TRUE)
> ##D   print(lambda.search.results$summary)
> ##D # note first result a.wght = 4.05 is the optimized result for the grid
> ##D # search given above.
> ## End(Not run)
> ########################################################################    
> # search over two multi-resolution levels varying the  levels of alpha's
> ########################################################################
> ## Not run: 
> ##D # NOTE: search ranges found largely by trial and error to make this
> ##D # example work also the grid is quite coarse ( and NC is small) to
> ##D # be quick as a help file example
> ##D   data(NorthAmericanRainfall)
> ##D # rename for less typing
> ##D   x<- cbind( NorthAmericanRainfall$longitude, NorthAmericanRainfall$latitude)
> ##D # total precip in 1/10 mm for JJA 
> ##D  y<- log10(NorthAmericanRainfall$precip)
> ##D # cut down the size of this data set so examples run quickly
> ##D # examples also work with  the full data set. Also try NC= 100 for a
> ##D # nontrivial model.
> ##D   ind<- x[,1] > -90 & x[,2] < 35 #
> ##D   x<- x[ind,]
> ##D   y<- y[ind]
> ##D   
> ##D   Ndes<- 10  # NOTE: this is set very small just to make example run fast
> ##D   set.seed(124)
> ##D   par.grid<- list()
> ##D # create grid of alphas to sum to 1 use a mixture model parametrization
> ##D #  alpha1 = (1/(1 + exp(gamma1)) ,
> ##D # alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
> ##D # 
> ##D   par.grid$gamma<- cbind(runif( Ndes, -3,2), runif( Ndes, -3,2))
> ##D   par.grid$a.wght<- matrix( 4.5, nrow=Ndes, ncol=3)
> ##D # log lambda grid search values
> ##D   par.grid$llambda<- runif( Ndes,-5,-3)  
> ##D    LKinfo1<- LKrigSetup( x, NC=5, nlevel=3, a.wght=5, alpha=c(1.0,.5,.25))
> ##D # NOTE: a.wght in call is not used. Also a better search is to profile over
> ##D #  llambda
> ##D 
> ##D  alpha.search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
> ##D                                     lambda.profile=FALSE)
> ##D 
> ##D ########################################################################
> ##D # Viewing the search results
> ##D ########################################################################
> ##D 
> ##D # this scatterplot is good for a quick look because  effective degrees
> ##D # of freedom is a useful summary of fit. 
> ##D   plot( alpha.search.results$summary[,1:2], 
> ##D          xlab="effective degrees freedom",
> ##D          ylab="ln profile likelihood")
> ##D #
> ## End(Not run)
> 
> ## Not run: 
> ##D # a two level model search 
> ##D # with profiling over lambda.
> ##D #  This takes a few minutes
> ##D   Ndes<- 40 
> ##D   nlevel<-2 
> ##D   par.grid<- list()
> ##D ## create grid of alphas to sum to 1 use a mixture model parametrization:
> ##D #    alpha1 = (1/(1 + exp(gamma1)) ,
> ##D #   alpha2 = exp( gamma1) / ( 1 + exp( gamma1))
> ##D   set.seed(123)
> ##D   par.grid$gamma<- runif( Ndes,-3,4)
> ##D ## values for range (a.wght)
> ##D   a.wght<-  4 + 1/exp(seq( 0,4,,Ndes))
> ##D   par.grid$a.wght<- cbind( a.wght, a.wght)
> ##D # log lambda grid search values (these are the starting values)
> ##D   par.grid$llambda<- rep(-4, Ndes)
> ##D 
> ##D   LKinfo1<- LKrigSetup( x, NC=15, nlevel=nlevel, 
> ##D                           a.wght=5, alpha=rep( NA,2) ) 
> ##D ##
> ##D ## the search over the parameter list in par.grid  maximizing over lambda 
> ##D   search.results<- LKrig.MLE( x,y,LKinfo=LKinfo1, par.grid=par.grid,
> ##D                                  lambda.profile=TRUE)
> ##D # plotting results
> ##D set.panel(1,2)
> ##D  plot( search.results$summary[,1:2], 
> ##D          xlab="effective degrees freedom",
> ##D          ylab="ln profile likelihood")
> ##D  xtemp<- matrix(NA, ncol=2, nrow=Ndes)
> ##D  for( k in 1:Ndes){
> ##D    xtemp[k,] <- c( (search.results$par.grid$alpha[[k]])[1],
> ##D                   (search.results$par.grid$a.wght[[k]])[1] )
> ##D }
> ##D  quilt.plot( xtemp,search.results$summary[,2])
> ##D #  fit using Tps
> ##D  tps.out<- Tps(  xtemp,search.results$summary[,2], lambda=0)
> ##D  contour( predictSurface(tps.out), lwd=3,add=TRUE)
> ##D  set.panel()
> ## End(Not run)
> ## Not run: 
> ##D # searching over nu 
> ##D data(ozone2)
> ##D x<- ozone2$lon.lat
> ##D y<- ozone2$y[16,]
> ##D good<- !is.na(y)
> ##D y<- y[good]
> ##D x<- x[good,]
> ##D par.grid<- expand.grid( nu=seq(.5,1.5,,4), a.wght=c(4.01,4.1, 4.2,4.5,5))
> ##D par.grid$llambda<- rep( NA, length(par.grid$nu))
> ##D LKinfo<- LKrigSetup(x,  nlevel=3,NC=5)
> ##D out<- LKrig.MLE( x,y, LKinfo=LKinfo, par.grid=par.grid, verbose=TRUE)
> ## End(Not run)
> ## Not run: 
> ##D # an MLE fit taking advantage of replicated fields
> ##D set.seed(222)
> ##D N<- 500
> ##D M<-100
> ##D sigma<- .1
> ##D x<- matrix( runif(N*2), N,2)
> ##D LKinfo<- LKrigSetup( x, NC=5, nlevel=4,nu=1.0, a.wght=4.2 )
> ##D # the replicate fields
> ##D y<- LKrig.sim(x,LKinfo=LKinfo, M=M ) + sigma*matrix( rnorm(N*M), N,M)
> ##D # with correct lambda
> ##D obj<- LKrig( x,y, LKinfo=LKinfo, lambda=.1)
> ##D print( obj$sigma.MLE)
> ##D print( obj$rho.MLE)
> ## End(Not run)
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>