# # 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) #options( echo=TRUE) # ##### generate test data # data( ozone2) # x is a two column matrix where each row is a location in lon/lat # coordinates x<- ozone2$lon.lat # y is a vector of ozone measurements at day 16 a the locations. y<- ozone2$y[16,] #ind<- !is.na( y) #x<- x[ind,] #y<- y[ind] x<-x[1:31,] y<-y[1:31] y[31]<-NA ################ test that optim results match the model evaluated ################ at the optimized parameters. optim.args = list(method = "BFGS", control = list(fnscale = -1, parscale = c(0.5, 0.5), ndeps = c(0.05,0.05))) MLEfit0 <- mKrigMLEJoint(x, y, cov.params.start= list(lambda=.5, aRange=1.2), cov.fun="stationary.cov", optim.args=optim.args, cov.args = list(Covariance = "Matern", smoothness=1.0), na.rm=TRUE, mKrig.args = list( m=1), verbose=FALSE) # check agreement with a fast return note cov.params.start and lambda.fixed # are the switches to indicate this case MLEfit0C <- mKrigMLEJoint(x, y, cov.params.start = NULL, cov.function = "stationary.cov", cov.args = list(Covariance = "Matern", lambda = MLEfit0$pars.MLE["lambda"], smoothness = 1.0, aRange = MLEfit0$pars.MLE["aRange"]), na.rm = TRUE, mKrig.args = list( m=1), verbose = FALSE ) test.for.zero( MLEfit0$summary["lnProfileLike.FULL"], MLEfit0C$summary["lnProfileLike.FULL"], tag="Likelihood Values optim and the fast return") test.for.zero( MLEfit0$summary["lnProfileLike.FULL"], MLEfit0$optimResults$value, tag="Likelihood Values summary and optim") obj0<- mKrig( x,y, cov.args = list(Covariance = "Matern", smoothness = 1.0), na.rm=TRUE, m=1, lambda= MLEfit0$pars.MLE["lambda"], aRange=MLEfit0$pars.MLE["aRange"]) test.for.zero( MLEfit0$summary["lnProfileLike.FULL"], obj0$summary["lnProfileLike.FULL"], tag="Likelihood Values summary and direct mKrig call") test.for.zero( MLEfit0$summary["sigma2"],obj0$summary["sigma2"], tag="... and sigma^2.MLE") # test that grid seraching is correct aRange.MLE<- MLEfit0$summary["aRange"] par.grid<- list( aRange= c(.5, 1.0, 1.5)*aRange.MLE ) MLEfit1<- mKrigMLEGrid(x, y, cov.fun = "stationary.cov", cov.args = list(Covariance = "Matern", smoothness = 1.0 ), par.grid = par.grid, mKrig.args = list( m=1), na.rm = TRUE, verbose = FALSE, cov.params.start = list( lambda = .2) ) hold<- (MLEfit1$summary[1,"lnProfileLike.FULL"] < MLEfit1$summary[2,"lnProfileLike.FULL"]) & (MLEfit1$summary[3,"lnProfileLike.FULL"] < MLEfit1$summary[2,"lnProfileLike.FULL"]) test.for.zero(as.numeric(hold), 1, relative=FALSE, tag="consistency of Likelihood values") ########################## ### now evaluate on the "grid" of lambdas found by profiling lambda.MLEs<- MLEfit1$summary[,"lambda"] par.grid<- list( lambda = lambda.MLEs, aRange = c(.5, 1.0, 1.5)*aRange.MLE ) MLEfit1B<- mKrigMLEGrid(x, y, cov.function = "stationary.cov", cov.args = list(Covariance = "Matern", smoothness = 1.0), par.grid = par.grid, mKrig.args = list( m=1), na.rm = TRUE, verbose = FALSE) tempCol<- c( "lnProfileLike.FULL", "lambda", "tau","sigma2") test.for.zero( as.matrix(MLEfit1$summary[,tempCol]), as.matrix(MLEfit1B$summary[,tempCol]), tag="grid search with and w/o profile") par.grid<- list( lambda = c(.999, 1.0, 1.001)*MLEfit0$summary["lambda"], aRange = rep(MLEfit0$summary["aRange"] ,3 ) ) MLEfit2 <- mKrigMLEGrid(x, y, cov.function = "stationary.cov", cov.args = list(Covariance = "Matern", smoothness = 1.0), mKrig.args = list( m=1), par.grid = par.grid, verbose = FALSE) hold<- (MLEfit2$summary[1,"lnProfileLike.FULL"] < MLEfit2$summary[2,"lnProfileLike.FULL"]) & (MLEfit2$summary[3,"lnProfileLike.FULL"] < MLEfit2$summary[2,"lnProfileLike.FULL"]) test.for.zero(as.numeric(hold), 1, relative=FALSE, tag="crude test of maxmimum") #MLEfit3<- MLESpatialProcess( x,y, # cov.args = list(Covariance = "Matern", # smoothness = 1.0), # mKrig.args = list( m=1), # cov.params.start = list( lambda =.2, aRange = NA) # ) MLEfit3<- spatialProcess( x,y, cov.args = list(Covariance = "Matern", smoothness = 1.0), mKrig.args = list( m=1), cov.params.start = list( lambda =.2) ) test.for.zero(MLEfit0$summary[1:2]/ (MLEfit3$summary[1:2]), 1, tol=1e-5, tag="Testing MLE from spatialProcess ") ######### making sure spatialProcess uses parameter information correctly obj<- spatialProcess( x, y, mKrig.args= list(m = 1), lambda= MLEfit0$summary["lambda"], aRange = MLEfit0$summary["aRange"] ) obj1<- spatialProcess( x, y, mKrig.args= list(m = 1), ) test.for.zero(MLEfit0$summary[1], obj$summary["lnProfileLike.FULL"], tag="spatialProcess finding MLE " ) test.for.zero(MLEfit0$summary[1], obj1$summary["lnProfileLike.FULL"], tol=1e-5, tag="spatialProcess given MLE " ) # testing Krig function #out1<- Krig( x,y, cov.fun="stationary.cov", # cov.args = list(Covariance = "Matern", # smoothness=1.0, aRange=.9), # na.rm=TRUE, # m=2) #generate observation locations set.seed( 22) n=100 x = matrix(runif(2*n), nrow=n) #generate observations at the locations trueARange = .1 trueLambda = .1 distanceMatrix<- rdist(x,x) Sigma<- Matern( distanceMatrix/trueARange, smoothness=1.0 ) U = chol(Sigma) M<- 2e3 # lots of replicated fields. set.seed( 332) y = t(U)%*%matrix( rnorm(n*M), n,M) + sqrt(trueLambda)*matrix( rnorm(n*M), n,M) out<- mKrig( x,y, lambda=trueLambda, aRange=trueARange, cov.function ="stationary.cov",cov.args = list(Covariance = "Matern", smoothness=1.0) ) optim.args = list(method = "BFGS", control = list(fnscale = -1, ndeps = c(0.09,0.09))) MLEfitA <- mKrigMLEJoint(x, y, cov.params.start= list(aRange=.2, lambda=.01), cov.function="stationary.cov", optim.args=optim.args, cov.args = list(Covariance = "Matern", smoothness=1.0), na.rm=TRUE, reltol = 1e-6, mKrig.args = list( m=0), verbose=FALSE) cat("Testing mKrigMLEJoint against true values", fill=TRUE) test.for.zero( MLEfitA$summary["lambda"],.1, tol=.006) test.for.zero( MLEfitA$summary["aRange"],.1, tol=.02) test.for.zero( MLEfitA$summary["sigma2"], 1.0, tol=.02) ### now test REML fitting MLEfitB <- mKrigMLEJoint(x, y, cov.params.start= list(aRange=.12, lambda=.5), cov.function="stationary.cov", optim.args=optim.args, cov.args = list(Covariance = "Matern", smoothness=1.0), na.rm=TRUE, mKrig.args = list( m=0), REML=TRUE, verbose=FALSE) cat("Testing mKrigMLEJoint with REML against true values", fill=TRUE) test.for.zero( MLEfitB$summary["lambda"],.1, tol=.007) test.for.zero( MLEfitB$summary["aRange"],.1, tol=.01) test.for.zero( MLEfitB$summary["sigma2"], 1.0, tol=.01) cat("Testing mKrigMLEJoint with REML FALSE against true values", fill=TRUE) MLEfitC <- mKrigMLEJoint(x, y, cov.params.start= list(aRange=.12, lambda=.5), cov.function ="stationary.cov", optim.args=optim.args, cov.args = list(Covariance = "Matern", smoothness=1.0), na.rm=TRUE, mKrig.args = list( m=2), REML=FALSE, verbose=FALSE ) test.for.zero( MLEfitC$summary["lambda"], .1, tol=.02) test.for.zero( MLEfitC$summary[ "aRange"], .1, tol=.02) test.for.zero( MLEfitC$summary["sigma2"], 1.0, tol=.01) cat("all done with mKrigMLEGrid tests", fill=TRUE) options( echo=TRUE)