R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # this is a test script to verify the likelihood computations are 
> # correct with the eigen decomposition format used in Krig
> # see Krig.flplike for the concise computation.
> #  
> 
> library(fields)
Loading required package: spam
Loading required package: grid
Spam version 1.4-0 (2016-08-29) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: 'spam'

The following objects are masked from 'package:base':

    backsolve, forwardsolve

Loading required package: maps
> 
> #options( echo=FALSE)
> test.for.zero.flag<- 1
> 
> 
> data( ozone2)
> x<- ozone2$lon.lat
> y<- ozone2$y[16,]
> is.good <- !is.na( y)
> x<- x[is.good,]
> y<- y[is.good]
> 
> theta<- 2.0
> 
> # check log likelihood calculation
>   nu<- 1.5
>   lambda<- .2
>   out<- mKrig( x,y, theta=theta,Covariance="Matern", smoothness=nu, lambda=lambda) 
> 
> # peg rho and sigma as MLEs from mKrig
>   rho <- out$rho.MLE
>   sigma2<- rho*lambda
>   N<- length( y)
>   dd<- rdist( x,x)
>   M<-  rho* Matern( dd, range= theta, smoothness=nu) + sigma2* diag( 1, N)
>   X<- fields.mkpoly( x, 2)
>   Mi<- solve( M)
>   betahat<-  solve(t(X)%*%Mi%*%X)%*% t(X)%*% Mi%*% y
>   res<- y - X%*%betahat
>   ccoef<- ( Mi%*% ( res))*rho
> 
> # sanity check that estimates are the same
>   test.for.zero( ccoef, out$c, tag="check ccoef")
Testing:  check ccoef
PASSED test at tolerance  1e-08
> 
> # find full log likelihood
>   chol(M)-> cM
>   lLike<-  -(N/2)*log(2*pi) - (1/2)* (2*sum( log( diag(cM)))) - (1/2)* t(res)%*% Mi %*% res
> # formula for full likelihood using peices from mKrig
>   lLike.test<- -(N/2)*log(2*pi) - (1/2)* out$lnDetCov - (1/2)*(N)*log( rho) - (1/2)*out$quad.form/rho
> 
>   test.for.zero( lLike, lLike.test, tag="llike full verses rhohat")
Testing:  llike full verses rhohat
PASSED test at tolerance  1e-08
>   test.for.zero( lLike, out$lnProfileLike, tag="llike profile from mKrig")
Testing:  llike profile from mKrig
PASSED test at tolerance  1e-08
> 
> # REML check
>   nu<- 1.5
>   theta<- .6
>   obj<- Krig( x,y, theta=theta,Covariance="Matern", smoothness=nu )
> 
> # sanity check that c coefficients agree with Krig
>   rho<- 500
>   lambda<- .2
>   sigma2<- lambda*rho
> 
>   hold<- REML.test( x,y,rho, sigma2, theta, nu=1.5)
>   ccoef2<- Krig.coef( obj, lambda)$c
>   test.for.zero( hold$ccoef, ccoef2, tag="ccoefs")
Testing:  ccoefs
PASSED test at tolerance  1e-08
> 
> # check RSS with Krig decomposition.
>   RSS1<- sum( (lambda*ccoef2)**2)
>   lD <- obj$matrices$D * lambda
>   RSS2 <- sum(((obj$matrices$u * lD)/(1 + lD))^2) 
>   test.for.zero( RSS2, RSS1, tag=" RSS using matrices")
Testing:   RSS using matrices
PASSED test at tolerance  1e-08
> 
> # check quadratic form with Krig
>   D.temp<- obj$matrices$D[  obj$matrices$D>0]
>   A3test<- (1/lambda)* obj$matrices$V %*% diag((D.temp*lambda)/ (1 +D.temp*lambda) )%*% t( obj$matrices$V)
>   test.for.zero(solve(A3test), hold$A/rho, tol=5e-8)
PASSED test at tolerance  5e-08
>   Quad3<-   sum( D.temp*(obj$matrices$u[obj$matrices$D>0])^2/(1+lambda*D.temp))
> 
>   test.for.zero( hold$quad.form, Quad3/rho, tag="quad form")
Testing:  quad form
PASSED test at tolerance  1e-08
> 
> # test determinants
>   N2<- length( D.temp)
>   det4<- -sum( log(D.temp/(1 + D.temp*lambda)) )
>   det1<- sum( log(eigen( hold$A/rho)$values)) 
>   test.for.zero( det1, det4, tag="det" )
Testing:  det
PASSED test at tolerance  1e-08
> 
> # test REML Likelihood
>   lLikeREML.test<--1*( (N2/2)*log(2*pi) - (1/2)*(sum( log(D.temp/(1 + D.temp*lambda)) ) - N2*log(rho)) +
+                                       (1/2)*sum( lD*(obj$matrices$u)^2/(1+lD)) /(lambda*rho) )
> 
> test.for.zero( hold$REML.like, lLikeREML.test, tag="REML using matrices")
Testing:  REML using matrices
PASSED test at tolerance  1e-08
> 
> 
> # profile likelihood
> 
> # lnProfileLike <- (-np/2 - log(2*pi)*(np/2)
> #                      - (np/2)*log(rho.MLE) - (1/2) * lnDetCov)
> #  test using full REML likelihood.
>   nu<- 1.5
>   rho<- 7000
>   lambda<- .02
>   sigma2<- lambda*rho
>   theta<- 2.0
>   obj<- Krig( x,y, theta=theta,Covariance="Matern", smoothness=nu )
>   hold<- REML.test(x,y,rho, sigma2, theta, nu=1.5)
>   np<- hold$N2
>   rho.MLE<- c(hold$rhohat)
>   lnDetCov<-sum( log(eigen( hold$A/rho)$values))  
> 
>   l0<- REML.test(x,y,rho.MLE, rho.MLE*lambda, theta, nu=1.5)$REML.like
>   l1<-   (-np/2 - log(2*pi)*(np/2)- (np/2)*log(rho.MLE) - (1/2) * lnDetCov)
>   l2<-  (-1)*Krig.flplike( lambda, obj)
> 
>   test.for.zero( l0,l2, tag="REML profile flplike")
Testing:  REML profile flplike
PASSED test at tolerance  1e-08
>   test.for.zero( l1,l2, tag="REML profile flplike")
Testing:  REML profile flplike
PASSED test at tolerance  1e-08
> 
> 
> cat("all done with likelihood  tests", fill=TRUE)
all done with likelihood  tests
> options( echo=TRUE)
>       
> 
> proc.time()
   user  system elapsed 
  1.010   0.066   1.087