# NOT RUN {
# }
# NOT RUN {
#perform joint likelihood maximization over lambda and theta.
# NOTE: optim can get a bad answer with poor initial starts.
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
obj<- mKrigMLEJoint(x,y,
cov.args=list(Covariance="Matern", smoothness=1.0),
cov.params.start=list(theta=.2, lambda = .1)
)
#
# check lnLikeihood evaluations that were culled from optim
# these are in obj$lnLike.eval
# funny ranges are set to avoid very low likelihood values
quilt.plot( log10(obj$lnLike.eval[,c("theta","lambda") ] ), obj$lnLike.eval[,"lnProfileLike.FULL"])
points( log10(obj$pars.MLE[1]), log10(obj$pars.MLE[2]), pch=16, col="grey" )
# some synthetic data with replicates
N<- 50
set.seed(123)
x<- matrix(runif(2*N), N,2)
theta<- .2
Sigma<- Matern( rdist(x,x)/theta , smoothness=1.0)
Sigma.5<- chol( Sigma)
sigma<- .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) -->
F.true<- t( Sigma.5)%*% matrix( rnorm(N*M), N,M)
Y<- F.true + sigma* matrix( rnorm(N*M), N,M)
# find MLE for lambda with grid of ranges
# and smoothness fixed in Matern
par.grid<- list( theta= seq( .1,.35,,8))
obj1b<- mKrigMLEGrid( x,Y,
cov.args = list(Covariance="Matern", smoothness=1.0),
par.grid = par.grid
)
obj$summary # take a look
# profile over theta
plot( par.grid$theta, obj1b$summary[,"lnProfileLike.FULL"],
type="b", log="x")
# }
# NOT RUN {
# }
# NOT RUN {
# 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(theta =.5, lambda = .5),
mKrig.args= list( m=0))
test.for.zero( obj1a$summary["sigmaMLE"], sigma, tol=.0075)
test.for.zero( obj1a$summary["theta"], theta, tol=.05)
# }
# NOT RUN {
##########################################################################
# A bootstrap example
# Here is a example of a more efficient (but less robust) bootstrap using
# mKrigMLEJoint and tuned starting values
##########################################################################
# }
# NOT RUN {
data( ozone2)
obj<- spatialProcess( ozone2$lon.lat,ozone2$y[16,] )
######### boot strap
set.seed(123)
M<- 50
# create M indepedent copies of the observation vector
ySynthetic<- simSpatialData( obj, M)
bootSummary<- NULL
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.fun = obj$cov.function.name,
cov.args = obj$cov.args,
cov.params.start = list( theta = obj$theta.MLE,
lambda = obj$lambda.MLE)
)
newSummary<- out$summary
bootSummary<- rbind( bootSummary, newSummary)
}
cat( " ", fill=TRUE )
obj$summary
stats( bootSummary)
# }
# NOT RUN {
#perform joint likelihood maximization over lambda, theta, 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( theta = .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(theta = .18,
smoothness = 1.1,
lambda = .08),
, REML=TRUE)
obj3$summary
# }
# NOT RUN {
#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(theta=.2, lambda=.1),
REML=TRUE)
#look at lnLikelihood evaluations
obj4$summary
# nails it!
# }
Run the code above in your browser using DataLab