-
Notifications
You must be signed in to change notification settings - Fork 8
/
Krig.html
848 lines (719 loc) · 26.4 KB
/
Krig.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
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
<!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: Kriging surface estimate</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 Krig {fields}"><tr><td>Krig {fields}</td><td style="text-align: right;">R Documentation</td></tr></table>
<h2>
Kriging surface estimate
</h2>
<h3>Description</h3>
<p>Fits a surface to irregularly spaced data. The Kriging model assumes
that the unknown function is a realization of a Gaussian
random spatial processes. The assumed model is additive
Y = P(x) + Z(X) + e, where P is a low order polynomial and Z is a
mean zero,
Gaussian stochastic process with a
covariance that is unknown up to a scale constant. The main advantages of
this function are the flexibility in specifying the covariance as an
R language function and also the supporting functions plot, predict,
predictSE, surface for
subsequent analysis. Krig also supports a correlation model where the mean
and marginal variances are supplied.
</p>
<h3>Usage</h3>
<pre>
Krig(x, Y, cov.function = "stationary.cov", lambda = NA, df
= NA, GCV = FALSE, Z = NULL, cost = 1, knots = NA,
weights = NULL, m = 2, nstep.cv = 200, scale.type =
"user", x.center = rep(0, ncol(x)), x.scale = rep(1,
ncol(x)), rho = NA, sigma2 = NA, method = "REML",
verbose = FALSE, mean.obj = NA, sd.obj = NA,
null.function = "Krig.null.function", wght.function =
NULL, offset = 0, na.rm = TRUE, cov.args = NULL,
chol.args = NULL, null.args = NULL, wght.args = NULL,
W = NULL, give.warnings = TRUE, ...)
## S3 method for class 'Krig'
fitted(object,...)
## S3 method for class 'Krig'
coef(object,...)
resid.Krig(object,...)
</pre>
<h3>Arguments</h3>
<table summary="R argblock">
<tr valign="top"><td><code>chol.args</code></td>
<td>
<p>Arguments to be passed to the cholesky decomposition in Krig.engine.fixed.
The default if NULL, assigned at the top level of this function, is
list( pivot=FALSE). This argument is useful when working with
the sparse matrix package. </p>
</td></tr>
<tr valign="top"><td><code>cov.args</code></td>
<td>
<p>A list with the arguments to call the covariance function. (in addition to the locations)
</p>
</td></tr>
<tr valign="top"><td><code>cov.function</code></td>
<td>
<p>Covariance function for data in the form of an R function (see
Exp.simple.cov as an example).
Default assumes that correlation is an exponential function of distance.
See also <code>stationary.cov</code> for more general choice of covariance
shapes. <code>exponential.cov</code> will be faster if only the exponential
covariance form is needed.
</p>
</td></tr>
<tr valign="top"><td><code>cost</code></td>
<td>
<p>Cost value used in GCV criterion. Corresponds to a penalty for
increased number of parameters. The default is 1.0 and corresponds to the
usual GCV function.
</p>
</td></tr>
<tr valign="top"><td><code>df</code></td>
<td>
<p>The effective number of parameters for the fitted surface. Conversely,
N- df, where N is the total number of observations is the degrees of
freedom associated with the residuals.
This is an alternative to specifying lambda and much more interpretable.
NOTE: GCV argument defaults to TRUE if this argument is used.
</p>
</td></tr>
<tr valign="top"><td><code>GCV</code></td>
<td>
<p> If TRUE matrix decompositions are done to allow estimating
lambda by GCV or REML and specifying smoothness by the effective degrees of
freedom. So the GCV switch does more than just supply a GCV estimate. Also if
lambda or df are passed the estimate will be evaluated at those values,
not at the GCV/REML estimates of lambda.
If FALSE Kriging estimate is found under a fixed lambda model. </p>
</td></tr>
<tr valign="top"><td><code>give.warnings</code></td>
<td>
<p> If TRUE warnings are given in gcv grid search limits.
If FALSE warnings are not given. Best to leave this TRUE!
This argument is set ot FALSE if warn is less than zero in the
top level, R options function. See options()$warn</p>
</td></tr>
<tr valign="top"><td><code>knots</code></td>
<td>
<p>A matrix of locations similar to x. These can define an alternative set of
basis functions for representing the estimate. One choice may be a
space-filling subset of the original x locations, thinning out the
design where locations cluster. The
default is to put a "knot" at all unique locations. (See details.)
</p>
</td></tr>
<tr valign="top"><td><code>lambda</code></td>
<td>
<p>Smoothing parameter that is the ratio of the error variance (sigma**2)
to the scale parameter of the
covariance function (rho). If omitted this is estimated by GCV ( see
method below).
</p>
</td></tr>
<tr valign="top"><td><code>method</code></td>
<td>
<p>Determines what "smoothing" parameter should be used. The default
is to estimate standard GCV
Other choices are: GCV.model, GCV.one, RMSE, pure error and REML. The
differences are explained below.
</p>
</td></tr>
<tr valign="top"><td><code>mean.obj</code></td>
<td>
<p>Object to predict the mean of the spatial process. This used in when
fitting a correlation model with varying spatial means and varying
marginal variances. (See details.)
</p>
</td></tr>
<tr valign="top"><td><code>m</code></td>
<td>
<p>A polynomial function of degree (m-1) will be
included in the model as the drift (or spatial trend) component.
The "m" notation is from thin-plate splines where m is the
derivative in the penalty function. With m=2 as the default a linear
model in the locations will be fit a fixed part of the model.
</p>
</td></tr>
<tr valign="top"><td><code>na.rm</code></td>
<td>
<p>If TRUE NAs will be removed from the <code>y</code> vector and the
corresponding rows of <code>x</code> – with a warning.
If FALSE Krig will just stop with a message. Once NAs are removed all
subsequent analysis in fields does not use those data. </p>
</td></tr>
<tr valign="top"><td><code>nstep.cv</code></td>
<td>
<p> Number of grid points for the coarse grid search to
minimize the GCV RMLE and other related criteria for finding lambda,
the smoothing parameter. Default is 200, fairly large to avoid some
cases of closely spaced local minima. Evaluations of the GCV and
related objective functions are cheap given the matrix decompositions
described below. </p>
</td></tr>
<tr valign="top"><td><code>null.args</code></td>
<td>
<p> Extra arguments for the null space function
<code>null.function</code>. If <code>fields.mkpoly</code> is passed as
<code>null.function</code> then this is set to a list with the value of
<code>m</code>. So the default is use a polynomial of degree m-1 for the
null space (fixed part) of the model. </p>
</td></tr>
<tr valign="top"><td><code>null.function</code></td>
<td>
<p>An R function that creates the matrices for the null space model.
The default is fields.mkpoly, an R function that creates a polynomial
regression matrix with all terms up to degree m-1. (See Details)
</p>
</td></tr>
<tr valign="top"><td><code>offset</code></td>
<td>
<p> The offset to be used in the GCV criterion. Default is
0. This would be used when Krig is part of a backfitting algorithm and
the offset is other model degrees of freedom from other regression
components. </p>
</td></tr>
<tr valign="top"><td><code>rho</code></td>
<td>
<p>Scale factor for covariance.
</p>
</td></tr>
<tr valign="top"><td><code>scale.type</code></td>
<td>
<p>This is a character string among: "range", "unit.sd", "user", "unscaled".
The independent variables and knots are scaled to the specified scale.type.
By default no scaling is done. This usuall makes sense for spatial locations.
Scale type of "range" scales the data to the interval (0,1) by forming
(x-min(x))/range(x) for each x. Scale type of "unit.sd"
Scale type of "user" allows specification of an x.center and x.scale by the
user. The default for "user" is mean 0 and standard deviation 1. Scale
type of "unscaled" does not scale the data.
</p>
</td></tr>
<tr valign="top"><td><code>sd.obj</code></td>
<td>
<p>Object to predict the marginal standard deviation of the spatial process.
</p>
</td></tr>
<tr valign="top"><td><code>sigma2</code></td>
<td>
<p>Variance of the errors, often called the nugget variance. If weights are
specified then the error variance is sigma2 divided by weights.
Note that lambda is defined as the ratio sigma2/rho.
</p>
</td></tr>
<tr valign="top"><td><code>verbose</code></td>
<td>
<p>If true will print out all kinds of intermediate stuff. Default is false,
of course as this is used mainly for debugging.
</p>
</td></tr>
<tr valign="top"><td><code>weights</code></td>
<td>
<p>Weights are proportional to the reciprocal variance of the measurement
error. The default is equal weighting i.e. vector of unit weights.
</p>
</td></tr>
<tr valign="top"><td><code>wght.function</code></td>
<td>
<p>An R function that creates a weights matrix to the observations.
This is only needed if the weight matirx has off diagonal elements.
The default is NULL indicating that the weight matrix is a diagonal, based
on the weights argument. (See details)
</p>
</td></tr>
<tr valign="top"><td><code>W</code></td>
<td>
<p>The observation weight matrix.</p>
</td></tr>
<tr valign="top"><td><code>wght.args</code></td>
<td>
<p>Optional arguments to be passed to the weight function (wght.function)
used to create the observation weight matrix.</p>
</td></tr>
<tr valign="top"><td><code>x</code></td>
<td>
<p>Matrix of independent variables. These could the locations for spatial
data or the indepedent variables in a regression.
</p>
</td></tr>
<tr valign="top"><td><code>x.center</code></td>
<td>
<p>Centering values to be subtracted from each column of the x matrix.
</p>
</td></tr>
<tr valign="top"><td><code>x.scale</code></td>
<td>
<p>Scale values that are divided into each column after centering.
</p>
</td></tr>
<tr valign="top"><td><code>Y</code></td>
<td>
<p>Vector of dependent variables. These are the values of the surface
(perhaps with measurement error) at the locations or the dependent
response in a regression.
</p>
</td></tr>
<tr valign="top"><td><code>Z</code></td>
<td>
<p>A vector of matrix of covariates to be include in the fixed part of the
model. If NULL (default) no addtional covariates are included.</p>
</td></tr>
<tr valign="top"><td><code>...</code></td>
<td>
<p>Optional arguments that appear are assumed to be additional arguments
to the covariance function. Or are included in methods functions (resid,
fitted, coef) as a
required argument.</p>
</td></tr>
<tr valign="top"><td><code>object</code></td>
<td>
<p> A Krig object</p>
</td></tr>
</table>
<h3>Details</h3>
<p>This function produces a object of class Krig. With this object it is
easy to subsequently predict with this fitted surface, find standard
errors, alter the y data ( but not x), etc.
</p>
<p>The Kriging model is: Y.k= f(x.k) = P(x.k) + Z(x.k) + e.k
</p>
<p>where ".k" means subscripted by k, Y is the dependent variable observed
at location x.k, P is a low order polynomial, Z is a mean zero, Gaussian
field with covariance function K and e is assumed to be independent
normal errors. The estimated surface is the best linear unbiased
estimate (BLUE) of f(x)= P(x) + Z(x) given the observed data. For this
estimate K, is taken to be rho*cov.function and the errors have variance
sigma**2. In more conventional geostatistical terms rho is the "sill" if
the covariance function is actually a correlation function and sigma**2
is the nugget variance or measure error variance (the two are confounded
in this model.) If the weights are given then the variance of e.k is
sigma**2/ weights.k . In the case that the weights are specified as a
matrix, W, using the wght.function option then the assumed covariance
matrix for the errors is sigma**2 Wi, where Wi is the inverse of W. It
is straightforward to show that the estimate of f only depends on sigma
and rho through the ratio lambda = sigma**2/ rho. This parameter, termed
the smoothing parameter plays a central role in the statistical
computations within <code>Krig</code>. See also the help for thin plate
splines, (<code>Tps</code>) to get another perspective on the smoothing
parameter.
</p>
<p>This function also supports a modest extension of the Generalized
Kriging model to include other covariates as fixed regression type
components. In matrix form Y = Zb + F + E where Z is a matrix of
covariates and b a fixed parameter vector, F the vector of function values at
the observations and E a vector of errors. The The <code>Z</code> argument in
the function is the way to specify this additional component.
</p>
<p>If the parameters rho and sigma2 are omitted in the call, then they are
estimated in the following way. If lambda is given, then sigma2 is
estimated from the residual sum of squares divided by the degrees of
freedom associated with the residuals. Rho is found as the difference
between the sums of squares of the predicted values having subtracted off
the polynomial part and sigma2. These estimates are the MLE's under Gaussian
assumptions on the process and errors.
If lambda is also omitted is it estimated
from the data using a variety of approaches and then the values for sigma
and rho are found in the same way from the estimated lambda.
</p>
<p>A useful extension of a stationary correlation to a nonstationary
covariance is what we term a correlation model.
If mean and marginal standard deviation objects are included in the call.
Then the observed data is standardized based on these functions. The
spatial process is then estimated with respect to the standardized scale.
However for predictions and standard errors the mean and standard
deviation surfaces are used to produce results in the original scale of
the observations.
</p>
<p>The GCV function has several alternative definitions when replicate
observations are present or if one uses a reduced set knots. Here are the
choices based on the method argument:
</p>
<p>GCV: leave-one-out GCV. But if
there are replicates it is leave one group out. (Wendy and Doug prefer
this one.)
</p>
<p>GCV.one: Really leave-one-out GCV even if there are replicate
points. This what the old tps function used in FUNFITS.
</p>
<p>rmse: Match the estimate of sigma**2 to a external value ( called rmse)
</p>
<p>pure error: Match the estimate of sigma**2 to the estimate based on
replicated data (pure error estimate in ANOVA language).
</p>
<p>GCV.model:
Only considers the residual sums of squares explained by the basis
functions.
</p>
<p>REML:
The process and errors are assumed to the Gaussian and the likelihood is
concentrated (or profiled) with respect to lambda. The MLE of lambda is
found from this criterion. Restricted means that the likelihood is formed from a linear transformation of the observations that is orthogonal to the
column space of P(x).
</p>
<p>WARNING: The covariance functions often have a nonlinear parameter(s) that
often control the strength of the correlations as a function of separation,
usually referred to as the range parameter. This parameter must be
specified in the call to Krig and will not be estimated.
</p>
<h3>Value</h3>
<p>A object of class Krig. This includes the predicted values in
fitted.values and the residuals in residuals. The results of the grid
search to minimize the generalized cross validation function are
returned in gcv.grid.
</p>
<p>The coef.Krig function only returns the coefficients, "d", associated with the
fixed part of the model (also known as the null space or spatial drift).
</p>
<table summary="R valueblock">
<tr valign="top"><td><code>call</code></td>
<td>
<p>Call to the function
</p>
</td></tr>
<tr valign="top"><td><code>y</code></td>
<td>
<p>Vector of dependent variables.
</p>
</td></tr>
<tr valign="top"><td><code>x</code></td>
<td>
<p>Matrix of independent variables.
</p>
</td></tr>
<tr valign="top"><td><code>weights</code></td>
<td>
<p>Vector of weights.
</p>
</td></tr>
<tr valign="top"><td><code>knots</code></td>
<td>
<p>Locations used to define the basis functions.
</p>
</td></tr>
<tr valign="top"><td><code>transform</code></td>
<td>
<p>List of components used in centering and scaling data.
</p>
</td></tr>
<tr valign="top"><td><code>np</code></td>
<td>
<p>Total number of parameters in the model.
</p>
</td></tr>
<tr valign="top"><td><code>nt</code></td>
<td>
<p>Number of parameters in the null space.
</p>
</td></tr>
<tr valign="top"><td><code>matrices</code></td>
<td>
<p>List of matrices from the decompositions (D, G, u, X, qr.T).
</p>
</td></tr>
<tr valign="top"><td><code>gcv.grid</code></td>
<td>
<p>Matrix of values from the GCV grid search. The first column
is the grid of lambda values used in the search, the second column
is the trace of the A matrix, the third column is the GCV values and
the fourth column is the estimated value of sigma conditional on the vlaue
of lambda.
</p>
</td></tr>
<tr valign="top"><td><code>lambda.est</code></td>
<td>
<p>A table of estimated smoothing parameters with corresponding degrees
of freedom and estimates of sigma found by different methods.
</p>
</td></tr>
<tr valign="top"><td><code>cost</code></td>
<td>
<p>Cost value used in GCV criterion.
</p>
</td></tr>
<tr valign="top"><td><code>m</code></td>
<td>
<p>Order of the polynomial space: highest degree polynomial is (m-1).
This is a fixed part of the surface often referred to as the drift
or spatial trend.
</p>
</td></tr>
<tr valign="top"><td><code>eff.df</code></td>
<td>
<p>Effective degrees of freedom of the model.
</p>
</td></tr>
<tr valign="top"><td><code>fitted.values</code></td>
<td>
<p>Predicted values from the fit.
</p>
</td></tr>
<tr valign="top"><td><code>residuals</code></td>
<td>
<p>Residuals from the fit.
</p>
</td></tr>
<tr valign="top"><td><code>lambda</code></td>
<td>
<p>Value of the smoothing parameter used in the fit.
Lambda is defined as sigma**2/rho. See discussion in details.
</p>
</td></tr>
<tr valign="top"><td><code>yname</code></td>
<td>
<p>Name of the response.
</p>
</td></tr>
<tr valign="top"><td><code>cov.function</code></td>
<td>
<p>Covariance function of the model.
</p>
</td></tr>
<tr valign="top"><td><code>beta</code></td>
<td>
<p>Estimated coefficients in the ridge regression format
</p>
</td></tr>
<tr valign="top"><td><code>d</code></td>
<td>
<p>Estimated coefficients for the polynomial basis functions that span the
null space
</p>
</td></tr>
<tr valign="top"><td><code>fitted.values.null</code></td>
<td>
<p>Fitted values for just the polynomial part of the estimate
</p>
</td></tr>
<tr valign="top"><td><code>trace</code></td>
<td>
<p>Effective number of parameters in model.
</p>
</td></tr>
<tr valign="top"><td><code>c</code></td>
<td>
<p>Estimated coefficients for the basis functions derived from the
covariance.
</p>
</td></tr>
<tr valign="top"><td><code>coefficients</code></td>
<td>
<p>Same as the beta vector.
</p>
</td></tr>
<tr valign="top"><td><code>just.solve</code></td>
<td>
<p>Logical describing if the data has been interpolated using the basis
functions.
</p>
</td></tr>
<tr valign="top"><td><code>shat</code></td>
<td>
<p>Estimated standard deviation of the measurement error (nugget effect).
</p>
</td></tr>
<tr valign="top"><td><code>sigma2</code></td>
<td>
<p>Estimated variance of the measurement error (shat**2).
</p>
</td></tr>
<tr valign="top"><td><code>rho</code></td>
<td>
<p>Scale factor for covariance. COV(h(x),h(x<code>)) = rho*cov.function(x,x</code>)
If the covariance is actually a
correlation function then rho is also the "sill".
</p>
</td></tr>
<tr valign="top"><td><code>mean.var</code></td>
<td>
<p>Normalization of the covariance function used to find rho.
</p>
</td></tr>
<tr valign="top"><td><code>best.model</code></td>
<td>
<p>Vector containing the value of lambda, the estimated variance of the
measurement error and the scale factor for covariance used in the fit.
</p>
</td></tr>
</table>
<h3>References</h3>
<p>See "Additive Models" by Hastie and Tibshirani, "Spatial Statistics" by
Cressie and the FIELDS manual.
</p>
<h3>See Also</h3>
<p>summary.Krig, predict.Krig, predictSE.Krig, predictSurfaceSE,
predictSurface, plot.Krig,
surface.Krig
</p>
<h3>Examples</h3>
<pre>
# a 2-d example
# fitting a surface to ozone
# measurements. Exponential covariance, range parameter is 20 (in miles)
fit <- Krig(ChicagoO3$x, ChicagoO3$y, theta=20)
summary( fit) # summary of fit
set.panel( 2,2)
plot(fit) # four diagnostic plots of fit
set.panel()
surface( fit, type="C") # look at the surface
# predict at data
predict( fit)
# predict using 7.5 effective degrees of freedom:
predict( fit, df=7.5)
# predict on a grid ( grid chosen here by defaults)
out<- predictSurface( fit)
surface( out, type="C") # option "C" our favorite
# predict at arbitrary points (10,-10) and (20, 15)
xnew<- rbind( c( 10, -10), c( 20, 15))
predict( fit, xnew)
# standard errors of prediction based on covariance model.
predictSE( fit, xnew)
# surface of standard errors on a default grid
predictSurfaceSE( fit)-> out.p # this takes some time!
surface( out.p, type="C")
points( fit$x)
## Not run:
# Using another stationary covariance.
# smoothness is the shape parameter for the Matern.
fit <- Krig(ChicagoO3$x, ChicagoO3$y, Covariance="Matern", theta=10, smoothness=1.0)
summary( fit)
#
# Roll your own: creating very simple user defined Gaussian covariance
#
test.cov <- function(x1,x2,theta,marginal=FALSE,C=NA){
# return marginal variance
if( marginal) { return(rep( 1, nrow( x1)))}
# find cross covariance matrix
temp<- exp(-(rdist(x1,x2)/theta)**2)
if( is.na(C[1])){
return( temp)}
else{
return( temp%*%C)}
}
#
# use this and put in quadratic polynomial fixed function
fit.flame<- Krig(flame$x, flame$y, cov.function="test.cov", m=3, theta=.5)
#
# note how range parameter is passed to Krig.
# BTW: GCV indicates an interpolating model (nugget variance is zero)
# This is the content of the warning message.
# take a look ...
surface(fit.flame, type="I")
## End(Not run)
#
# Thin plate spline fit to ozone data using the radial
# basis function as a generalized covariance function
#
# p=2 is the power in the radial basis function (with a log term added for
# even dimensions)
# If m is the degree of derivative in penalty then p=2m-d
# where d is the dimension of x. p must be greater than 0.
# In the example below p = 2*2 - 2 = 2
#
out<- Krig( ChicagoO3$x, ChicagoO3$y,cov.function="Rad.cov",
m=2,p=2,scale.type="range")
# See also the Fields function Tps
# out should be identical to Tps( ChicagoO3$x, ChicagoO3$y)
#
# A Knot example
data(ozone2)
y16<- ozone2$y[16,]
# there are some missing values -- remove them
good<- !is.na( y16)
y<- y16[good]
x<- ozone2$lon.lat[ good,]
#
# the knots can be arbitrary but just for fun find them with a space
# filling design. Here we select 50 from the full set of 147 points
#
xknots<- cover.design( x, 50, num.nn= 75)$design # select 50 knot points
out<- Krig( x, y, knots=xknots, cov.function="Exp.cov", theta=300)
summary( out)
# note that that trA found by GCV is around 17 so 50>17 knots may be a
# reasonable approximation to the full estimator.
#
## Not run:
# the plot
surface( out, type="C")
US( add=TRUE)
points( x, col=2)
points( xknots, cex=2, pch="O")
## End(Not run)
## Not run:
## A quick way to deal with too much data if you intend to smooth it anyway
## Discretize the locations to a grid, then apply Krig
## to the discretized locations:
##
RM.approx<- as.image(RMprecip$y, x=RMprecip$x, nx=20, ny=20)
# take a look:
image.plot( RM.approx)
# discretized data (observations averaged if in the same grid box)
# 336 locations -- down form the full 806
# convert the image format to locations, obs and weight vectors
yd<- RM.approx$z[RM.approx$ind]
weights<- RM.approx$weights[RM.approx$ind] # takes into account averaging
xd<- RM.approx$xd
obj<- Krig( xd, yd, weights=weights, theta=4)
# compare to the full fit:
# Krig( RMprecip$x, RMprecip$y, theta=4)
## End(Not run)
## Not run:
# A correlation model example
# fit krig surface using a mean and sd function to standardize
# first get stats from 1987 summer Midwest O3 data set
data(ozone2)
stats.o3<- stats( ozone2$y)
mean.o3<- Tps( ozone2$lon.lat, c( stats.o3[2,]))
sd.o3<- Tps( ozone2$lon.lat, c( stats.o3[3,]))
#
# Now use these to fit particular day ( day 16)
# and use great circle distance
fit<- Krig( ozone2$lon.lat, ozone2$y[16,],
theta=350, mean.obj=mean.o3, sd.obj=sd.o3,
Covariance="Matern", Distance="rdist.earth",
smoothness=1.0,
na.rm=TRUE) #
# the finale
surface( fit, type="I")
US( add=TRUE)
points( fit$x)
title("Estimated ozone surface")
## End(Not run)
## Not run:
#
#
# explore some different values for the range and lambda using REML
theta <- seq( 100,500,,40)
PLL<- matrix( NA, 40,80)
# the loop
for( k in 1:40){
# call to Krig with different ranges
# also turn off warnings for GCV search
# to avoid lots of messages. (not recommended in general!)
PLL[k,]<- Krig( ozone2$lon.lat,ozone2$y[16,],
cov.function="stationary.cov",
theta=theta[k], mean.obj=mean.o3, sd.obj=sd.o3,
Covariance="Matern",smoothness=.5,
Distance="rdist.earth", nstep.cv=80,
give.warnings=FALSE, na.rm=TRUE)$gcv.grid[,7]
#
# gcv.grid is the grid search output from
# the optimization for estimating different estimates for lambda including
# REML
# default grid is equally spaced in eff.df scale ( and should the same across theta)
# here
}
# get lambda grid from looping
k<- 1
lam<- Krig( ozone2$lon.lat,ozone2$y[16,],
cov.function="stationary.cov",
theta=theta[k], mean.obj=mean.o3, sd.obj=sd.o3,
Covariance="Matern",smoothness=.5,
Distance="rdist.earth", nstep.cv=80,
give.warnings=FALSE, na.rm=TRUE)$gcv.grid[,1]
# see the 2 column of $gcv.grid to get the effective degress of freedom.
contour( theta,log(lam) , PLL)
## 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>