| LKrig.basis {LatticeKrig} | R Documentation |
Functions for generating a multi-resolution, compactly supported basis, multi-resolution covariance functions and simulating from these processes.
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,
theta = NULL)
LKrig.cov.plot( LKinfo, NP=200, center = NULL, xlim = NULL, ylim = NULL)
LKrigCovWeightedObs(x1, wX, LKinfo)
LKrig.basis(x1, LKinfo, 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)
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 |
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 |
choleskyMemory |
A list giving the memory requirements for
a sparse cholesky decomposition. See |
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
|
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 |
PHI |
A sparse matrix of basis functions (rows index points for evaluation, columns index basis functions). |
return.B |
If TRUE B is returned instead of the precision matrix t(B)%*%B. |
Q |
A sparse (spam format) precision matrix. |
theta |
Currently |
wX |
The wX matrix constructed in the |
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 |
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. |
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
LKrig.normalize.basis.slow takes advantage of the sparsity of the
precision matrix for the coefficients.
LKrig.normalize.basis.fast is even 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.
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.
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$rho.object
is NULL, which it usually is.) .
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.
LKrig.normalize.basis, LKrig.normalize.basis.fast, LKrig.normalize.basis.slow: 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.
Author(s)
Doug Nychka
See Also
LKrig, mKrig, Krig, fastTps, Wendland, LKrigSAR
Examples
# 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)
#
#
## Not run:
# 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.
#
## End(Not run)
## Not run:
#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()
## End(Not run)
## Not run:
# 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()
## End(Not run)
#
# example of basis function indexing
#
## Not run:
# 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()
## End(Not run)
## Not run:
# 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!
## End(Not run)