LKrig.basis: Functions for generating a multi-resolution, compactly...
In LatticeKrig: Multi-Resolution Kriging Based on Markov Random Fields
Description
Usage
Arguments
Details
Value
Author(s)
See Also
Examples
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
1
2
3
4
5
6
7
8
9
10
#
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 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.)
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).
return.B
If TRUE B is returned instead of the precision matrix
t(B)%*%B.
Q
A sparse (spam format) precision matrix.
theta
Currently theta is a place holder argument
for future development.
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.
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
1
2
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:
1
2
3
4
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:
1
2
3
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
2
3
4
5
6
7
8
__________
. -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
1
2
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
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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
# 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)
LatticeKrig documentation built on Nov. 9, 2019, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Details Value Author(s) See Also Examples
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
1 2 3 4 5 6 7 8 9 10 | #
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
1 2 |
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:
1 2 3 4 | 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:
1 2 3 |
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 2 3 4 5 6 7 8 | __________
. -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
1 2 | 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
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 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 | # 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)
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.