# # 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) y<- ozone2$y[16,] x<- ozone2$lon.lat # # Omit the NAs good<- !is.na( y) x<- x[good,] y<- y[good] #source("~/Home/Src/fields/R/mKrig.family.R") # now look at mKrig w/o sparse matrix look<- mKrig( x,y, cov.function="stationary.cov", aRange=10, lambda=.3, chol.args=list( pivot=FALSE)) lookKrig<- Krig( x,y, cov.function="stationary.cov", aRange=10) test.df<-Krig.ftrace(look$lambda,lookKrig$matrices$D) test<- Krig.coef( lookKrig, lambda=look$lambda) test.for.zero( look$d, test$d, tag="Krig mKrig d coef") test.for.zero( look$c, test$c, tag="Krig mKrig c coef") # test of trace calculation look<- mKrig( x,y, cov.function="stationary.cov", aRange=10, lambda=.3, find.trA=TRUE, NtrA= 1000, iseed=243) test.for.zero( look$eff.df, test.df,tol=.01, tag="Monte Carlo eff.df") # lookKrig<-Krig( x,y, cov.function="stationary.cov", aRange=350, Distance="rdist.earth",Covariance="Wendland", cov.args=list( k=2, dimension=2) ) look<- mKrig( x,y, cov.function="stationary.cov", aRange=350, Distance="rdist.earth",Covariance="Wendland", cov.args=list( k=2, dimension=2), lambda=lookKrig$lambda, find.trA=TRUE, NtrA= 1000, iseed=243) test.for.zero( look$c, lookKrig$c, tag="Test of wendland and great circle") test.for.zero(look$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D) ,tol=.01, tag="eff.df") # same calculation using sparse matrices. look4<- mKrig( x,y, cov.function="wendland.cov", aRange=350, Dist.args=list( method="greatcircle"), cov.args=list( k=2), lambda=lookKrig$lambda, find.trA=TRUE, NtrA=500, iseed=243) test.for.zero( look$c.coef, look4$c.coef,tol=8e-7, tag="Test of sparse wendland and great circle") test.for.zero(look4$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D), tol=.01, tag="sparse eff.df") # great circle distance switch has been a big bug -- test some options look<- mKrig( x,y, cov.function="wendland.cov", aRange=350, Dist.args=list( method="greatcircle"), cov.args=list( k=2),lambda=lookKrig$lambda, find.trA=TRUE, NtrA=1000, iseed=243) test.for.zero(look$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D), tol=1e-2, tag="exact sparse eff.df") # compare to fast Tps look3<- fastTps( x,y,aRange=350,lambda=lookKrig$lambda, NtrA=200, iseed=243, lon.lat=TRUE) #look3$c<- lookKrig$c #look3$d<- lookKrig$d object<- look3 np<- object$np Ey <- diag(1, np) NtrA <- np hold <- predict.mKrig(object, ynew = Ey, collapseFixedEffect=FALSE) hold2<- matrix( NA, np,np) for( k in 1:np){ hold2[,k] <- predict.Krig(lookKrig, y = Ey[,k]) } #plot( diag(hold), diag(hold2)) test.for.zero( look3$c, lookKrig$c, tol=5e-7) test.for.zero( look3$d, lookKrig$d, tol=2e-8) test.for.zero( look3$fitted.values, lookKrig$fitted.values, tol=1e-7) test.for.zero( predict( look3, xnew= look3$x), predict( lookKrig, xnew= lookKrig$x), tol=5e-7) test.for.zero( hold[,1], hold2[,1], tol=1e-7, relative=FALSE) test.for.zero(diag(hold),diag(hold2), tol=2E-7, relative=FALSE, tag="exact sparse eff.df by predict -- fastTps") #plot( diag(hold), ( 1- diag(hold2)/ diag(hold)) ) test.for.zero(look3$eff.df,sum( diag(hold)) , tag="fastTps ef.df exact" ) test.for.zero(look3$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D), tol=2e-7, tag="exact sparse eff.df through mKrig-- fastTps") # calculations of likelihood, sigma and tau lam<-.2 out<- mKrig( x,y, cov.function =Exp.cov, aRange=4, lambda=lam) out2<- Krig( x,y, cov.function =Exp.cov, aRange=4, lambda=lam) Sigma<- Exp.cov( x,x,aRange=4) X<- cbind( rep(1, nrow(x)), x) Sinv<- solve( Sigma + lam* diag( 1, nrow( x))) #checks on likelihoods # quadratic form: betaHat<- c(solve( t(X)%*%Sinv%*%(X) ) %*% t(X) %*%Sinv%*%y) test.for.zero( betaHat, out$beta, tag="initial check on d for likelihood") r<- y -X%*%betaHat N<- nrow(x) look<- t( r)%*%(Sinv)%*%r/N test.for.zero( look, out$summary["sigma2"], tag="sigma2 hat from likelihood") test.for.zero( look, out2$sigma.MLE, tag="sigma2 hat from likelihood compared to Krig") # check determinant lam<- .2 Sigma<- Exp.cov( x,x,aRange=4) M<- Sigma + lam * diag( 1, nrow(x)) chol( M)-> Mc look2<- sum( log(diag( Mc)))*2 out<-mKrig( x,y,cov.function =Exp.cov, aRange=4, lambda=lam) test.for.zero( out$lnDetCov, look2) test.for.zero( out$lnDetCov, determinant(M, log=TRUE)$modulus) # weighted version lam<- .2 Sigma<- Exp.cov( x,x,aRange=4) set.seed( 123) weights<- runif(nrow( x)) M<- Sigma + diag(lam/ weights) chol( M)-> Mc look2<- sum( log(diag( Mc)))*2 out<-mKrig( x,y,weights=weights, cov.function =Exp.cov, aRange=4, lambda=lam) test.for.zero( out$lnDetCov, look2) test.for.zero( look2, determinant(M, log=TRUE)$modulus) test.for.zero( out$lnDetCov, determinant(M, log=TRUE)$modulus) # check profile likelihood by estimating MLE lam.true<- .2 N<- nrow( x) Sigma<- Exp.cov( x,x,aRange=4) M<- Sigma + lam.true * diag( 1, nrow(x)) chol( M)-> Mc t(Mc)%*%Mc -> test ##D set.seed( 234) ##D NSIM<- 100 ##D hold2<-rep( NA, NSIM) ##D temp.fun<- function(lglam){ ##D out<-mKrig( x,ytemp, ##D cov.function =Exp.cov, aRange=4, lambda=exp(lglam)) ##D return(-1* out$lnProfileLike)} ##D hold1<-rep( NA, NSIM) ##D yt<- rep( 1, N) ##D obj<- Krig( x,yt, aRange=4) ##D E<- matrix( rnorm( NSIM*N), ncol=NSIM) ##D for ( j in 1:NSIM){ ##D cat( j, " ") ##D ytemp <- x%*%c(1,2) + t(Mc)%*%E[,j] ##D out<- optim( log(.2), temp.fun, method="BFGS") ##D hold2[j]<- exp(out$par) ##D hold1[j]<- gcv.Krig(obj, y=ytemp)$lambda.est[6,1] ##D } ##D test.for.zero( median( hold1), .2, tol=.08) ##D test.for.zero( median( hold2), .2, tol=.12)