https://github.com/cran/fields
Tip revision: 04279c16ce718b82025aae97f6fcd30ba3a8b1c5 authored by Doug Nychka on 05 September 2011, 20:18:33 UTC
version 6.6
version 6.6
Tip revision: 04279c1
MLEfast.R
MLE.objective.fn<- function( ltheta,info, value=TRUE){
# marginal process covariance matrix
y<- as.matrix(info$y)
x<- info$x
smoothness<- info$smoothness
ngrid<-info$ngrid
# number of reps
M<- ncol(y)
Tmatrix <- fields.mkpoly(x,2)
qr.T <- qr(Tmatrix)
N<- nrow(y)
Q2 <- qr.yq2(qr.T, diag( 1,N))
ys<- t(Q2)%*%y
N2 <- length( ys)
theta<- exp( ltheta)
K<- Matern(rdist(x,x)/ theta ,
smoothness= smoothness)
Ke<- eigen( t(Q2)%*%K%*%Q2, symmetric=TRUE)
u2<- t(Ke$vectors)%*%ys
# mean over replicates -- mean square for coefficients for
# a particular eigenfunction.
u2.MS<- c(rowMeans(u2**2))
D2<- Ke$values
N2<- length( D2)
# grid of lambda based on spacing of eigenvalues
ngrid<- min( ngrid,N2)
lambda.grid<- exp(seq( log(D2[1]), log(D2[N2]),,ngrid))
trA<-minus.pflike<- rep( NA, ngrid)
#grid search followed by golden section
# -log likelihood
temp.fn<- function(llam,info){
lam.temp<- exp( llam)
u2<- info$u2.MS
D2<- info$D2
N2<- length(u2.MS)
# MLE of rho
rho.MLE<- (sum( (u2.MS)/(lam.temp + D2 )))/N2
# ln determinant
lnDetCov<- sum( log(lam.temp + D2) )
-1*M*(-N2/2 - log(2*pi)*(N2/2) - (N2/2)*log(rho.MLE) - (1/2) * lnDetCov ) }
# information list for calling golden section search.
info<- list( D2=D2, u2=u2.MS, M=M)
out<- golden.section.search(
f=temp.fn, f.extra=info, gridx= log(lambda.grid),
tol=1e-7)
minus.LogProfileLike<- out$fmin
lambda.MLE<- exp( out$x )
rho.MLE<- (sum( (u2.MS)/(lambda.MLE + D2)))/N2
sigma.MLE<- sqrt(lambda.MLE* rho.MLE)
trA<- sum( D2/(lambda.MLE +D2))
pars<- c( rho.MLE, theta, sigma.MLE,trA )
names(pars)<- c( "rho", "theta", "sigma", "trA")
if(value){
return(minus.LogProfileLike)}
else{
return( list( minus.lPlike=minus.LogProfileLike,
lambda.MLE= lambda.MLE,pars=pars,
mle.grid= out$coarse.search) )}
}
MLE.Matern.fast<- function(x,y,smoothness, theta.grid=NULL,
ngrid=20, verbose=FALSE,m=2,... ){
# 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] }
# list to pass to the objective function
# NOTE: large ngrid here is very cheap after the eigen decomposition
# has been done.
info<- list( x=x, y=y, smoothness= smoothness, ngrid=80)
# if grid for ranges is missing use some quantiles of
# pairwise distances among data.
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)}
else{
ngrid<- length( theta.grid)}
# grid search/golden section search
# note that search is in log scale.
out<- golden.section.search(
f=MLE.objective.fn, f.extra=info, gridx= log(theta.grid))
theta.MLE<- exp(out$x)
REML<- -out$fmin
# one final call with the theta.MLE value to recover MLEs for rho and sigma
out2<- MLE.objective.fn( log(theta.MLE), info,value=FALSE)
return( list(smoothness=smoothness,
pars= out2$pars[1:3], REML=REML,trA= out2$pars[4],
REML.grid= cbind( theta.grid, -1*out$coarse.search[,2]) ) )
}