#
# Midwest ozone data 'day 16' stripped of missings
data( ozone2)
y<- ozone2$y[16,]
good<- !is.na( y)
y<-y[good]
x<- ozone2$lon.lat[good,]
# nearly interpolate using defaults (Exponential covariance range = 2.0)
# see also mKrigMLEGrid to choose lambda by maxmimum likelihood
out<- mKrig( x,y, aRange = 2.0, lambda=.01)
out.p<- predictSurface( out)
surface( out.p)
#
# NOTE this should be identical to
# Krig( x,y, aRange=2.0, lambda=.01)
##############################################################################
# an example using a "Z" covariate and the Matern family
# again see mKrigMLEGrid to choose parameters by MLE.
data(COmonthlyMet)
yCO<- CO.tmin.MAM.climate
good<- !is.na( yCO)
yCO<-yCO[good]
xCO<- CO.loc[good,]
Z<- CO.elev[good]
out<- mKrig( xCO,yCO, Z=Z, cov.function="stationary.cov", Covariance="Matern",
aRange=4.0, smoothness=1.0, lambda=.1)
set.panel(2,1)
# quilt.plot with elevations
quilt.plot( xCO, predict(out))
# Smooth surface without elevation linear term included
surface( out)
set.panel()
#########################################################################
# here is a series of examples with bigger datasets
# using a compactly supported covariance directly
set.seed( 334)
N<- 1000
x<- matrix( 2*(runif(2*N)-.5),ncol=2)
y<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( 1000)*.1
look2<-mKrig( x,y, cov.function="wendland.cov",k=2, aRange=.2,
lambda=.1)
# take a look at fitted surface
predictSurface(look2)-> out.p
surface( out.p)
# this works because the number of nonzero elements within distance aRange
# are less than the default maximum allocated size of the
# sparse covariance matrix.
# see options() for the default values. The names follow the convention
# spam.arg where arg is the name of the spam component
# e.g. spam.nearestdistnnz
# The following will give a warning for aRange=.9 because
# allocation for the covariance matirx storage is too small.
# Here aRange controls the support of the covariance and so
# indirectly the number of nonzero elements in the sparse matrix
if (FALSE) {
look2<- mKrig( x,y, cov.function="wendland.cov",k=2, aRange=.9, lambda=.1)
}
# The warning resets the memory allocation for the covariance matrix
# according the to values options(spam.nearestdistnnz=c(416052,400))'
# this is inefficient becuase the preliminary pass failed.
# the following call completes the computation in "one pass"
# without a warning and without having to reallocate more memory.
options( spam.nearestdistnnz=c(416052,400))
look2<- mKrig( x,y, cov.function="wendland.cov",k=2,
aRange=.9, lambda=1e-2)
# as a check notice that
# print( look2)
# reports the number of nonzero elements consistent with the specifc allocation
# increase in spam.options
# new data set of 1500 locations
set.seed( 234)
N<- 1500
x<- matrix( 2*(runif(2*N)-.5),ncol=2)
y<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( N)*.01
if (FALSE) {
# the following is an example of where the allocation (for nnzR)
# for the cholesky factor is too small. A warning is issued and
# the allocation is increased by 25
#
look2<- mKrig( x,y,
cov.function="wendland.cov",k=2, aRange=.1, lambda=1e2 )
}
# to avoid the warning
look2<-mKrig( x,y,
cov.function="wendland.cov", k=2, aRange=.1,
lambda=1e2, chol.args=list(pivot=TRUE, memory=list(nnzR= 450000)))
###############################################################################
# fiting multiple data sets
#
#\dontrun{
y1<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( N)*.01
y2<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( N)*.01
Y<- cbind(y1,y2)
look3<- mKrig( x,Y,cov.function="wendland.cov",k=2, aRange=.1,
lambda=1e2 )
# note slight difference in summary because two data sets have been fit.
print( look3)
#}
if (FALSE) {
##################################################################
# finding a good choice for aRange as a taper
# Suppose the target is a spatial prediction using roughly 50 nearest neighbors
# (tapering covariances is effective for roughly 20 or more in the situation of
# interpolation) see Furrer, Genton and Nychka (2006).
# take a look at a random set of 100 points to get idea of scale
# and saving computation time by not looking at the complete set
# of points
# NOTE: This could also be done directly using the FNN package for finding nearest
# neighbors
set.seed(223)
ind<- sample( 1:N,100)
hold<- rdist( x[ind,], x)
dd<- apply( hold, 1, quantile, p= 50/N )
dguess<- max(dd)
# dguess is now a reasonable guess at finding cutoff distance for
# 50 or so neighbors
# full distance matrix excluding distances greater than dguess
hold2<- nearest.dist( x, x, delta= dguess )
# here is trick to find the number of nonsero rows for a matrix in spam format.
hold3<- diff( hold2@rowpointers)
# min( hold3) = 43 which we declare close enough. This also counts the diagonal
# So there are a minimum of 42 nearest neighbors ( median is 136)
# see table( hold3) for the distribution
# now the following will use no less than 43 - 1 nearest neighbors
# due to the tapering.
mKrig( x,y, cov.function="wendland.cov",k=2, aRange=dguess,
lambda=1e2) -> look2
}
###############################################################################
# use precomputed distance matrix
#
if (FALSE) {
y1<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( N)*.01
y2<- sin( 1.8*pi*x[,1])*sin( 2.5*pi*x[,2]) + rnorm( N)*.01
Y<- cbind(y1,y2)
#precompute distance matrix in compact form
distMat = rdist(x, compact=TRUE)
look3<- mKrig( x,Y,cov.function="stationary.cov", aRange=.1,
lambda=1e2, distMat=distMat )
#precompute distance matrix in standard form
distMat = rdist(x)
look3<- mKrig( x,Y,cov.function="stationary.cov", aRange=.1,
lambda=1e2, distMat=distMat )
}
Run the code above in your browser using DataLab