if (FALSE) {
#perform joint likelihood maximization over lambda and aRange.
# NOTE: optim can get a bad answer with poor initial starts.
data(ozone2)
s<- ozone2$lon.lat
z<- ozone2$y[16,]
gridList<- list( aRange = seq( .4,1.0,length.out=20),
lambda = 10**seq( -1.5,0,length.out=20)
)
par.grid<- make.surface.grid( gridList)
out<- mKrigMLEGrid( s,z, par.grid=par.grid,
cov.args= list(smoothness=1.0,
Covariance="Matern" )
)
outP<- as.surface( par.grid, out$summary[,"lnProfileLike.FULL"])
image.plot( outP$x, log10(outP$y),outP$z,
xlab="aRange", ylab="log10 lambda")
}
if (FALSE) {
N<- 50
set.seed(123)
x<- matrix(runif(2*N), N,2)
aRange<- .2
Sigma<- Matern( rdist(x,x)/aRange , smoothness=1.0)
Sigma.5<- chol( Sigma)
tau<- .1
# 250 independent spatial data sets but a common covariance function
# -- there is little overhead in
# MLE across independent realizations and a good test of code validity.
M<-250
F.true<- t( Sigma.5) %*% matrix( rnorm(N*M), N,M)
Y<- F.true + tau* matrix( rnorm(N*M), N,M)
# find MLE for lambda with grid of ranges
# and smoothness fixed in Matern
par.grid<- list( aRange= seq( .1,.35,,8))
obj1b<- mKrigMLEGrid( x,Y,
cov.args = list(Covariance="Matern", smoothness=1.0),
cov.params.start=list( lambda = .5),
par.grid = par.grid
)
obj1b$summary # take a look
# profile over aRange
plot( par.grid$aRange, obj1b$summary[,"lnProfileLike.FULL"],
type="b", log="x")
}
if (FALSE) {
# m=0 is a simple switch to indicate _no_ fixed spatial drift
# (the default and highly recommended is linear drift, m=2).
# However, m=0 results in MLEs that are less biased, being the correct model
# -- in fact it nails it !
obj1a<- mKrigMLEJoint(x,Y,
cov.args=list(Covariance="Matern", smoothness=1.0),
cov.params.start=list(aRange =.5, lambda = .5),
mKrig.args= list( m=0))
test.for.zero( obj1a$summary["tau"], tau, tol=.007)
test.for.zero( obj1a$summary["aRange"], aRange, tol=.015)
}
##########################################################################
# A bootstrap example
# Here is a example of a more efficient (but less robust) bootstrap using
# mKrigMLEJoint and tuned starting values
##########################################################################
if (FALSE) {
data( ozone2)
obj<- spatialProcess( ozone2$lon.lat,ozone2$y[16,] )
######### boot strap
set.seed(123)
M<- 25
# create M indepedent copies of the observation vector
ySynthetic<- simSpatialData( obj, M)
bootSummary<- NULL
aRangeMLE<- obj$summary["aRange"]
lambdaMLE<- obj$summary["lambda"]
for( k in 1:M){
cat( k, " " )
# here the MLEs are found using the easy top level level wrapper
# see mKrigMLEJoint for a more efficient strategy
out <- mKrigMLEJoint(obj$x, ySynthetic[,k],
weights = obj$weights,
mKrig.args = obj$mKrig.args,
cov.function = obj$cov.function.name,
cov.args = obj$cov.args,
cov.params.start = list( aRange = aRangeMLE,
lambda = lambdaMLE)
)
newSummary<- out$summary
bootSummary<- rbind( bootSummary, newSummary)
}
cat( " ", fill=TRUE )
obj$summary
stats( bootSummary)
}
if (FALSE) {
#perform joint likelihood maximization over lambda, aRange, and smoothness.
#note: finding smoothness is not a robust optimiztion
# can get a bad answer with poor initial guesses.
obj2<- mKrigMLEJoint(x,Y,
cov.args=list(Covariance="Matern"),
cov.params.start=list( aRange = .18,
smoothness = 1.1,
lambda = .08),
)
#look at lnLikelihood evaluations
obj2$summary
#compare to REML
obj3<- mKrigMLEJoint(x,Y,
cov.args=list(Covariance="Matern"),
cov.params.start=list(aRange = .18,
smoothness = 1.1,
lambda = .08),
, REML=TRUE)
obj3$summary
}
if (FALSE) {
#look at lnLikelihood evaluations
# check convergence of MLE to true fit with no fixed part
#
obj4<- mKrigMLEJoint(x,Y,
mKrig.args= list( m=0),
cov.args=list(Covariance="Matern", smoothness=1),
cov.params.start=list(aRange=.2, lambda=.1),
REML=TRUE)
#look at lnLikelihood evaluations
obj4$summary
# nails it!
}
Run the code above in your browser using DataLab