tests/Krig.test.R

In fields: Tools for Spatial Data

#
# 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))
#
#
#  test of fixed lambda case
#  Check against linear algebra
#

options( echo=FALSE)
test.for.zero.flag<-1

fit<- Krig( ChicagoO3$x, ChicagoO3$y, aRange=50)

x<- ChicagoO3$x
K<- Exp.cov(x, x,aRange=50)
T<- fields.mkpoly(x, 2)
W<- diag( 20)
 lambda<- fit$lambda
M<- (lambda* diag(20) + K) 
###########################
test.d<- c(solve( t(T) %*% solve( M)%*%T) %*% t(T)%*% solve( M) %*% fit$yM)
test.c<- solve( M)%*% ( fit$yM - T%*% test.d)

#compare to  fit$d
test.for.zero( test.d, fit$d, tag="Compare d coef" )
#compare to  fit$d
test.for.zero( test.c, fit$c,tag="Compare c coef" )

fit2<- Krig( ChicagoO3$x, ChicagoO3$y, aRange=50, lambda= fit$lambda)
#compare to  fit$d
test.for.zero( test.d, fit2$d, tag="Compare d coef fixed lambda" )
#compare to  fit$d
test.for.zero( test.c, fit2$c,tag="Compare c coef fixed lambda" )

# test of Krig.coef

Krig.coef( fit)->test
test.for.zero( test.d, test$d, tag="d coef Krig.coef" )
test.for.zero( test.c, test$c, tag= "c coef Krig.coef" )

Krig.coef( fit2)->test
test.for.zero( test.d, test$d,tag="d coef Krig.coef fixed" )
test.for.zero( test.c, test$c, tag="c coef Krig.coef fixed" )
# checking A matrix in the case of noreps

set.seed( 222)
weights<-  10+ runif( length(ChicagoO3$y))
#weights<- rep( 1, 20)
test2<- Krig( ChicagoO3$x, ChicagoO3$y, aRange=50, weights= weights)
Atest<- Krig.Amatrix( test2)
K<-Exp.cov(ChicagoO3$x, ChicagoO3$x,aRange=50)
H<- matrix(0, 23,23)
H[(1:20)+3 , (1:20)+3]<- K
X<- cbind( fields.mkpoly( ChicagoO3$x, 2), K)
lambda<- test2$lambda
 Alam <-  X%*%solve(
                 t(X)%*%diag(weights)%*%X + lambda*H
                 )%*% t(X)%*%diag(weights) 
 test.for.zero( Alam, Atest, tag="Amatrix no reps", tol=5e-8)

# test for new y fixed case
set.seed( 123)
ynew<- rnorm( fit2$N)

test.d<- c(solve( t(T) %*% solve( M)%*%T) %*% t(T)%*% solve( M) %*% ynew)
test.c<- solve( M)%*% ( ynew - T%*% test.d)

test<- Krig.coef( fit, y= ynew)
test.for.zero( test.d, test$d, tag= "d coef new y" )
test.for.zero( test.c, test$c, tag="c coef new y" )


Krig.coef( fit2, y= ynew)->test
test.for.zero( test.d, test$d, tag= "d coef new y fixed" )
test.for.zero( test.c, test$c, tag=" c coef new y fixed"  )

# test for multiple new y's
Krig.coef( fit2, y= cbind( ynew+ rnorm(fit2$N), ynew))->test2
test.for.zero( test.d, test2$d[,2], tag= "d coef several new y fixed" )
test.for.zero( test.c, test2$c[,2], tag=" c coef several new y fixed"  )


#cat("done with simple Krig data", fill=TRUE)


# These tests are about whether decompositions 
# handle just a fixed lambda or are more general 

# checking passing lambda or df to Krig

out<- Tps( ChicagoO3$x, ChicagoO3$y,lambda=.001 )
out2<- predict( out, lambda=.001)
test.for.zero( out2, predict( out), tag="Tps with fixed lam")

out<- Tps( ChicagoO3$x, ChicagoO3$y, df=5)
out2<- predict( out, df=5)
test.for.zero( out2, predict( out), tag="Tps with fixed df")

# same for Krig

out0<- Krig( ChicagoO3$x, ChicagoO3$y, aRange=50,lambda=.5)
out<- Krig( ChicagoO3$x, ChicagoO3$y, aRange=50,lambda=.5,GCV=TRUE)
test.for.zero( 
      predict(out0), predict( out), tag="Krig with fixed lam argument")

#A very nasty case with knots and weights

set.seed(123)
x<- matrix( runif( 30), 15,2)
y<- rnorm( 15)*.01 + x[,1]**2 +  x[,2]**2

weights<- runif(15)*10

# compare to 
Krig( x,y, cov.function=Exp.cov, weights=weights)-> out.new
Krig( x,y, cov.function=Exp.cov, weights=weights, 
          lambda=1)-> out.new2

# compute test using linear algebra
K<- Exp.cov( x, x)
H<- matrix(0, 18,18)
H[4:18, 4:18]<- K
X<- cbind( fields.mkpoly( x, 2), Exp.cov( x, x))
lambda<-1


c(   solve(t(X)%*%(weights*X) + lambda*H)%*% t(X)%*% (weights*y) )-> temp
temp.c<- temp[4:18]
temp.d<- temp[1:3]


# test for d coefficients
test.for.zero( out.new2$d, temp.d, tag=" d coef")
# test for c coefficents
test.for.zero( out.new2$c, temp.c, tag="c coef" )

# and 
test<- Krig.coef( out.new2, lambda=1)

# test for d coefficients
test.for.zero( temp.d, test$d, tag= "d fixed case")
# test for c coefficents 
test.for.zero( temp.c, test$c, tag=" c fixed case" )


ynew<- 1:15

#compare 
test<- Krig.coef( out.new, lambda=.5, y=ynew) 
test2<- Krig( x,ynew, cov.function=Exp.cov,
              lambda= .5, weights=weights)
# test for d coefficients
test.for.zero(  test2$d,test$d, tag=" d new y")
# test for c coefficents 
test.for.zero( test2$c, test$c,tag= "c new y" )




#cat("test with reps" , fill=TRUE)
#

##################################
#cat( "test  A matrix",fill=TRUE)
##################################
 
set.seed(133)
x<- matrix( runif( 30), 15,2)*2  
y<- rnorm( nrow( x))*.5 + + x[,1]**2 +  x[,2]**2
# perturb so that this example does not generate (harmless) warnings in gcv search
weights<- runif( nrow( x))*10

out.new<- Krig( x,y, weights= weights)


testY<- predict( out.new)
testY2<-  Krig.Amatrix(out.new)%*% y
test.for.zero( testY, testY2, tag="testing A matrix")

set.seed(333)
yNew<- rnorm( 15)
testY<- predict( out.new, y=yNew)
testY2<-  Krig.Amatrix(out.new)%*% yNew
test.for.zero( testY, testY2, tag="testing A matrix new y")

Alam<- Krig.Amatrix(out.new)
trA<- sum( diag( Alam))
test.for.zero( trA, out.new$eff.df,  tag="checking trace")

####### checking GCV
#######

MSE<- mean( out.new$residuals^2 * weights)
n<- length( y)
GCV<- MSE/ (1- trA/n )^2
 
test.for.zero(GCV, out.new$lambda.est[6,3], tag="GCV" )