KrigFindLambda.RdThis is a secondary function that will use the computed Krig object and find various estimates of the smoothing parameter lambda. These are several different flavors of cross-validation, a moment matching strategy and the profile likelihood. This function can also be used independently with different data sets (the y's) if the covariates ( the x's) are the same and thus reduce the computation.
KrigFindLambda(
out, lambda.grid = NA, cost = 1, nstep.cv = 200, rmse
= NA, verbose = FALSE, tol = 1e-05, offset = 0, y =
NULL, give.warnings = TRUE)
gcv.sreg (
out, lambda.grid = NA, cost = 1, nstep.cv = 80, rmse =
NA, offset = 0, trmin = NA, trmax = NA, verbose =
FALSE, tol = 1e-05, give.warnings = TRUE)A Krig or sreg object.
Grid of lambdas for coarse search. The default is equally spaced on effective degree of freedom scale.
Cost used in GCV denominator
Number of grid points in coarse search.
Target root mean squared error to match with the estimate of tau**2
If true prints intermediate results.
Tolerance in delcaring convergence of golden section search or bisection search.
Additional degrees of freedom to be added into the GCV denominator.
A new data vector to be used in place of the one associated with the Krig object (obj)
If FALSE will suppress warnings about grid search being out of range for various estimates based on GCV and REML.
Minimum value of lambda for grid search specified in terms of effective degrees of freedom.
Maximum value for grid search.
This function finds several estimates of the smoothing parameter using first a coarse grid search followed by a refinement using a minimization ( in the case of GCV or maximum likelihood) or bisection in the case of mathcing the rmse. Details of the estimators can be found in the help file for the Krig function.
The Krig object passed to this function has some matrix decompostions that facilitate rapid computation of the GCV and ML functions and do not depend on the independent variable. This makes it possible to compute the Krig object once and to reuse the decompostions for multiple data sets. (But keep in mind if the x values change then the object must be recalculated.) The example below show show this can be used for a simulation study on the variability for estimating the smoothing parameter.
A list giving a summary of estimates and diagonostic details with the following components:
A matrix describing results of the coarse search rows are values of lambda and the columns are lambda= value of smoothing parameter, trA=effective degrees of freedom, GCV=Usual GCV criterion, GCV.one=GCV criterion leave-one-out, GCV.model= GCV based on average response in the case of replicates, tauHat= Implied estimate of tau , -Log Profile= negative log of profiel likelihood for the lambda.
Summary table of all estimates Rows index different types of estimates: GCV, GCV.model, GCV.one, RMSE, pure error, -Log Profile and the columns are the estimated values for lambda, trA, GCV, tauHat.
#
Tps( ChicagoO3$x, ChicagoO3$y)-> obj # default is to find lambda by GCV
summary( obj)
#> CALL:
#> Tps(x = ChicagoO3$x, Y = ChicagoO3$y)
#>
#> Number of Observations: 20
#> Number of unique points: 20
#> Number of parameters in the null space 3
#> Parameters for fixed spatial drift 3
#> Effective degrees of freedom: 4.5
#> Residual degrees of freedom: 15.5
#> MLE tau 3.779
#> GCV tau 4.073
#> MLE sigma 347.7
#> Scale passed for covariance (sigma) <NA>
#> Scale passed for nugget (tau^2) <NA>
#> Smoothing parameter lambda 0.04107
#>
#> Residual Summary:
#> min 1st Q median 3rd Q max
#> -6.8060 -1.4390 -0.5064 1.4440 7.7890
#>
#> Covariance Model: Rad.cov
#> Names of non-default covariance arguments:
#> p
#>
#> DETAILS ON SMOOTHING PARAMETER:
#> Method used: GCV Cost: 1
#> lambda trA GCV GCV.one GCV.model tauHat
#> 0.04107 4.50304 21.40938 21.40938 NA 4.07296
#>
#> Summary of all estimates found for lambda
#> lambda trA GCV tauHat -lnLike Prof converge
#> GCV 0.04107 4.503 21.41 4.073 49.00 5
#> GCV.model NA NA NA NA NA NA
#> GCV.one 0.04107 4.503 21.41 4.073 NA 5
#> RMSE NA NA NA NA NA NA
#> pure error NA NA NA NA NA NA
#> REML 0.02972 4.886 21.49 4.030 48.98 4
KrigFindLambda( obj)-> out
print( out$lambda.est) # results agree with Tps summary
#> lambda trA GCV tauHat -lnLike Prof converge
#> GCV 0.04107313 4.503038 21.40938 4.072962 49.00042 5
#> GCV.model NA NA NA NA NA NA
#> GCV.one 0.04107313 4.503038 21.40938 4.072962 NA 5
#> RMSE NA NA NA NA NA NA
#> pure error NA NA NA NA NA NA
#> REML 0.02971891 4.886273 21.48837 4.029698 48.98323 4
sreg( rat.diet$t, rat.diet$trt)-> out
gcv.sreg( out, tol=1e-10) # higher tolerance search for minimum
#> $gcv.grid
#> lambda trA GCV GCV.one GCV.model tauHat
#> 1 8.183514e-04 37.049991 8.075493 8.075493 NA 0.6354341
#> 2 1.028159e-03 36.645399 8.050161 8.050161 NA 0.6971539
#> 3 1.291758e-03 36.177994 8.017769 8.017769 NA 0.7616816
#> 4 1.622937e-03 35.645737 7.976059 7.976059 NA 0.8282481
#> 5 2.039024e-03 35.048882 7.922064 7.922064 NA 0.8958739
#> 6 2.561787e-03 34.390175 7.852044 7.852044 NA 0.9633881
#> 7 3.218575e-03 33.674681 7.761566 7.761566 NA 1.0294736
#> 8 4.043749e-03 32.909227 7.645807 7.645807 NA 1.0927368
#> 9 5.080481e-03 32.101546 7.500139 7.500139 NA 1.1518031
#> 10 6.383009e-03 31.259303 7.320957 7.320957 NA 1.2054290
#> 11 8.019478e-03 30.389253 7.106600 7.106600 NA 1.2526192
#> 12 1.007550e-02 29.496749 6.858138 6.858138 NA 1.2927272
#> 13 1.265865e-02 28.585697 6.579728 6.579728 NA 1.3255215
#> 14 1.590406e-02 27.658936 6.278397 6.278397 NA 1.3511980
#> 15 1.998152e-02 26.718887 5.963220 5.963220 NA 1.3703358
#> 16 2.510436e-02 25.768238 5.644109 5.644109 NA 1.3838029
#> 17 3.154059e-02 24.810474 5.330515 5.330515 NA 1.3926315
#> 18 3.962693e-02 23.850097 5.030364 5.030364 NA 1.3978879
#> 19 4.978644e-02 22.892527 4.749417 4.749417 NA 1.4005594
#> 20 6.255063e-02 21.943732 4.491090 4.491090 NA 1.4014757
#> 21 7.858729e-02 21.009719 4.256659 4.256659 NA 1.4012676
#> 22 9.873541e-02 20.096021 4.045690 4.045690 NA 1.4003629
#> 23 1.240491e-01 19.207263 3.856567 3.856567 NA 1.3990108
#> 24 1.558526e-01 18.346901 3.687000 3.687000 NA 1.3973233
#> 25 1.958100e-01 17.517133 3.534459 3.534459 NA 1.3953247
#> 26 2.460115e-01 16.718988 3.396497 3.396497 NA 1.3929989
#> 27 3.090837e-01 15.952539 3.270965 3.270965 NA 1.3903276
#> 28 3.883262e-01 15.217195 3.156109 3.156109 NA 1.3873155
#> 29 4.878848e-01 14.511992 3.050601 3.050601 NA 1.3840034
#> 30 6.129682e-01 13.835842 2.953499 2.953499 NA 1.3804710
#> 31 7.701203e-01 13.187710 2.864191 2.864191 NA 1.3768351
#> 32 9.675628e-01 12.566692 2.782335 2.782335 NA 1.3732454
#> 33 1.215625e+00 11.972029 2.707798 2.707798 NA 1.3698802
#> 34 1.527286e+00 11.403079 2.640610 2.640610 NA 1.3669421
#> 35 1.918850e+00 10.859261 2.580903 2.580903 NA 1.3646500
#> 36 2.410803e+00 10.340010 2.528864 2.528864 NA 1.3632278
#> 37 3.028882e+00 9.844737 2.484663 2.484663 NA 1.3628873
#> 38 3.805423e+00 9.372811 2.448394 2.448394 NA 1.3638090
#> 39 4.781053e+00 8.923550 2.420010 2.420010 NA 1.3661223
#> 40 6.006814e+00 8.496230 2.399287 2.399287 NA 1.3698897
#> 41 7.546834e+00 8.090092 2.385799 2.385799 NA 1.3750976
#> 42 9.481682e+00 7.704359 2.378923 2.378923 NA 1.3816557
#> 43 1.191259e+01 7.338243 2.377862 2.377862 NA 1.3894041
#> 44 1.496672e+01 6.990946 2.381688 2.381688 NA 1.3981267
#> 45 1.880387e+01 6.661671 2.389389 2.389389 NA 1.4075698
#> 46 2.362479e+01 6.349617 2.399931 2.399931 NA 1.4174613
#> 47 2.968168e+01 6.053988 2.412314 2.412314 NA 1.4275325
#> 48 3.729144e+01 5.774001 2.425636 2.425636 NA 1.4375387
#> 49 4.685218e+01 5.508892 2.439160 2.439160 NA 1.4472801
#> 50 5.886408e+01 5.257928 2.452380 2.452380 NA 1.4566240
#> 51 7.395559e+01 5.020415 2.465096 2.465096 NA 1.4655265
#> 52 9.291624e+01 4.795705 2.477490 2.477490 NA 1.4740559
#> 53 1.167380e+02 4.583193 2.490207 2.490207 NA 1.4824182
#> 54 1.466672e+02 4.382320 2.504464 2.504464 NA 1.4909878
#> 55 1.842695e+02 4.192570 2.522173 2.522173 NA 1.5003448
#> 56 2.315123e+02 4.013463 2.546106 2.546106 NA 1.5113200
#> 57 2.908672e+02 3.844543 2.580099 2.580099 NA 1.5250435
#> 58 3.654394e+02 3.685346 2.629251 2.629251 NA 1.5429833
#> 59 4.591303e+02 3.535377 2.700100 2.700100 NA 1.5669507
#> 60 5.768417e+02 3.394090 2.800645 2.800645 NA 1.5990342
#> 61 7.247317e+02 3.260896 2.940108 2.940108 NA 1.6414254
#> 62 9.105376e+02 3.135210 3.128290 3.128290 NA 1.6961152
#> 63 1.143980e+03 3.016522 3.374422 3.374422 NA 1.7644891
#> 64 1.437273e+03 2.904481 3.685565 3.685565 NA 1.8469127
#> 65 1.805759e+03 2.798959 4.064815 4.064815 NA 1.9424445
#> 66 2.268717e+03 2.700060 4.509777 4.509777 NA 2.0487936
#> 67 2.850368e+03 2.608079 5.011838 5.011838 NA 2.1625627
#> 68 3.581143e+03 2.523411 5.556621 5.556621 NA 2.2797131
#> 69 4.499272e+03 2.446432 6.125644 6.125644 NA 2.3961195
#> 70 5.652790e+03 2.377396 6.698779 6.698779 NA 2.5080728
#> 71 7.102046e+03 2.316354 7.256855 7.256855 NA 2.6126316
#> 72 8.922860e+03 2.263129 7.783785 7.783785 NA 2.7077851
#> 73 1.121049e+04 2.217331 8.267813 8.267813 NA 2.7924451
#> 74 1.408463e+04 2.178399 8.701848 8.701848 NA 2.8663206
#> 75 1.769563e+04 2.145660 9.083044 9.083044 NA 2.9297309
#> 76 2.223241e+04 2.118387 9.411967 9.411967 NA 2.9834094
#> 77 2.793233e+04 2.095851 9.691596 9.691596 NA 3.0283282
#> 78 3.509359e+04 2.077357 9.926394 9.926394 NA 3.0655601
#> 79 4.409085e+04 2.062267 10.121545 10.121545 NA 3.0961801
#> 80 5.539481e+04 2.050011 10.282391 10.282391 NA 3.1212022
#>
#> $lambda.est
#> lambda trA GCV tauHat converge
#> GCV 11.10732 7.44848 2.377623 1.386913 16
#> GCV.model NA NA NA NA NA
#> GCV.one 11.10732 7.44848 2.377623 1.386913 16
#> RMSE NA NA NA NA NA
#> pure error NA NA NA NA NA
#>
#> $warningTable
#> Warning Refine indexMIN leftEndpoint rightEndpoint lambda effdf
#> GCV FALSE TRUE 43 FALSE FALSE 11.91259 7.338243
#> GCV.model FALSE FALSE NA NA NA NA NA
#> GCV.one FALSE TRUE 43 FALSE FALSE 11.91259 7.338243
#> RMSE FALSE FALSE NA NA NA NA NA
#> pure error FALSE FALSE NA NA NA NA NA
#>
if (FALSE) {
# a simulation example
x<- seq( 0,1,,150)
f<- x**2*( 1-x)
f<- f/sqrt( var( f))
set.seed(123) # let's all use the same seed
tau<- .1
y<- f + rnorm( 150)*tau
Tps( x,y)-> obj # create Krig object
hold<- hold2<- matrix( NA, ncol=6, nrow=200)
for( k in 1:200){
# look at GCV estimates of lambda
# new data simulated
y<- f + rnorm(150)*tau
# save GCV estimates
lambdaTable<- KrigFindLambda(obj, y=y, give.warnings=FALSE)$lambda.est
hold[k,]<- lambdaTable[1,]
hold2[k,]<- lambdaTable[6,]
}
matplot( cbind(hold[,2], hold2[,2]),cbind( hold[,4],hold2[,4]),
xlab="estimated eff. df", ylab="tau hat", pch=16, col=c("orange3", "green2"), type="p")
yline( tau, col="grey", lwd=2)
}