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
|
#
# fields is a package for analysis of spatial data written for
# the R software environment.
# Copyright (C) 2022 Colorado School of Mines
# 1500 Illinois St., Golden, CO 80401
# Contact: Douglas Nychka, douglasnychka@gmail.edu,
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with the R software environment if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
# or see http://www.r-project.org/Licenses/GPL-2
##END HEADER
##END HEADER
suppressMessages(library(fields))
options( echo=FALSE)
set.seed( 123)
x<- matrix( runif( 20),10, 2)
y<- rnorm(10)
lambda<- .1
aRange<- .2
out<- mKrig( x,y, aRange= aRange, lambda=lambda)
test.for.zero( out$lnDetOmega,
2*log( prod(diag(chol(out$Omega))))
)
Mc<- exp( -rdist( x,x)/aRange) + lambda* diag( 1,10)
OmegaTest<- solve(t(out$Tmatrix)%*%solve( Mc)%*% out$Tmatrix)
test.for.zero( OmegaTest, out$Omega,tag= "mKrigOmega")
test.for.zero( log(det(OmegaTest)), out$lnDetOmega,
tag="lnDetOmega")
test.for.zero( log( det( Mc)), out$lnDetCov, tag="lnDetMc" )
# check that det adjustment really works.
set.seed( 323)
x<- matrix( runif( 20), 10, 2)
temp<- matrix( NA, 50,8)
aRangeGrid<- seq( .1,.5, ,50)
lambdaGrid<- 10**(runif( 50, -2,0))
Q<- qr.qy( qr( cbind( rep(1,10),x) ), diag( 1,10))
Q2<- Q[,4:10]
y<- rnorm(10)
testDet<- function(lambda, obj)
{
D2 <- obj$matrices$D[obj$matrices$D > 0]
u2 <- obj$matrices$u[obj$matrices$D > 0]
lD <- D2 * lambda
N2 <- length(D2)
sigma.MLE <- (sum((D2 * (u2)^2)/(1 + lD)))/N2
lnDetCov <- -sum(log(D2/(1 + lD)))
# -1 * (-N2/2 - log(2 * pi) * (N2/2) - (N2/2) * log(sigma.MLE) -
# (1/2) * lnDetCov)
return( c(lnDetCov, sigma.MLE) )
}
for ( k in 1:50) {
out<- mKrig( x,y, aRange = aRangeGrid[k],
lambda = lambdaGrid[k]
)
# turn off warnings for lambda search because all we want are
# matrix decompositions independent of lambda
out2<- Krig( x,y, aRange= aRangeGrid[k],
cov.args=list( Covariance = "Exponential"),
give.warnings=FALSE)
Mc<- exp( -rdist( x,x)/aRangeGrid[k] ) + lambdaGrid[k]* diag( 1,10)
X<- out$Tmatrix
temp[k,]<-c(
out$lnDetCov,
out$lnDetOmega,
log( det( solve(t( Q2)%*%Mc%*%Q2) ) ),
log( det(Mc) ),
-1*log( det( t(X)%*%solve(Mc)%*%X ) ),
testDet( lambdaGrid[k], out2 ),
out$summary["sigma2"]
)
}
test.for.zero( temp[,2], temp[,5], tag="testing det Omega formula")
resid<- temp[,1] - temp[,2] + temp[,3]
test.for.zero( mean(resid), resid, relative=FALSE,
tag="REML Det shortcut")
#### testing Krig verses mKrig
#
test.for.zero( temp[,3], -temp[,6],
tag="Q2 Det and Eigen Det")
###### testing sigma.MLE from mKrig and Krig
test.for.zero( (7/10)*temp[,7], temp[,8],
tag="sigma.MLE Krig verses mKrig")
|
