Simple function to search over covariance parameters for Lattice Krig.

Share:

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
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

See Also

LKrig LatticeKrig

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
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
# 
# 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)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.