-
Notifications
You must be signed in to change notification settings - Fork 8
/
exp.cov.html
451 lines (385 loc) · 16.6 KB
/
exp.cov.html
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
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns="http://www.w3.org/1999/xhtml"><head><title>R: Exponential family, radial basis functions,cubic spline,...</title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<link rel="stylesheet" type="text/css" href="R.css" />
</head><body>
<table width="100%" summary="page for Covariance functions {fields}"><tr><td>Covariance functions {fields}</td><td style="text-align: right;">R Documentation</td></tr></table>
<h2>
Exponential family, radial basis
functions,cubic spline, compactly supported Wendland family and
stationary covariances. </h2>
<h3>Description</h3>
<p>Given two sets of locations these functions compute the cross covariance matrix for
some covariance families. In addition these functions can take advantage
of spareness, implement more efficient multiplcation of the
cross covariance by a vector or matrix and also return a marginal
variance to be consistent with calls by the Krig function.
</p>
<p><code>stationary.cov</code> and <code>Exp.cov</code> have additional arguments for
precomputed distance matrices and for calculating only the upper triangle
and diagonal of the output covariance matrix to save time. Also, they
support using the <code>rdist</code> function with <code>compact=TRUE</code> or input
distance matrices in compact form, where only the upper triangle of the
distance matrix is used to save time.
</p>
<p>Note: These functions have been been renamed from the previous fields functions
using 'Exp' in place of 'exp' to avoid conflict with the generic exponential
function (<code>exp(...)</code>)in R.
</p>
<h3>Usage</h3>
<pre>
Exp.cov(x1, x2=NULL, theta = 1, p=1, distMat = NA,
C = NA, marginal = FALSE, onlyUpper=FALSE)
Exp.simple.cov(x1, x2, theta =1, C=NA,marginal=FALSE)
Rad.cov(x1, x2, p = 1, m=NA, with.log = TRUE, with.constant = TRUE,
C=NA,marginal=FALSE, derivative=0)
cubic.cov(x1, x2, theta = 1, C=NA, marginal=FALSE)
Rad.simple.cov(x1, x2, p=1, with.log = TRUE, with.constant = TRUE,
C = NA, marginal=FALSE)
stationary.cov(x1, x2=NULL, Covariance = "Exponential", Distance = "rdist",
Dist.args = NULL, theta = 1, V = NULL, C = NA, marginal = FALSE,
derivative = 0, distMat = NA, onlyUpper = FALSE, ...)
stationary.taper.cov(x1, x2, Covariance="Exponential",
Taper="Wendland",
Dist.args=NULL, Taper.args=NULL,
theta=1.0,V=NULL, C=NA, marginal=FALSE,
spam.format=TRUE,verbose=FALSE,...)
wendland.cov(x1, x2, theta = 1, V=NULL, k = 2, C = NA,
marginal =FALSE,Dist.args = list(method = "euclidean"),
spam.format = TRUE, derivative = 0, verbose=FALSE)
</pre>
<h3>Arguments</h3>
<table summary="R argblock">
<tr valign="top"><td><code>x1</code></td>
<td>
<p> Matrix of first set of locations where each row gives the
coordinates of a particular point.</p>
</td></tr>
<tr valign="top"><td><code>x2</code></td>
<td>
<p> Matrix of second set of locations where each row gives the
coordinatesof a particular point. If this is missing x1 is used. </p>
</td></tr>
<tr valign="top"><td><code>theta</code></td>
<td>
<p> Range (or scale) parameter. This should be a scalar (use
the V argument for other scaling options). Any distance calculated for
a covariance function is divided by theta before the covariance function
is evaluated.</p>
</td></tr>
<tr valign="top"><td><code>V</code></td>
<td>
<p> A matrix that describes the inverse linear transformation of
the coordinates before distances are found. In R code this
transformation is: <code>x1 %*% t(solve(V))</code> Default is NULL, that
is the transformation is just dividing distance by the scalar value
<code>theta</code>. See Details below. If one has a vector of "theta's"
that are the scaling for each coordinate then just express this as
<code>V = diag(theta)</code> in the call to this function.</p>
</td></tr>
<tr valign="top"><td><code>C</code></td>
<td>
<p> A vector with the same length as the number of rows of x2.
If specified the covariance matrix will be multiplied by this vector.</p>
</td></tr>
<tr valign="top"><td><code>marginal</code></td>
<td>
<p>If TRUE returns just the diagonal elements of the
covariance matrix using the <code>x1</code> locations. In this case this is
just 1.0. The marginal argument will trivial for this function is a
required argument and capability for all covariance functions used
with Krig.</p>
</td></tr>
<tr valign="top"><td><code>p</code></td>
<td>
<p> Exponent in the exponential covariance family. p=1 gives an
exponential and p=2 gives a Gaussian. Default is the exponential
form. For the radial basis function this is the exponent applied to
the distance between locations.</p>
</td></tr>
<tr valign="top"><td><code>m</code></td>
<td>
<p>For the radial basis function p = 2m-d, with d being the dimension of the
locations, is the exponent applied to
the distance between locations. (m is a common way of parametrizing this exponent.)</p>
</td></tr>
<tr valign="top"><td><code>with.constant</code></td>
<td>
<p> If TRUE includes complicated constant for radial
basis functions. See the function <code>radbad.constant</code> for more
details. The default is TRUE, include the constant. Without the usual
constant the lambda used here will differ by a constant from spline
estimators ( e.g. cubic smoothing splines) that use the
constant. Also a negative value for the constant may be necessary to
make the radial basis positive definite as opposed to negative
definite. </p>
</td></tr>
<tr valign="top"><td><code>with.log</code></td>
<td>
<p>If TRUE include a log term for even dimensions. This
is needed to be a thin plate spline of integer order. </p>
</td></tr>
<tr valign="top"><td><code>Covariance</code></td>
<td>
<p>Character string that is the name of the covariance
shape function for the distance between locations. Choices in fields
are <code>Exponential</code>, <code>Matern</code></p>
</td></tr>
<tr valign="top"><td><code>Distance</code></td>
<td>
<p>Character string that is the name of the distance
function to use. Choices in fields are <code>rdist</code>,
<code>rdist.earth</code></p>
</td></tr>
<tr valign="top"><td><code>Taper</code></td>
<td>
<p>Character string that is the name of the taper function
to use. Choices in fields are listed in help(taper).</p>
</td></tr>
<tr valign="top"><td><code>Dist.args</code></td>
<td>
<p> A list of optional arguments to pass to the Distance
function.</p>
</td></tr>
<tr valign="top"><td><code>Taper.args</code></td>
<td>
<p> A list of optional arguments to pass to the Taper
function. <code>theta</code> should always be the name for the range (or
scale) paremeter.</p>
</td></tr>
<tr valign="top"><td><code>spam.format</code></td>
<td>
<p>If TRUE returns matrix in sparse matrix format
implemented in the spam package. If FALSE just returns a full
matrix. </p>
</td></tr>
<tr valign="top"><td><code>k</code></td>
<td>
<p>The order of the Wendland covariance function. See help on
Wendland.</p>
</td></tr>
<tr valign="top"><td><code>derivative</code></td>
<td>
<p> If nonzero evaluates the partials of the
covariance function at locations x1. This must be used with the "C" option
and is mainly called from within a predict function. The partial
derivative is taken with respect to <code>x1</code>. </p>
</td></tr>
<tr valign="top"><td><code>verbose</code></td>
<td>
<p>If TRUE prints out some useful information for
debugging.</p>
</td></tr>
<tr valign="top"><td><code>distMat</code></td>
<td>
<p>If the distance matrix between <code>x1</code> and <code>x2</code> has already been
computed, it can be passed via this argument so it won't need to be
recomputed.
</p>
</td></tr>
<tr valign="top"><td><code>onlyUpper</code></td>
<td>
<p>For internal use only, not meant to be set by the user. Automatically
set to <code>TRUE</code> by <code>mKrigMLEJoint</code> or <code>mKrigMLEGrid</code> if
<code>lambda.profile</code> is set to <code>TRUE</code>, but set to <code>FALSE</code>
for the final parameter fit so output is compatible with rest of
<code>fields</code>.
</p>
<p>If <code>TRUE</code>, only the upper triangle and diagonal of the covariance
matrix is computed to save time (although if a non-compact distance
matrix is used, the onlyUpper argument is set to FALSE). If <code>FALSE</code>,
the entire covariance matrix is computed. In general, it should
only be set to <code>TRUE</code> for <code>mKrigMLEJoint</code> and <code>mKrigMLEGrid</code>,
and the default is set to <code>FALSE</code> so it is compatible with all of
<code>fields</code>.
</p>
</td></tr>
<tr valign="top"><td><code>...</code></td>
<td>
<p> Any other arguments that will be passed to the
covariance function. e.g. <code>smoothness</code> for the Matern.</p>
</td></tr> </table>
<h3>Details</h3>
<p> For purposes of illustration, the function
<code>Exp.cov.simple</code> is provided in fields as a simple example and
implements the R code discussed below. List this function out as a
way to see the standard set of arguments that fields uses to define a
covariance function. It can also serve as a template for creating new
covariance functions for the <code>Krig</code> and <code>mKrig</code>
functions. Also see the higher level function <code>stationary.cov</code> to
mix and match different covariance shapes and distance functions.
</p>
<p>A common scaling for stationary covariances: If <code>x1</code> and
<code>x2</code> are matrices where <code>nrow(x1)=m</code> and <code>nrow(x2)=n</code>
then this function will return a mXn matrix where the (i,j) element
is the covariance between the locations <code>x1[i,]</code> and
<code>x2[j,]</code>. The exponential covariance function is computed as
exp( -(D.ij)) where D.ij is a distance between <code>x1[i,]</code> and
<code>x2[j,]</code> but having first been scaled by theta. Specifically if
<code>theta</code> is a matrix to represent a linear transformation of the
coordinates, then let <code>u= x1%*% t(solve( theta))</code> and <code>v=
x2%*% t(solve(theta))</code>. Form the mXn distance matrix with
elements:
</p>
<p><code>D[i,j] = sqrt( sum( ( u[i,] - v[j,])**2 ) )</code>.
</p>
<p>and the cross covariance matrix is found by <code>exp(-D)</code>. The
tapered form (ignoring scaling parameters) is a matrix with i,j entry
<code>exp(-D[i,j])*T(D[i,j])</code>. With T being a positive definite
tapering function that is also assumed to be zero beyond 1.
</p>
<p>Note that if theta is a scalar then this defines an isotropic
covariance function and the functional form is essentially
<code>exp(-D/theta)</code>.
</p>
<p>Implementation: The function <code>r.dist</code> is a useful FIELDS function
that finds the cross Euclidean distance matrix (D defined above) for
two sets of locations. Thus in compact R code we have
</p>
<p>exp(-rdist(u, v))
</p>
<p>Note that this function must also support two other kinds of calls:
</p>
<p>If marginal is TRUE then just the diagonal elements are returned (in R
code <code>diag( exp(-rdist(u,u)) )</code>).
</p>
<p>If C is passed then the returned value is <code> exp(-rdist(u, v))
%*% C</code>.
</p>
<p>Some details on particular covariance functions:
</p>
<dl>
<dt>Radial basis functions (<code>Rad.cov</code>:</dt><dd><p>The
functional form is Constant* rdist(u, v)**p for odd dimensions and
Constant* rdist(u,v)**p * log( rdist(u,v) ) For an m th order thin plate
spline in d dimensions p= 2*m-d and must be positive. The constant,
depending on m and d, is coded in the fields function
<code>radbas.constant</code>. This form is only a generalized covariance
function – it is only positive definite when restricted to linear
subspace. See <code>Rad.simple.cov</code> for a coding of the radial basis
functions in R code.</p>
</dd>
<dt>Stationary covariance <code>stationary.cov</code>:</dt><dd><p>Here the
computation is to apply the function Covariance to the distances found
by the Distance function. For example
</p>
<p><code>Exp.cov(x1,x2, theta=MyTheta)</code>
</p>
<p>and
</p>
<p><code>stationary.cov( x1,x2, theta=MyTheta, Distance= "rdist",
Covariance="Exponential")</code>
</p>
<p>are the same. This also the same as
</p>
<p><code>stationary.cov( x1,x2, theta=MyTheta, Distance= "rdist",
Covariance="Matern",smoothness=.5)</code>. </p>
</dd>
<dt>Stationary tapered covariance <code>stationary.taper.cov</code>:</dt><dd><p>The
resulting cross covariance is the direct or Shure product of the
tapering function and the covariance. In R code given location
matrices, <code>x1</code> and <code>x2</code> and using Euclidean distance.
</p>
<p><code>Covariance(rdist( x1, x2)/theta)*Taper( rdist( x1,
x2)/Taper.args$theta)</code>
</p>
<p>By convention, the <code>Taper</code> function is assumed to be identically
zero outside the interval [0,1]. Some efficiency is introduced within
the function to search for pairs of locations that are nonzero with
respect to the Taper. This is done by the SPAM function
<code>nearest.dist</code>. This search may find more nonzero pairs than
dimensioned internally and SPAM will try to increase the space. One
can also reset the SPAM options to avoid these warnings. For
spam.format TRUE the multiplication with the <code>C</code> argument is done
with the spam sparse multiplication routines through the "overloading"
of the <code>%*%</code> operator. </p>
</dd>
</dl>
<p>About the FORTRAN: The actual function <code>Exp.cov</code> and
<code>Rad.cov</code> call FORTRAN to
make the evaluation more efficient this is especially important when the
C argument is supplied. So unfortunately the actual production code in
Exp.cov is not as crisp as the R code sketched above. See
<code>Rad.simple.cov</code> for a R coding of the radial basis functions.
</p>
<h3>Value</h3>
<p> If the argument C is NULL the cross covariance matrix is
returned. In general if nrow(x1)=m and nrow(x2)=n then the returned
matrix will be mXn. Moreover, if x1 is equal to x2 then this is the
covariance matrix for this set of locations.
</p>
<p>If C is a vector of length n, then returned value is the
multiplication of the cross covariance matrix with this vector.
</p>
<h3>See Also</h3>
<p>Krig, rdist, rdist.earth, gauss.cov, Exp.image.cov, Exponential, Matern,
Wendland.cov, mKrig</p>
<h3>Examples</h3>
<pre>
# exponential covariance matrix ( marginal variance =1) for the ozone
#locations
out<- Exp.cov( ChicagoO3$x, theta=100)
# out is a 20X20 matrix
out2<- Exp.cov( ChicagoO3$x[6:20,],ChicagoO3$x[1:2,], theta=100)
# out2 is 15X2 matrix
# Kriging fit where the nugget variance is found by GCV
# Matern covariance shape with range of 100.
#
fit<- Krig( ChicagoO3$x, ChicagoO3$y, Covariance="Matern", theta=100,smoothness=2)
data( ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
# Omit the NAs
good<- !is.na( y)
x<- x[good,]
y<- y[good]
# example of calling the taper version directly
# Note that default covariance is exponential and default taper is
# Wendland (k=2).
stationary.taper.cov( x[1:3,],x[1:10,] , theta=1.5, Taper.args= list(k=2,theta=2.0,
dimension=2) )-> temp
# temp is now a tapered 3X10 cross covariance matrix in sparse format.
is.spam( temp) # evaluates to TRUE
# should be identical to
# the direct matrix product
temp2<- Exp.cov( x[1:3,],x[1:10,], theta=1.5) * Wendland(rdist(x[1:3,],x[1:10,]),
theta= 2.0, k=2, dimension=2)
test.for.zero( as.matrix(temp), temp2)
# Testing that the "V" option works as advertized ...
x1<- x[1:20,]
x2<- x[1:10,]
V<- matrix( c(2,1,0,4), 2,2)
Vi<- solve( V)
u1<- t(Vi%*% t(x1))
u2<- t(Vi%*% t(x2))
look<- exp(-1*rdist(u1,u2))
look2<- stationary.cov( x1,x2, V= V)
test.for.zero( look, look2)
# Here is an example of how the cross covariance multiply works
# and lots of options on the arguments
Ctest<- rnorm(10)
temp<- stationary.cov( x,x[1:10,], C= Ctest,
Covariance= "Wendland",
k=2, dimension=2, theta=1.5 )
# do multiply explicitly
temp2<- stationary.cov( x,x[1:10,],
Covariance= "Wendland",
k=2, dimension=2, theta=1.5 )%*% Ctest
test.for.zero( temp, temp2)
# use the tapered stationary version
# cov.args is part of the argument list passed to stationary.taper.cov
# within Krig.
# This example needs the spam package.
#
## Not run:
Krig(x,y, cov.function = "stationary.taper.cov", theta=1.5,
cov.args= list(Taper.args= list(k=2, dimension=2,theta=2.0) )
) -> out2
# NOTE: Wendland is the default taper here.
## End(Not run)
# BTW this is very similar to
## Not run:
Krig(x,y, theta= 1.5)-> out
## End(Not run)
</pre>
<hr /><div style="text-align: center;">[Package <em>fields</em> version 9.9 <a href="00Index.html">Index</a>]</div>
</body></html>