tests/Tps.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))
options(echo=FALSE)
test.for.zero.flag<- 1

data(ozone2)

x<- ozone2$lon.lat
y<- ozone2$y[16,]

temp<- Rad.cov( x,x, p=2)
temp2<- RadialBasis( rdist( x,x), M=2, dimension=2)

temp3<-  rdist( x,x)
temp3 <- ifelse( abs(temp3) < 1e-14, 0,log( temp3)*(temp3^2) )
temp3<- radbas.constant( 2,2)*temp3

test.for.zero( temp, temp2, tag="Tps radial basis function 2d")
test.for.zero( temp, temp3, tag="Tps radial basis function 2d")
test.for.zero( temp2,temp3, tag="Tps radial basis function 2d")


set.seed( 123)
xtemp<- matrix( runif( 40*3), ncol=3) 
temp<- Rad.cov( xtemp,xtemp, p= 2*4-3)
temp2<- RadialBasis( rdist( xtemp,xtemp), M=4, dimension=3)

temp3<-  rdist( xtemp,xtemp)
temp3 <- ifelse( abs(temp3) < 1e-14, 0, temp3^(2*4 -3) )
temp3<- radbas.constant( 4,3)*temp3

test.for.zero( temp, temp2, tag="Tps radial basis function 3d")
test.for.zero( temp, temp3, tag="Tps radial basis function 3d")
test.for.zero( temp2,temp3, tag="Tps radial basis function 3d")

#### testing multiplication of a vector
#### mainly to make the FORTRAN has been written correctly
#### after replacing the ddot call with an explicit  do loop
set.seed( 123)
C<- matrix( rnorm( 10*5),10,5 )
x<- matrix( runif( 10*2), 10,2)
temp3<- rdist( x,x)
K<-  ifelse( abs(temp3) < 1e-14, 0,log( temp3)*(temp3^2) )
K<- K * radbas.constant( 2,2)
test.for.zero( Rad.cov( x,x,m=2, C=C) , K%*%C, tol=1e-10)

set.seed( 123)
C<- matrix( rnorm( 10*5),10,5 )
x<- matrix( runif( 10*3), 10,3)
temp3<- rdist( x,x)
K<-  ifelse( abs(temp3) < 1e-14, 0,(temp3^(2*4-3)) )
K<- K * radbas.constant( 4,3)
test.for.zero( Rad.cov( x,x,m=4, C=C) , K%*%C,tol=1e-10)


#####  testing derivative formula

set.seed( 123)
C<- matrix( rnorm( 10*1),10,1 )
x<- matrix( runif( 10*2), 10,2)
temp0<-  Rad.cov( x,x, p=4, derivative=1, C=C)

eps<- 1e-6
temp1<- (
           Rad.cov( cbind(x[,1]+eps, x[,2]),x, p=4, derivative=0, C=C) 
         - Rad.cov( cbind(x[,1]-eps, x[,2]),x, p=4, derivative=0, C=C) )/ (2*eps)
temp2<- (
           Rad.cov( cbind(x[,1], x[,2]+eps),x, p=4, derivative=0, C=C) 
         - Rad.cov( cbind(x[,1], x[,2]-eps),x , p=4,derivative=0,C=C) )/ (2*eps)

test.for.zero( temp0[,1], temp1, tag=" der of Rad.cov", tol=1e-6)
test.for.zero( temp0[,2], temp2, tag=" der of Rad.cov", tol=1e-6)



# comparing  Rad.cov used by Tps with simpler function called 
# by stationary.cov
set.seed( 222)
x<- matrix( runif( 10*2), 10,2)
C<- matrix( rnorm( 10*3),10,3 ) 
temp<- Rad.cov( x,x, p=2, C=C)
temp2<- RadialBasis( rdist( x,x), M=2, dimension=2)%*%C
test.for.zero( temp, temp2)

#### Basic matrix form for Tps as sanity check
data("ozone2")
s<- ozone2$lon.lat
y<- ozone2$y[16,]

good<- !is.na( y)
s<- s[good,]
y<- y[good]
data(ozone2)


obj<-Tps( s,y, scale.type="unscaled", with.constant=FALSE)

# now work out the matrix expressions explicitly
lam.test<- obj$lambda
N<-length(y)

Tmatrix<- cbind( rep( 1,N), s)
D<- rdist( s,s)
R<- ifelse( D==0, 0, D**2 * log(D))
A<- rbind(
          cbind( R+diag(lam.test,N), Tmatrix),
          cbind( t(Tmatrix), matrix(0,3,3)))

 hold<-solve( A, c( y, rep(0,3)))
 c.coef<- hold[1:N]
 d.coef<- hold[ (1:3)+N]
 zhat<-  R%*%c.coef + Tmatrix%*% d.coef
  test.for.zero( zhat, obj$fitted.values, tag="Tps 2-d m=2 sanity check")
# out of sample prediction
snew<- rbind( c( -87,41),
              c( -81,44)
              )
T1<- cbind(  1, snew)
D<- rdist( snew,s)
R1<- ifelse( D==0, 0, D**2 * log(D))
z1<-  R1%*%c.coef + T1%*% d.coef
  test.for.zero( z1, predict( obj, x=snew), tag="Tps 2-d m=2 sanity predict")

#### test Tps verses Krig note scaling must be the same
   out<- Tps( s,y)
   out2<- Krig( s,y, Covariance="RadialBasis", 
           M=2, dimension=2, scale.type="range", method="GCV")
   test.for.zero( predict(out), predict(out2), tag="Tps vs.  Krig w/ GCV")

# test for fixed lambda
   test.for.zero( 
   predict(out,lambda=.1), predict(out2, lambda=.1),
   tag="Tps vs. radial basis w Krig")

#### testing derivative using predict function 
   set.seed( 233)
   x<- matrix( (rnorm( 1000)*2 -1), ncol=2)
   y<- (x[,1]**2 + 2*x[,1]*x[,2] -  x[,2]**2)/2

   out<- Tps( x, y, scale.type="unscaled")

   xg<- make.surface.grid( list(x=seq(-.7,.7,,10),  y=seq(-.7,.7,,10)) )
   test<- cbind( xg[,1] + xg[,2], xg[,1] - xg[,2])
#   test<- xg
   look<- predictDerivative.Krig( out, x= xg) 
   test.for.zero( look[,1], test[,1], tol=1e-3)
   test.for.zero( look[,2], test[,2], tol=1e-3)

############################################################
### testing Tps version using spatialProcess and Tps.cov
############################################################
   
   set.seed(222)
   n<- 50
   x1<- cbind( runif(n), runif(n))*100
   x2<-  cbind( runif(5), runif(5))
   #x2<- x1
   cardinalX<- cbind( runif(3), runif(3))
   m<- 2
   
   # simple check of marginal variances 
   look<- Tps.cov( x1,x1,cardinalX, m=m)
   look2<- Tps.cov( x1,cardinalX=cardinalX, m=m, marginal=TRUE)
   test.for.zero(diag(look), look2, tag="Tps.cov marginal" )
   
   
## comparing with the Tps function   
   data("ozone2")
   s<- ozone2$lon.lat
   y<- ozone2$y[16,]
   
   good<- !is.na( y)
   s<- s[good,]
   y<- y[good]
##### Tps used as benchmark
   
   out0<- Tps( s,y, scale.type ="unscaled", method="REML")
   lambdaHat<- out0$lambda.est[6,1]
   fHat<- predict( out0)
   
   cardinalX<- s[1:3,]
   out2<- mKrig( s,y, cov.function="Tps.cov",
                 cov.args= list( cardinalX=cardinalX,
                                 aRange=NA),
                 m=2,lambda=lambdaHat
   )
   
   # should be invariant to cardinal points and not need aRange =NA
   # defaults supplied by spatialProcessSetDefaults when detecting Tps.cov
   out3<- spatialProcess( s, y, cov.function="Tps.cov",
                             mKrig.args = list(m=2,  NtrA = 200),
                            lambda=lambdaHat )
   
   out4<- spatialProcess( s, y, cov.function="Tps.cov",
                          cov.args= list( 
                            aRange=NA),
                          #REML=TRUE,
                          verbose=FALSE, gridN=50
   )
   
   test.for.zero(fHat, predict( out2, s), tag="Tps vs spatialProcess" )
   # other parameters and likelihood
   test.for.zero(out0$tauHat.MLE, 
                 out3$summary["tau"],tag="Tps vs spatialProcess tau" )
   test.for.zero(fHat, predict( out3, s),
                 tag="Tps vs spatialProcess default" )
   # eff.df exact for spatialProcess because nTrA is larger than 
   # sample size.
   test.for.zero( out0$eff.df, out3$summary["eff.df"],
                  tag="Eff.df Tps vs spatialProcess")
   # look at prediction standard error computation.
   SE0<- predictSE( out0, s)
   SE3<- predictSE( out3,s)
   test.for.zero(SE0, SE3, tag="Tps vs spatialProcess SE")
# compare REML for Tps and spatialProcess  
   # test.for.zero(out0$lambda.est["REML", "-lnLike Prof"], 
   #               out4$summary["lnProfileREML.FULL"],
   #               tag="Tps vs spatialProcess  REML log like" )
   # test.for.zero(out0$lambda.est["REML", "lambda"], 
   #               out4$summary["lambda"],
   #               tag="Tps vs spatialProcess w REML lambda Hat" )
   # 
   gridList<- fields.x.to.grid( s, nx=20, ny=20)
   
   sGrid<- make.surface.grid( gridList)
   
   fHatGrid0<- predict( out0, sGrid)
   fHatGrid3<- predict( out3,sGrid)
   test.for.zero(fHatGrid0,fHatGrid3,
                 tag="Tps vs spatialProcess predictions on a grid")
   
   fHatGrid0SE<- predict( out0, sGrid)
   fHatGrid3SE<- predict( out3,sGrid)
   test.for.zero(fHatGrid0SE,fHatGrid3SE,
                 tag="Tps vs spatialProcess SE predictions on a grid")
   
options( echo=TRUE)
cat("all done testing Tps and spatialProcess with Tps.cov", fill=TRUE)
# 
# lambda0<- out0$lambda.est["REML", "lambda"]
# 
# Krig.flplike(out0, lambda0)