R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: i386-apple-darwin9.8.0/i386 (32-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
Package 'spam' is loaded. Spam version 0.23-0 (2010-09-01).
Type demo( spam) for some demos, help( Spam) for an overview
of this package.
Help for individual functions is optained by adding the
suffix '.spam' to the function name, e.g. 'help(chol.spam)'.

Attaching package: 'spam'

The following object(s) are masked from 'package:base':

    backsolve, forwardsolve, norm

 Use help(fields) for an overview of this library

library( fields, keep.source=TRUE) retains comments in the source code.

Copyright 2004-2011, Licensed under GPL, www.gpl.org/licenses/gpl.html 
> 
> #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)
>       
>