#
# examples with doMKrig==TRUE
#
#generate observation locations
n=200
x = matrix(runif(2*n), nrow=n)
#generate observations at the locations
trueTheta = .2
trueLambda = .1
Sigma = exp( -rdist(x,x) /trueTheta )
# y = t(chol(Sigma))%*% (rnorm(n)) + trueLambda* rnorm( n)
y = t(chol(Sigma))%*% (rnorm(n)) + trueLambda* rnorm( n)
#Use exponential covariance, assume the true range parameter is known
out = MLESpatialProcess(x, y,
cov.args=list(Covariance="Exponential", range=trueTheta),
doMKrig=TRUE)
#Use exponential covariance, use a range to determine MLE of range parameter
## Not run:
# testThetas = seq(from=trueTheta/2, to=2*trueTheta, length=6)
# par.grid=list(theta=testThetas)
# out = MLESpatialProcess(x, y,
# cov.args=list(Covariance="Exponential"),
# par.grid=par.grid, doMKrig=TRUE)
#
# #Use exponential covariance, use a range to determine MLE of range
# #parameter, set custom lambda.grid
# testLambdas= seq(from=trueLambda/2, to=2*trueLambda, length=6)
# out = MLESpatialProcess(x, y, cov.args=list(Covariance="Exponential"),
# lambda.grid=testLambdas, par.grid=par.grid, doMKrig=TRUE)
#
# #Use Matern covariance, compute joint MLE of range, smoothness, and lambda.
# #This may take a few seconds
# testSmoothness = c(.5, 1, 2)
# par.grid=list(range=testThetas, smoothness=testSmoothness)
# out = MLESpatialProcess(x, y, cov.args=list(Covariance="Matern"),
# par.grid=par.grid, doMKrig=TRUE)
# ## End(Not run)
#
# examples with doMKrig==FALSE
#
N<- 100
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
# F.true<- t( Sigma.5)%*% rnorm(N)
F.true<- t( Sigma.5) %*% rnorm(N)
Y<- F.true + sigma*rnorm(N)
# find MLE for sigma rho and theta smoothness fixed first
# data set
obj<- MLESpatialProcess( x,Y)
obj$pars
# profile likelihood over theta
plot(obj$eval.grid[,1], obj$eval.grid[,6], xlab="theta",
ylab= "log Profile likelihood", type="p" )
xline( obj$pars["theta"], col="red")
# log likelihood surface over theta and log lambda
image.plot( obj$logLikelihoodSurface$x,
obj$logLikelihoodSurface$y, obj$logLikelihoodSurface$z,
xlab="theta (range)", ylab="log lambda" )
# MLE
points( obj$pars[1], log(obj$pars[2]), pch=16, col="magenta", cex=1.2)
# using "fast" version
obj.fast<- MLESpatialProcess.fast( x,Y)
obj.fast$pars
# points where likelihood evaluated:
quilt.plot( log( obj.fast$eval.grid[,1:2] ), obj.fast$eval.grid[,7],
xlab="log(theta)",ylab="log(lambda)")
# parameters are slightly different due to the differences of REML and the full likelihood.
# example with a covariate
## Not run:
# data(COmonthlyMet)
# obj2<- MLESpatialProcess( CO.loc, CO.tmean.MAM.climate)
# obj3<- MLESpatialProcess( CO.loc, CO.tmean.MAM.climate, Z= CO.elev)
# ind<- !is.na( CO.tmean.MAM.climate)
# obj4<- MLESpatialProcess.fast( CO.loc[ind,], CO.tmean.MAM.climate[ind],
# Z= CO.elev[ind])
# # elevation makes a difference
# obj2$pars
# obj3$pars
# obj4$pars
# ## End(Not run)
## Not run:
# # fits for ozone data
# data( ozone2)
# NDays<- nrow( ozone2$y)
# O3MLE<- matrix( NA, nrow= NDays, ncol=4)
# dimnames(O3MLE)<- list(NULL, c("theta", "lambda", "rho", "sigma"))
# for( day in 1: NDays){
# cat( day, " ")
# O3MLE[day,]<- MLESpatialProcess( ozone2$lon.lat, ozone2$y[day,],
# Distance="rdist.earth")$pars
# }
# plot( log(O3MLE[,1]), log(O3MLE[,3]))
# ## End(Not run)
Run the code above in your browser using DataLab