https://github.com/cran/fields
Tip revision: 5a1a5a6622a7bae2a2bf0968c65f85a638fbf534 authored by Doug Nychka on 02 June 2011, 00:00:00 UTC
version 6.5.2
version 6.5.2
Tip revision: 5a1a5a6
REMLtest.R
##################################
######### testing functions
##################################
MaternGLS.test<- function(x,y, smoothness=1.5,init= log( c(1.0,.2,.1))){
# some simulations within fields to
# study the variability in estimates of the covariance parameters.
N<- length( y)
Tmatrix <- fields.mkpoly(x,2)
qr.T <- qr(Tmatrix)
Q2 <- qr.yq2(qr.T, diag( 1,N))
nu<- smoothness
loglmvn <- function(pars,nu,y,x,d){
lrho=pars[1]; ltheta=pars[2]; lsig2=pars[3]
# print( pars)
N <-length(y)
M <- ( exp(lrho)*Matern(d,range=exp(ltheta),smoothness=nu) +
exp(lsig2)*diag(N) )
X<- fields.mkpoly( x, 2)
Mi<- solve( M)
betahat<- solve(t(X)%*%Mi%*%X)%*% t(X)%*% Mi%*% y
res<- y - X%*%betahat
chol(M)-> cM
lLike<- (N/2)*log(2*pi) - (1/2)* (2*sum( log( diag(cM)))) - (1/2)* t(res)%*% Mi %*% res
ycept<- -lLike
if((abs(lrho)>20)| (abs(ltheta)>10)| (abs(lsig2)>40)) {
return( ycept + 1000*sum(abs(abs(pars)-100))) }
else{
return( ycept)}
}
d<- rdist( x,x)
temp<- optim(init,loglmvn,method="L-BFGS-B",nu=nu,y=y,x=x,d=d)
out<- exp(temp$par)
return( list( smoothness=smoothness, pars= out, optim=temp))
}
MaternGLSProfile.test<- function(x,y, smoothness=1.5,init= log( c( .05,1.0))){
# some simulations within fields to
# study the variability in estimates of the covariance parameters.
N<- length( y)
nu<- smoothness
loglmvn <- function(pars,nu,y,x,d){
llam=pars[1]; ltheta=pars[2]
# print( pars)
N <-length(y)
theta<- exp( ltheta)
lambda<- exp(llam)
lLike<- mKrig( x,y, theta=theta,Covariance="Matern", smoothness=nu, lambda=lambda)$lnProfileLike
ycept<- -lLike
if((abs(llam)>20)| (abs(ltheta)>10) ) {
return( ycept + 1000*sum(abs(abs(pars)-100))) }
else{
return( ycept)}
}
d<- rdist( x,x)
temp<- optim(init,loglmvn,method="L-BFGS-B",nu=nu,y=y,x=x,d=d)
out<- exp(temp$par)
rho.MLE<-mKrig( x,y, theta=out[2],Covariance="Matern", smoothness=nu, lambda=out[1])$rho.MLE
sigma2.MLE<- out[1]*rho.MLE
return( list( smoothness=smoothness,pars=c(rho.MLE, out[2], sigma2.MLE),
optim=temp))
}
MaternQR.test<- function(x,y, smoothness=1.5,init= log( c(1.0,.2,.1))){
# some simulations within fields to
# study the variability in estimates of the covariance parameters.
nu<- smoothness
loglmvn <- function(pars,nu,x,y){
N<- length( y)
Tmatrix <- fields.mkpoly(x,2)
qr.T <- qr(Tmatrix)
Q2 <- qr.yq2(qr.T, diag( 1,N))
ys<- t(Q2)%*%y
N2<- length( ys)
lrho=pars[1]; ltheta=pars[2]; lsig2=pars[3]
d<- rdist( x,x)
A <- ( Matern(d,scale=exp(lrho),range=exp(ltheta),smoothness=nu) + exp(lsig2)*diag(N) )
A <- t(Q2)%*%A%*%Q2
A<- chol( A)
w = backsolve(A,ys,transpose=TRUE)
ycept<- (N2/2)*log(2*pi) + sum(log(diag(A))) + (1/2)*t(w)%*%w
if((abs(lrho)>100)| (abs(ltheta)>100)| (abs(ltheta)>100)) {
return( ycept + 1000*sum(abs(abs(pars)-100))) }
else{
return( ycept)}
}
temp<- optim(init,loglmvn,method="L-BFGS-B",nu=nu,x=x,y=y)
out<- exp(temp$par)
llike<- loglmvn( temp$par, nu,x,y)
return( list( smoothness=smoothness,pars= out, llike=llike, optim=temp))
}
MaternQRProfile.test<- function(x,y, smoothness=1.5,init= log( c(1.0)) ){
# some simulations within fields to
# study the variability in estimates of the covariance parameters.
nu<- smoothness
loglmvn <- function(pars,nu,x,y){
ltheta=pars[1]
# print( exp(ltheta))
ycept<- Krig( x,y, Covariance="Matern", theta= exp( ltheta),
smoothness=nu, method="REML")$lambda.est[6,5]
# print( c(exp(ltheta),ycept))
if((abs(ltheta)>100)){
return( ycept + 1000*sum(abs(abs(pars)-100))) }
else{
return( ycept)}
}
temp<- optim(init,loglmvn,method="L-BFGS-B",nu=nu,x=x,y=y)
theta.est<- exp(temp$par[1])
out2<-Krig( x,y, Covariance="Matern", theta= theta.est,smoothness=nu, method="REML")
# MLE based on reduced degrees of freedom:
offset<- (out2$N/( out2$N-3))
out3<- c( out2$rho.MLE*offset, theta.est, out2$shat.MLE^2*offset)
return( list(obj=out2, smoothness=smoothness,pars= out3, trA= out2$eff.df, optim=temp))
}
# this function has correct formula for REML likelihood
REML.test<-function (x,y,rho,sigma2, theta, nu=1.5){
Tmatrix <- fields.mkpoly(x,2)
qr.T <- qr(Tmatrix)
N<- length( y)
Q2 <- qr.yq2(qr.T, diag( 1,N))
ys<- t(Q2)%*%y
N2 <- length( ys)
A <- (rho* Matern(rdist(x,x),range=theta,smoothness=nu) + sigma2*diag(1,N) )
A <- t(Q2)%*%A%*%Q2
Ac<- chol( A)
w <- backsolve(Ac,ys,transpose=TRUE)
REML.like<- (N2/2)*log(2*pi) + (1/2)*2*sum(log(diag(Ac))) + (1/2)*t(w)%*%w
REML.like<- -1* REML.like
ccoef<- rho* Q2%*%solve( A)%*%ys
return( list(REML.like= REML.like, A=A, ccoef=ccoef,
quad.form=t(w)%*%w, rhohat= (t(w)%*%w/N2)*rho, det=2*sum(log(diag(Ac))), N2=N2 ) )
}
MLE.Matern<- function(x,y,smoothness, theta.grid=NULL,
ngrid=20, verbose=FALSE,m=2,niter=25, tol=1e-5,... ){
# remove missing values and print out a warning
bad<- is.na( y)
if( sum(bad)>0){
cat("removed ",sum(bad)," NAs", fill=TRUE)
x<- x[!bad,]
y<- y[!bad] }
objective.fn<-
function( ltheta, info){
minus.lPLike<- Krig( info$x,info$y,Covariance="Matern",
smoothness= info$smoothness, theta=exp(ltheta),
method="REML", nstep.cv=80, give.warnings=FALSE,m=m,...)$lambda.est[6,5]
return(minus.lPLike)}
# list to pass to the objective function
info<- list( x=x, y=y, smoothness= smoothness)
#
# if grid for ranges is missing use some quantiles of pairwise distances among data.
# this will only work if the likelihood at endpoints is smaller than middle.
# (i.e. convex)
if( is.null(theta.grid)){
theta.range<- quantile( rdist( x,x), c( .03,.97))
theta.grid <- seq( theta.range[1], theta.range[2],,ngrid)}
if( length( theta.grid)==2){
theta.grid <- seq( theta.grid[1], theta.grid[2],,ngrid)}
ngrid<- length( theta.grid)
sighat<- rhohat<- trA<-theta<- rep( NA,ngrid)
minus.REML<- rep( NA, ngrid)
# grid search
for( j in 1:ngrid){
minus.REML[j] <- objective.fn( log(theta.grid[j]), info)
}
IMIN<- ( 1: ngrid)[ min(minus.REML)==minus.REML]
# best point for theta from grid search
if( IMIN ==1 |IMIN ==ngrid){
cat("REML at end of search interval:", fill=TRUE)
temp<- cbind(theta.grid,-minus.REML)
dimnames(temp)<- list( NULL, c("theta", "logProfileLike"))
return( list(smoothness=smoothness, pars= rep( NA,3) ,
REML=NA, trA=NA,REML.grid= temp))}
lstart<- log( theta.grid)[ IMIN + c( -1,0,1)]
# golden section search -- this assumes convex minus log likelihood
# note that search is in log scale.
out<- golden.section.search( lstart[1], lstart[2], lstart[3], f=objective.fn, f.extra=info,
niter=niter, tol=tol)$x
theta.MLE<- exp(out)
# one final call to Krig with the theta.MLE value to recover MLEs for rho and sigma
hold<- Krig( info$x,info$y,Covariance="Matern",
smoothness= info$smoothness, theta=theta.MLE,
method="REML")
sig.MLE<- hold$shat.MLE
rhohat<- hold$rhohat
trA<- hold$lambda.est[6,2]
REML<- hold$lambda.est[6,5] # actually minus log REML
out<- c( rhohat, theta.MLE,sig.MLE)
names(out)<- c( "rho","theta", "sigma")
temp<- cbind(theta.grid,-minus.REML)
dimnames(temp)<- list( NULL, c("theta", "logProfileLike"))
return( list(smoothness=smoothness, pars= out,REML=-REML, trA=trA,REML.grid= temp) )
}