# LatticeKrig # Copyright 2004-2011, Institute for Mathematics Applied Geosciences # University Corporation for Atmospheric Research # Licensed under the GPL -- www.gpl.org/licenses/gpl.html suppressMessages(library( LatticeKrig)) options( echo=FALSE) # tests for computing the determinant and quad form test.for.zero.flag<- 1 alpha<- c(1,.5,.5) nlevel<-3 # in order to give different a.wght by level this should be in the # form of a list with each value (not just a vector) a.wght<-list(5,5,10) lnDet<- function(A){ sum( log( eigen( A, symmetric=TRUE)$values))} data( ozone2) x<-ozone2$lon.lat[1:20,] y<- ozone2$y[16,1:20] good <- !is.na( y) x<- x[good,] y<- y[good] #x<- transformx(x, "range") N<- length( y) lambda <- .8 set.seed(123) weights<- runif(N) W<- diag(weights) # a micro sized lattice so determinant is not too big or small # obj<- LKrig( x,y,NC=5, weights= weights, lambda=lambda, # nlevel=nlevel,alpha=alpha,a.wght=a.wght, # NtrA=5,iseed=122) LKinfo<- LKrigSetup( x,NC=5, lambda=lambda, nlevel=nlevel,alpha=alpha,a.wght=a.wght, fixedFunctionArgs = list(m = 2)) obj<- LKrig( x,y,LKinfo=LKinfo, weights= weights, NtrA=5,iseed=122) # The equivalent Kriging using mKrig and the LatticeKrig covariance function obj0<- mKrig( x,y, weights= weights, lambda=lambda, m=2, cov.function="LKrig.cov", cov.args=list(LKinfo=LKinfo), NtrA=5, iseed=122) # quick up front test of likelihood computation # This is only valid if mKrig is computating the likelihood correctly ;-). test.for.zero( obj$lnProfileLike.FULL, obj0$summary["lnProfileLike.FULL"], tag="lnProfileFULL from LKrig and identical compuation using mKrig") ###### check of formula with weights PHI<- as.matrix(LKrig.basis( x,LKinfo)) Q <- as.matrix(LKrig.precision(LKinfo)) M1<- PHI%*%solve( Q)%*%t(PHI) + lambda*solve( W) B1<- (t(PHI)%*%(W/lambda)%*%PHI + Q) B2<- (1/lambda) * ( t(PHI)%*%(W)%*%PHI + lambda*Q) B3<- t(PHI)%*%(W)%*%PHI + lambda*Q N2<- nrow(Q) # lnDet( M1) # lnDet( B1) - lnDet(Q) - lnDet( W/lambda) # lnDet( B2) - lnDet(Q) - lnDet( W/lambda) test.for.zero( lnDet( B3) - lnDet(Q) - sum( log( weights)) + (N-N2)*log(lambda), lnDet( M1), tag="Direct formula") test.for.zero( obj$lnDetCov, obj0$lnDetCov, tag= "lnDetCov for mKrig and LatticeKrig") test.for.zero( obj$sigma2.MLE, obj0$summary["sigma2"], tag= "quadratic forms for sigma2 MLE") test.for.zero( obj$lnProfileLike, obj0$summary["lnProfileLike.FULL"], tag="Profile Likelihood concentrated on lambda" ) # with replicated fields -- first pass including if LKrigSetup works and no fixed part. yRep<- cbind( y, y+5, y*1.2) LKinfo<- LKrigSetup( x,NC=5, weights= weights, lambda=lambda, nlevel=nlevel,alpha=alpha,a.wght=a.wght, fixedFunction=NULL) # fixedFunctionArgs = list(m = 1)) objR<- LKrig( x,yRep, LKinfo= LKinfo, weights= weights, NtrA=5,iseed=122) # The equivalent Kriging using mKrig and the LatticeKrig covariance function obj0R<- mKrig( x,yRep, weights= weights, lambda=LKinfo$lambda, m=0, cov.function="LKrig.cov", cov.args=list(LKinfo=LKinfo), collapseFixedEffect = FALSE, NtrA=5, iseed=122) # quick up front test of likelihood computation # This is only valid if mKrig is computating the likelihood correctly ;-). test.for.zero( objR$sigma2.MLE, obj0R$replicateInfo$sigma2.MLE, tag="individual sigma2 estimates from LKrig and identical compuation using mKrig") test.for.zero( objR$sigma2.MLE.FULL,obj0R$summary["sigma2"], tag=" full sigma2 estimate from LKrig and identical compuation using mKrig") test.for.zero( objR$lnProfileLike,obj0R$replicateInfo$lnProfileLike, tag=" individual lnProfile from LKrig and identical computation using mKrig") # try this out with replicated data and fixed part Y<- cbind( y,2*y, y*3) obj<- LKrig( x,Y,NC=5, weights= weights, lambda=lambda,nlevel=nlevel,alpha=alpha,a.wght=a.wght, NtrA=5,iseed=122) LKinfo<- obj$LKinfo # now check these formulas as implemented in LatticeKrig obj0<- mKrig( x,Y, weights= weights, lambda=lambda, m=2, cov.function="LKrig.cov", cov.args=list(LKinfo=LKinfo), NtrA=20, iseed=122, collapseFixedEffect = FALSE) test.for.zero( obj0$repInfo$sigma2.MLE, obj$sigma2.MLE, tag="MLEs for sigma2 with replicate fields" ) test.for.zero( obj0$repInfo$lnProfileLike, obj$lnProfileLike, tag="Profile Likelihood concentrated on lambda with replicate fields" ) # test of full likelihood # t(y- T%*%d.coef) %*% solve( tau^2/w + sigma2*K) %*% ( y- T%*%d.coef) # (1/sigma2) t(y- T%*%d.coef) %*% solve( lambda/w + K) %*% ( y- T%*%d.coef) # (1/sigma2) * quad.form # NOTE d.coef the MLE is GLS estimate and only depends on _lambda_ and not sigma2 # small data set with reps and weights # data( ozone2) x<-ozone2$lon.lat[1:20,] y<- ozone2$y[16,1:20] good <- !is.na( y) x<- x[good,] y<- y[good] N<- length( y) set.seed(123) weights<- runif(N) W<- diag(weights) tau<- .1 sigma2<- 2.2 lambda<- tau^2/sigma2 obj.test<- LKrig( x,Y,NC=5, weights= weights, nlevel=nlevel,alpha=alpha,a.wght=a.wght, NtrA=5,iseed=122, tau=tau, sigma2=sigma2) Tmatrix <- cbind(fields.mkpoly(x, 2)) res<- Y - Tmatrix%*%obj.test$d.coef # the M matrix you have all been waiting for! M<- (diag(lambda/weights) + LKrig.cov(x,x,LKinfo=obj.test$LKinfo)) qtest<- diag( t(res)%*% solve( M)%*% res ) test.for.zero( qtest, obj.test$quad.form, tag="quad form for arbitrary tau and sigma2") n<- nrow( Y) lnDetCov.test<- sum( log( eigen(M)$values)) test.for.zero( obj.test$lnDetCov, lnDetCov.test, tag="lnDetCov with arbitrary sigma2 tau") lnLike.test<- (-qtest/(2*sigma2) - log(2*pi)*(n/2) -(n/2)*log(sigma2) - (1/2) * lnDetCov.test) test.for.zero( lnLike.test, obj.test$lnLike, tag="lnLike arbitrary tau and sigma2") cat("all done with lnPLike with weights", fill=TRUE) options( echo=FALSE)