https://github.com/cran/fields
Tip revision: 6c8b30169bba182a68765ee3cb9b4e2ef7d38332 authored by Doug Nychka on 16 November 2011, 00:00:00 UTC
version 6.6.3
version 6.6.3
Tip revision: 6c8b301
Likelihood.test.Rout.save
R version 2.13.2 (2011-09-30)
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.27-0 (2011-08-17).
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
>
> #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)
>
>