# NOTE: CRAN requires that the "run" examples execute within 5 seconds.
# so to meet this constraint many of these have been
# severely limited to be quick, typically making NC and nlevel
# very small.
# 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]
# fairly arbitrary choices for covariance parameters and lambda
# just to show a basic level call
obj1<- LKrig( x,y, a.wght=5, nlevel=3, nu=1.0, NC=10, lambda=.1)
print(obj1)
# estimates of fixed model parameters with their SEs
summary( obj1)$coefficients
#
# this is the same as:
# LKinfoEX<- LKrigSetup( x, a.wght=5, nlevel=3, nu=1.0, NC=4)
# obj1<- LKrig( x,y, LKinfo= LKinfoEX, lambda=.1)
# thin plate spline-like model with the lambda parameter estimated by
# maximum likelihood. Default choices are made for a.wght, nlevel, NC
# and alpha.
if (FALSE) {
obj<- LatticeKrig( x, y)
# summary of fit and a plot of fitted surface
print( obj)
surface( obj )
US(add=TRUE)
points(x)
# prediction standard errors
out.se<- predictSE( obj, xnew= x)
}
# breaking down the LatticeKrig function into several steps.
# also use a different covariance model that has fewer basis functions
# (to make the example run more quickly)
if (FALSE) {
LKinfo<- LKrigSetup( x, nlevel=3, nu=1, NC=5, a.wght=5,
lambda=.01)
# maximize likelihood over lambda see help( LKrig.MLE) for details
# this assumes the value of 5 for a.wght. In many cases the fit is not
# very sensitive to the range parameter such as a.wght in this case --
# but very sensitive to lambda when varied on a log scale.
MLE.fit<- LKrig.MLE(x,y, LKinfo=LKinfo)
MLE.fit$summary # summary of optimization over lambda.
# fit using MLE for lambda MLE function has added MLE value of lambda to
# the LKinfo object.
obj<- LKrig( x,y, LKinfo=MLE.fit$LKinfo.MLE)
print( obj)
# find prediction standard errors at locations based on fixing covariance
# at MLE's
out.se<- predictSE( obj, xnew= x)
# one could evaluate the SE on a grid to get the surface of predicted SE's
# for large grids it is better to use LKrig.sim.conditional to estimate
# the variances by Monte Carlo
}
##########################################################################
# Use multi-resolution model that approximates an exponential covariance
# Note that a.wght, related to a range/scale parameter
# is specified at a (seemingly) arbitrary value.
##########################################################################
LKinfo<- LKrigSetup( x, NC=4, nu=1, nlevel=2, a.wght= 5)
# take a look at the implied covariance function solid=along x
# and dashed=along y
check.cov<- LKrig.cov.plot( LKinfo)
matplot( check.cov$d, check.cov$cov, type="l", lty=c(1,2))
############################################################################
# Search over lambda to find MLE for ozone data with approximate exponential
# covariance
###########################################################################
if (FALSE) {
LKinfo.temp<- LKrigSetup( x, NC=6, nu=1, nlevel=3, a.wght= 5, lambda=.5)
# starting value for lambda optimization
MLE.search<- LKrig.MLE(x,y, LKinfo=LKinfo.temp)
# this function returns an LKinfo object with the MLE for lambda included.
MLE.ozone.fit<- LKrig( x,y, LKinfo= MLE.search$LKinfo.MLE)
}
###########################################################################
# Including a covariate (linear fixed part in spatial model)
##########################################################################
data(COmonthlyMet)
y.CO<- CO.tmin.MAM.climate
good<- !is.na( y.CO)
y.CO<-y.CO[good]
x.CO<- as.matrix(CO.loc[good,])
Z.CO<- CO.elev[good]
# single level with large range parameter -- similar to a thin plate spline
# lambda specified
# fit with elevation
# note how for convenience the LKrig parameters are
# just included here and not passed as a separate LKinfo object.
obj.CO.elev<- LKrig( x.CO,y.CO,Z=Z.CO, nlevel=1, NC=50, alpha=1, lambda=.005,
a.wght=4.1)
# BTW the coefficient for the linear term for elevation is obj.CO$d.coef[4]
# fitted surface without the elevation term
if (FALSE) {
LKinfo<- LKrigSetup( x.CO, nlevel=1, NC=20,alpha=1, a.wght=4.1, lambda=1.0)
# lambda given here is just the starting value for MLE optimization
CO.MLE<- LKrig.MLE( x.CO,y.CO,Z=Z.CO, LKinfo=LKinfo)
obj.CO.elev<- LKrig( x.CO,y.CO,Z=Z.CO, LKinfo= CO.MLE$LKinfo.MLE)
CO.surface2<- predictSurface( obj.CO.elev, drop.Z=TRUE, nx=50, ny=50)
# pull off CO elevations at same locations on grid as the surface
data( RMelevation)
# a superset of elevations at 4km resolution
elev.surface<- interp.surface.grid( RMelevation, CO.surface2)
CO.full<- predictSurface( obj.CO.elev, ZGrid= elev.surface, nx=50, ny=50)
# for comparison fit without elevation as a linear covariate:
CO.MLE2<- LKrig.MLE( x.CO,y.CO, LKinfo=LKinfo)
obj.CO<- LKrig( x.CO,y.CO, LKinfo= CO.MLE2$LKinfo.MLE)
# surface estimate
CO.surface<- predictSurface( obj.CO, nx=50, ny=50)
set.panel( 2,1)
coltab<- topo.colors(256)
zr<- range( c( CO.full$z), na.rm=TRUE)
image.plot( CO.surface, col=coltab, zlim =zr)
US( add=TRUE,lwd=2)
title( "MAM min temperatures without elevation")
image.plot( CO.full, col=coltab, zlim=zr)
title( "Including elevation")
US( add=TRUE,lwd=2)
}
if (FALSE) {
#### a 3-d example using alternative geometry
set.seed( 123)
N<- 500
x<- matrix( runif(3* N,-1,1), ncol=3, nrow=N)
y<- 10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)
LKinfo<- LKrigSetup( x,
LKGeometry = "LKBox",
nlevel = 1,
a.wght = 6.01,
NC = 5,
NC.buffer = 2,
normalize = TRUE,
choleskyMemory = list(nnzR= 2e6),
lambda = .1,
mean.neighbor=200
)
# NOTE: memory for the spam routines has been increased to
# avoid warnings
out1<- LKrig( x,y, LKinfo=LKinfo)
}
if (FALSE) {
# MLE search over lambda and a bigger problem
# but no normalization
N<- 5e4
x<- matrix( runif(3* N,-1,1), ncol=3, nrow=N)
y<- 10*exp( -rdist( x, rbind( c(.5,.5,.6) ) )/.5)
LKinfo<- LKrigSetup( x,
LKGeometry = "LKBox",
nlevel = 1,
a.wght = 6.01,
NC = 20,
NC.buffer = 2,
normalize = FALSE,
choleskyMemory = list(nnzR= 25e6),
lambda = .1,
mean.neighbor=200
)
# add in timing
system.time( out1<- LKrig( x,y, LKinfo=LKinfo) ) # about 25 seconds
# LKinfo$normalize<- TRUE
# system.time( out2<- LatticeKrig( x,y, LKinfo=LKinfo) )# ~250 seconds
}
########################################################################
# for a more elaborate search over a.wght, alpha and lambda to find
# joint MLEs see help(LKrig.MLE)
########################################################################
########################################################################
# A bigger problem: 26K observations and 4.6K basis functions
# fitting takes about 15 seconds on a laptop for a fixed covariance
# LKrig.MLE to find the MLE (not included) for lambda takes about
# 8 minutes
#######################################################################
if (FALSE) {
data(CO2)
# the Ring geometry will be periodic in the first dimension and rectangular on
# second. A useful approximation for spherical data omitting the polar caps.
LKinfo.CO2<- LKrigSetup(CO2$lon.lat, NC=100,nlevel=1, lambda=.2,
a.wght=5, alpha=1,
LKGeometry="LKRing", choleskyMemory = list(nnzR=2e6) )
print(LKinfo.CO2)
obj1<- LKrig( CO2$lon.lat,CO2$y,LKinfo=LKinfo.CO2)
# 5700+ basis functions 101X57 lattice (number of basis functions
# reduced in y direction because of a rectangular domain
obj1$trA.est # about 2900+ effective degrees of freedom
#
glist<- list( x= seq( -180,180,,200),y=seq( -80,80,,100) )
fhat<- predictSurface( obj1,grid.list=glist)
#Plot data and gap-filled estimate
set.panel(2,1)
quilt.plot(CO2$lon.lat,CO2$y,zlim=c(373,381))
title("Simulated CO2 satellite observations")
world(add=TRUE,col="magenta")
image.plot( fhat,zlim=c(373,381))
world( add=TRUE, col="magenta")
title("Gap-filled global predictions")
}
set.panel()
#########################################################################
# Comparing LKrig to ordinary Kriging
########################################################################
# Here is an illustration of using the fields function mKrig with the
# LKrig covariance to reproduce the computations of LKrig. The
# difference is that mKrig can not take advantage of any sparsity in
# the precision matrix because its inverse, the covariance matrix, is
# not sparse. This example reinforces the concept that LKrig finds the
# the standard geostatistical estimate but just uses a particular
# covariance function defined via basis functions and the precision
# matrix.
# Load ozone data set (AGAIN)
if (FALSE) {
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]
a.wght<- 5
lambda <- 1.5
obj1<- LKrig( x,y,NC=16,nlevel=1, alpha=1, lambda=lambda, a.wght=5,
NtrA=20,iseed=122)
# in both calls iseed is set the same so we can compare
# Monte Carlo estimates of effective degrees of freedom
obj1$trA.est
# The covariance "parameters" are all in the list LKinfo
# to create this special list outside of a call to LKrig use
LKinfo.test <- LKrigSetup( x, NC=16, nlevel=1, alpha=1.0, a.wght=5)
# this call to mKrig should be identical to the LKrig results
# because it uses the LKrig.cov covariance with all the right parameters.
obj2<- mKrig( x,y, lambda=lambda, m=2, cov.function="LKrig.cov",
cov.args=list( LKinfo=LKinfo.test), NtrA=20, iseed=122)
# compare the two results this is also an
# example of how tests are automated in fields
# set flag to have tests print results
test.for.zero.flag<- TRUE
test.for.zero( obj1$fitted.values, obj2$fitted.values,
tag="comparing predicted values LKrig and mKrig")
# compare standard errors.
se1<- predictSE.LKrig( obj1)
se2<- predictSE.mKrig(obj2)
test.for.zero( se1, se2,
tag="comparing standard errors for LKrig and mKrig")
}
if (FALSE) {
##################################################################
# a linear inverse problem
##################################################################
# NOTE the settings in the model are small to make for a fast example
#
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
good<- !is.na(y)
x<- x[good,]
y<- y[good]
LKinfo<- LKrigSetup(x, a.wght=4.5, nlevel=3, nu=1, NC=4, lambda=.01)
# now observe a linear combination
NNDist<- LKDist(x,x, delta=1.5)
A<- NNDist
A$ra<- exp(-NNDist$ra)
# A is a weight matrix based on neighbors close by and
# in spind sparse matrix format
# now convert to spam format
A<- spind2spam(A)
TMatrix <- get(LKinfo$fixedFunction)(x = x)
# Tmatrix is a 3 column matrix of constant and the two spatial coordinates
# i.e. a linear function of the spatial variables
PHI<- LKrig.basis(x, LKinfo)
X<- A%*% PHI
U<- A%*%TMatrix
yIndirect<- A%*%y
#
# X<- A
out0<- LatticeKrig(x,y, LKinfo=LKinfo)
out1<- LatticeKrig(x,yIndirect, U=U, X=X, LKinfo=LKinfo)
# the predict function evaluates f in this case -- not the fitted values that
# are related to the
# observations
# partial agreement but not exact -- this is because the observational models
# assume different errors
#
plot( predict(out0,x), predict( out1,x))
out<- LKrig.sim.conditional( out1,M=5, nx=10, ny=10 )
}
Run the code above in your browser using DataLab