LKrig.basis: High level functions for generating and simulating a basis function,Gaussian Process.

Description

These functions support the LKrig function. Their main function is to create and evaluate radial basis functions of varying support on a nested set of regular grids. This series of grids forms a multi-resolution basis. The Gaussian process model is an expansion in these basis functions where the basis coefficients follow a Markov random field model for each resolution level. This family of functions generate the basis using sparse matrices, evaluate the covariance function of the process, and also simulate realizations of the process. LKrig.cov.plot is a useful function to get a quick plot of the covariance function implied by a LatticeKrig specification.

Usage

#
LKrig.cov(x1, x2 = NULL, LKinfo, C = NA, marginal = FALSE,
                 aRange = NA)
LKrig.cov.plot( LKinfo, NP=200, center = NULL, xlim = NULL, ylim = NULL)
LKrigCovWeightedObs(x1, wX, LKinfo) 
LKrig.basis(x1, LKinfo, Level = NULL, raw = FALSE, verbose = FALSE)
LKrig.precision(LKinfo, return.B = FALSE, verbose=FALSE)  
LKrig.quadraticform( Q, PHI, choleskyMemory = NULL) 
LKrig.spind2spam(obj, add.zero.rows=TRUE)
LKrigMarginalVariance(x1, LKinfo, verbose = FALSE)
LKFindSigma2VarianceWeights(x1, LKinfo)

Value

LKrig.basis: A matrix with number of rows equal to the rows of x1 and columns equal to the number of basis functions (LKinfo$m). Attached to the matrix is an info attribute that contains the list passed as LKinfo. Usually this value is in spam sparse matrix format.

LKrig.precision: For return.B ==FALSE a sparse, square matrix with dimensions of the number of basis functions. For return.B == TRUE the "B" SAR matrix is returned. This is useful for checking this function.

LKrig.cov: If C=NA a cross covariance matrix with dimensions nrow(x1) and nrow(x2) is used. If C is passed the result of multiplying the cross covariance matrix times C is used.

LKrigMarginalVariance: Gives the marginal variance of the LatticeKrig process at each level and at the locations in x1. Returned value is a matrix with nlevel columns indexing the levels and the number of rows equal to nrow(x1). If varLevels is a row of this matrix then sum( varLevels* LKinfo$alpha) is the marginal variance for the full process when the different levels are weighted by alpha. This is weighted sum is the marginal variance returned by LKrig.cov and marginal=TRUE ( Also assuming that LKinfo$sigma2.object

is NULL, which it usually is.) .

LKFindSigma2VarianceWeights Return either the marginal variance of the process as a scalar or the variance at arbitrary locations x1. This second role is part of a nonstationary specification of the model where the process marginal variance can vary over space.

LKrig.sim: A matrix with dimensions of nrow(x1) by M. Each column are vectors of simulated values at the locations x1.

LKrig.cov.plot: Evaluates the covariance specified in the list LKinfo with respect to the point center along a transects in the x and y directions intersecting this point. Note the rectangular extent of the spatial domain is part of the grid information in LKinfo. Returns components u, d and cov. Each of these are two column matrices with the first column being results in the x direction and second column in the y direction. d are the distances of the points from the center and u are the actual x or y spatial coordinates. cov are the values of the covariance function. If normalize is TRUE these will in fact be the correlation functions. To plot the returned list use


 out<- LKrig.cov.plot(LKinfo)
 matplot( out$d, out$cov, type="l")

LKrig.quadraticform: Returns a vector that is

diag(t(PHI)%*% solve(Q) %*% PHI)) closely related to the marginal variance of the process.

LKrigNormalizeBasis, LKrigNormalizeBasisFast, LKrigNormalizeBasisFFTInterpolate: A vector of variances corresponding to the unnormalized process at the locations.

LKrig.spind2spam: This converts a matrix in spind sparse format to spam sparse format. Although useful for checking, LatticeKrig now uses the the spam function directly to do the conversion. If obj is a list in spind format then the

obj2<- LKrig.spind2spam(obj)

and

obj2<- spam(obj[c("ind","ra")], nrow=obj$da[1], ncol=obj$da[2])

give the same result.

Arguments

add.zero.rows: If TRUE the conversion of the sparse matrix to spam format will have at least one element in each row. If there are no elements explicitly given in obj then an element with value zero is added. This technical detail is needed to accommodate the spam format for sparse matrices.
C: If passed the covariance matrix will be multiplied by this vector or matrix.
center: The point in the spatial domain that is used to evaluate the covariance function. The evaluation is done on x and y transects through the spatial domain intersecting at center and finding the covariance with respect to this point. If NULL defaults to the center of the spatial domain.
choleskyMemory: A list giving the memory requirements for a sparse cholesky decomposition. See chol.spam for details and also LKrig.
LKinfo: A list with components that give the information describing a multi-resolution basis with a Markov random field used for the covariance of the basis coefficients. This list is created in LKrig or by LKrigSetup and returned in the output object. (See section on returned Value below for this list's description.)
Level: If NULL the entire basis is evaluated. If an integer just that level is evaluated.
marginal: If TRUE returns the marginal variance. Currently not implemented!
NP: Number of points to evaluate the covariance function along each transect of the spatial domain.
obj: An object returned by LKrig or a sparse matrix in row/column format passed to LKrig.spind2spam.
PHI: A sparse matrix of basis functions (rows index points for evaluation, columns index basis functions).
raw: Do not normalize the basis functions even if it is indicated in the LKinfo object.
return.B: If TRUE B is returned instead of the precision matrix t(B)%*%B.
Q: A sparse (spam format) precision matrix.
aRange: Currently aRange is a place holder argument for future development where the correlation range can be specified.
wX: The wX matrix constructed in the LKrig function.
x1: A two column matrix of 2-dimension locations to evaluate basis functions or the first set of locations to evaluate the covariance function or the locations for the simulated process. Rows index the different locations: to be precise x1[i,1:2] are the "x" and "y" coordinates for the i th location.
x2: Second set of locations to evaluate covariance function.
xlim: Limits in x coordinate for evaluating the covariance model. Default is the spatial domain.
ylim: Limits in y coordinate for evaluating the covariance model.Default is the spatial domain.
verbose: If TRUE intermediate steps and other debugging information are printed.

Author

Doug Nychka

Details

The basis functions are two-dimensional radial basis functions based on the compactly supported stationary covariance function (Wendland covariance) and centered on regular grid points with the scaling tied to the grid spacing.

For a basis at the coarsest level, the grid centers are generated by expanding the two sequences


seq(grid.info$xmin,grid.info$xmax,grid.info$delta)
seq(grid.info$ymin,grid.info$ymax,grid.info$delta)

into a regular grid of center points. The same spacing delta is used in both directions. The unnormalized basis functions are evaluated at locations x1 by finding the pairwise, radial distances among centers and x1, scaling by grid.info$delta * overlap and then evaluating with the function name passed as BasisFunction. By default this is the 2-d Wendland covariance of order 2. Perhaps the most important point about the LKrig.basis is that it is designed to return a matrix of all basis functions as a sequence of points. There is no need to have a function that evaluates basis functions individually. In R code for a set of locations x1 and a rectangular spatial domain with ranges xmin, xmax, ymin ,ymax:


 centers<- expand.grid(seq(xmin,xmax,delta),
                       seq(ymin,ymax,delta) )
 bigD<- rdist( x1, centers)/(delta*2.5)
 PHI<- Wendland.function( bigD)

Note that there will be nrow(centers) basis functions generated where the precise number depends on the range of the domain and the choice of delta. The basis functions are indexed by the columns in PHI and this is a convention throughout this package. There will be nrow(x1) rows in PHI as each basis function will be evaluated at each 2-d location.

The basis functions are then normalized by scaling the basis functions at each location so that resulting marginal variance of the process is 1. This is done to coax the covariance model closer to a stationary representation. It is easiest to express this normalization by pseudo R code:

If Q is the precision matrix of the basis coefficients then in R/LatticeKrig code:


Omega<-  solve(Q)
process.variance <- diag(PHI%*% Omega %*%t(PHI) )
PHI.normalized <-  diag(1/sqrt(process.variance)) %*% PHI

where Omega is the unnormalized covariance matrix of the basis function coefficients.

Although correct, the code above is not an efficient algorithm to compute the unnormalized process variance. First the normalization can be done level by level rather than dealing with the entire multi-resolution process at once. Also it is important to work with the precision matrix rather than the covariance. The function LKrigNormalizeBasis takes advantage of the sparsity of the precision matrix for the coefficients. LKrigNormalizeBasisFast is a more efficient version for an isotropic type model when a.wght is constant for a given level and takes advantage of the Kronecker structure on the precision matrix at each level. LKrigNormalizeBasisFFTInterpolate is the fastest method, which provides an accurate but approximate method for the normalization, which takes advantage of situations where the number of locations significantly exceeds the number of basis functions.

The precision matrix for the basis coefficients at each resolution has the form t(B)%*% B. These matrices for the individual levels are assembled by LKrig.precision as the block diagonals of a larger precision matrix for the entire vector of coefficients. Note these matrices are created in a sparse format. The specific entries in B, the object created by LKrig.MRF.precision, are a first order Markov random field: without edge adjustments the diagonal elements have the value a.wght and the first order neighbors have the value -1.

Below we give more details on how the weights are determined. Following the notation in Lindgren and Rue a.wght= 4 + k2 with k2 greater than or equal to 0. Some schematics for filling in the B matrix are given below (values are weights for the SAR on the lattice with a period indicating zero weights).


                                                 __________ 
   .   -1     .         |  -e      .            |  a.wght  -e
                        |                       |
  -1  a.wght  -1        | a.wght   -e           |   -e      .
                        |                       |
   .  -1      .         |  -e      .            |    .      .
Interior point        Left edge              Upper left corner

To adjust for edges and corner lattice points we take two strategies. The first is add extra lattice points around the edges so the actual spatial domain has lattice points with the complete number of neighbors. The default for the LKRectangle geometry is to add a buffer of 5 points ( NC.buffer = 5 ) around the edges. In addition the neighbor weights are inflated sum to a fixed value. For example in the case of LKRectangle and following the schematic above we want the neighbor weights to sum to -4. Thus "e" for the middle figure will be set to -4/3 and for the right figure -2. Empirically these adjustments were found to improve how a.wght parameters chosen close to 4 could generate long range spatial correlations in the lattice coefficients. Note that these simple boundary adjustments happen to the buffer points and so are less likely to introduce artifacts into the spatial domain.

Examples

Run this code

# Load ozone data set
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find location that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]
  LKinfo<- LKrigSetup( x,NC=20,nlevel=1, alpha=1, lambda= .3 , a.wght=5)
# BTW lambda is close to MLE 

# What does the  LatticeKrig covariance function look like?
# set up LKinfo object
# NC=10 sets the grid for the first level of basis functions
# NC^2 = 100 grid points in first level if square domain.
# given four levels the number of basis functions
# = 10^2 + 19^2 +37^2 + 73^2 = 5329
# effective range scales as roughly kappa where a.wght =  4 + kappa^2    
# or exponential decreasing marginal variances for the components.
    NC<- 10 
    nlevel <- 4
    a.wght <-  4 + 1/(.5)^2
    alpha<-  1/2^(0:(nlevel-1)) 
    LKinfo2<- LKrigSetup( cbind( c( -1,1), c(-1,1)), NC=NC,
                   nlevel=nlevel, a.wght=a.wght,alpha=alpha)
# evaluate covariance  along the  horizontal line through
# midpoint of region -- (0,0) in this case. 
    look<- LKrig.cov.plot( LKinfo2)
# a plot of the covariance function in x and y with respect to (0,0)
    set.panel(2,1)  
    plot(look$u[,1], look$cov[,1], type="l")
    title("X transect")
    plot(look$u[,2], look$cov[,2], type="l")
    title("Y transect")
    set.panel(1,1)
#
#
if (FALSE) {
# full 2-d view of the covariance (this example follows the code
# in LKrig.cov.plot)
 x2<- cbind( 0,0)
 x1<- make.surface.grid( list(x=seq( -1,1,,40),  y=seq( -1,1,,40)))
 look<- LKrig.cov( x1,x2, LKinfo2)
 contour( as.surface( x1, look))
# Note nearly circular contours.
# of  course  plot(look[,80/2]) should look like plot above.
#
}

if (FALSE) {
#Some correlation functions from different models
set.panel(2,1)
# a selection of ranges:
  hold<- matrix( NA, nrow=150, ncol=4)
  kappa<- seq( .25,1,,4)
  x2<- cbind( 0,0)
  x1<-  cbind( seq(-1,1,,150), rep( 0,150))
  for( k in 1:4){
    LKtemp<-  LKrigSetup( cbind( c( -1,1), c(-1,1)), NC=NC,
                   nlevel=nlevel,
                   a.wght= 4  + 1/(kappa[k]^2),
                   alpha=alpha)
    hold[,k]<-  LKrig.cov( x1,x2, LKinfo=LKtemp)
  }
  matplot( x1[,1], hold, type="l", lty=1, col=rainbow(5), pch=16 )
# a selection of smoothness parameters
  ktemp<- .5 # fix range
  alpha.power<- seq( 1,4,,4)
  LKtemp<- LKinfo2
  for( k in 1:4){
   LKtemp<-  LKrigSetup( cbind( c( -1,1), c(-1,1)), NC=NC,
                   nlevel=nlevel,
                   a.wght= 4  + 1/(ktemp^2),
                   alpha=alpha^alpha.power[k])
    hold[,k]<-  LKrig.cov( x1,x2, LKinfo=LKtemp)
  }
  matplot( x1[,1], hold, type="l", lty=1, col=rainbow(5) )
 set.panel()
}
 
if (FALSE) {
# generating a basis on the domain [-1,1] by [-1,1] with 1 level
# Default number of buffer points are added to each side. 
  LKinfo<- LKrigSetup(cbind( c(-1,1), c(-1,1)), NC=6,
                                 nlevel=1, a.wght=4.5,alpha=1, NC.buffer=0 )
# evaluate the basis functions on a grid to look at them
  xg<- make.surface.grid( list(x=seq(-1,1,,50), y= seq(-1,1,,50)))
  PHI<- LKrig.basis( xg,LKinfo)
  dim(PHI) # should be  2500=50^2  by  36=6^2
# plot the 9th basis function  as.surface is a handy function to
# reformat the vector as an image object
# using the grid information in an attribute of the grid points
  set.panel(1,3) 
  image.plot(as.surface(xg, PHI[,9]))
  points(  make.surface.grid( LKrigLatticeCenters(LKinfo, 1)) , col="grey", cex=.5)
  title("A radial basis function")
# compare to the tensor product basis type
  LKinfo2<- LKrigSetup(cbind( c(-1,1), c(-1,1)), NC=6,
                                 nlevel=1, a.wght=4.5,alpha=1, NC.buffer=0,
                                 BasisType="Tensor" )
  PHI2<- LKrig.basis( xg,LKinfo2)
  image.plot(as.surface(xg, PHI2[,9]))
  points(  make.surface.grid( LKrigLatticeCenters(LKinfo, 1)), col="grey", cex=.5)
  title("Tensor product basis function")
  
  image.plot(as.surface(xg, PHI[,9] - PHI2[,9]))
  points(  make.surface.grid( LKrigLatticeCenters(LKinfo, 1)), col="grey", cex=.5)
  title(" Radial - Tensor for 9th basis function")                       
set.panel()
}
#
# example of basis function indexing
#
if (FALSE) {
# generating a basis on the domain [-1,1]X[-1,1] with 3 levels
# note that there are no buffering grid points.
  set.panel(3,2)
  LKinfo<-LKrigSetup(cbind( c(-1,1), c(-1,1)), NC=6,
                    a.wght=5, alpha=c(1,.5,.25), nlevel=3,
                    NC.buffer=0)
# evaluate the basis functions on a grid to look at them
  xtemp<- seq(-1,1,,40)
  xg<- make.surface.grid( list(x=xtemp, y= xtemp) )
  PHI<- LKrig.basis( xg,LKinfo)
# coerce to dense matrix format to make plotting easier.
  PHI<- spam2full(PHI)
# first tenth, and last basis function in each resolution level
# basis functions centers are added
 set.panel(3,3)
    for(  j in 1:3){
      id1<- LKinfo$latticeInfo$offset[j]+ 1
      id2<-  LKinfo$latticeInfo$offset[j]+ 10
      idlast<- LKinfo$latticeInfo$offset[j] +
                  LKinfo$latticeInfo$mx[j,1]*LKinfo$latticeInfo$mx[j,2]
   
      centers<-  make.surface.grid(LKrigLatticeCenters(LKinfo, j) )
      image.plot( as.surface(xg, PHI[,id1]))
      points( centers, cex=.2, col="grey")
      image.plot(as.surface(xg, PHI[,id2]))
      points( centers, cex=.2, col="grey")
      image.plot( as.surface(xg, PHI[,idlast]))
      points( centers, cex=.2, col="grey")
}
  set.panel()
}
if (FALSE) {
# examining the stationarity of covariance model
  temp.fun<- 
     function( NC.buffer=0, NC=4,  a.wght=4.01){
        LKinfo<- LKrigSetup(cbind( c(-1,1), c(-1,1)),nlevel=1, alpha=1,
                                 a.wght=a.wght, NC=NC,   
                                 NC.buffer=NC.buffer,
                                  choleskyMemory=list(nnzR=2e6))
        cov1y<- cov1x<- cov0x<- cov0y<-  matrix( NA, nrow=200, ncol=20)
        cov1dx<- cov1dy<- cov0dx<- cov0dy<- matrix( NA, nrow=200, ncol=20)
        cgrid<- seq( 0,1,,20)
        for( k in 1:20){
            hold<- LKrig.cov.plot( LKinfo,
                            center=rbind( c(cgrid[k], cgrid[k])), NP=200)
            cov1x[,k] <- hold$cov[,1]
            cov1y[,k] <- hold$cov[,2]
            cov1dx[,k] <- hold$d[,1]
            cov1dy[,k] <- hold$d[,2]
            hold<- LKrig.cov.plot( LKinfo,
                             center=rbind( c(cgrid[k],0) ), NP=200)
            cov0x[,k] <- hold$cov[,1]
            cov0y[,k] <- hold$cov[,2]
            cov0dx[,k] <- hold$d[,1]
            cov0dy[,k] <- hold$d[,2]
                }
         matplot( cov1dx, cov1x, type="l", col= rainbow(20),
                         xlab="", ylab="correlation")
         mtext( side=1, line=-1, text="diagonal X")
         title( paste(  " buffer=",NC.buffer), cex=.5)
         matplot( cov1dy, cov1y, type="l", col= rainbow(20),
                        xlab="", ylab="")
         mtext( side=1, line=-1, text="diagonal Y")
         matplot(cov0dx, cov0x, type="l", col= rainbow(20),
                        xlab="",       ylab="")
         mtext( side=1, line=-1, text="middle X")
         matplot( cov0dy, cov0y, type="l", col= rainbow(20),
                         xlab="",   ylab="")
         mtext( side=1, line=-1, text="middle Y")
         title( paste( NC, a.wght), cex=.5)
}


 set.panel(3,4)
par(mar=c(3,4,1,0), oma=c(1,1,1,1))
temp.fun(  NC.buffer=5, NC=4, a.wght=4.05)
temp.fun(  NC.buffer=5, NC=16, a.wght=4.05)
temp.fun(  NC.buffer=5, NC=64, a.wght=4.05)

set.panel(4,4)
par(mar=c(3,4,1,0), oma=c(1,1,1,1))
temp.fun( NC.buffer=0, NC=8)
temp.fun( NC.buffer=2, NC=8)
temp.fun( NC.buffer=4, NC=8)
# this next one takes a while
temp.fun( NC.buffer=8,  NC=8)
# stationary == curves should all line up!

}

Run the code above in your browser using DataLab