# NOT RUN {
## 10 member ensemble for the O3 data
# }
# NOT RUN {
data( "ozone2")
mKrigObject<- spatialProcess( ozone2$lon.lat, ozone2$y[16,],
smoothness=.5)
nx<- 65
ny<- 65
xGridList<- fields.x.to.grid( mKrigObject$x, nx=nx, ny=ny)
xGrid<- make.surface.grid( xGridList)
allTime0<- system.time(
look0<- sim.spatialProcess(mKrigObject, xp= xGrid, M=10)
)
print( allTime0)
## Local simulation with extra refinement of the grid for embedding
## and same grid size for prediction
## this runs much faster compared to exact method above
## as nx, ny are increased e.g. nx= 128, ny=128 is dramatic difference
allTime<- system.time(
look<- simLocal.mKrig(mKrigObject, M=10,nx=nx, ny=ny,
gridRefinement = 3,
np=3)
)
print( allTime)
print( look$timing)
# }
# NOT RUN {
# }
# NOT RUN {
## A simple example for setting up a bootstrap
## M below should be
## set to much larger sample size ( e.g. M <- 200) for better
## statistics
data( ozone2)
obj<- spatialProcess( ozone2$lon.lat,ozone2$y[16,] )
aHat<- obj$summary["aRange"]
lambdaHat<- obj$summary["lambda"]
######### boot strap
# create M independent copies of the observation vector
# here we just grab the model information from the
# spatialProcess object above.
#
# However, one could just create the list
# obj<- list( x= ozone2$lon.lat,
# cov.function.name="stationary.cov",
# summary= c( tau= 9.47, sigma2= 499.79, aRange= .700),
# cov.args= list( Covariance="Matern", smoothness=1.0),
# weights= rep( 1, nrow(ozone2$lon.lat) )
# )
# Here summary component has the parameters
# tau, sigma2 and aRange
# and cov.args component has the remaining ones.
set.seed(223)
M<- 25
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
newSummary<- spatialProcess(obj$x,ySynthetic[,k],
cov.params.start= list(
aRange = aHat,
lambda = lambdaHat)
)$summary
bootSummary<- rbind( bootSummary, newSummary)
}
cat( fill= TRUE)
# the results and 95<!-- % confidence interval -->
stats( bootSummary )
obj$summary
tmpBoot<- bootSummary[,c("lambda", "aRange") ]
confidenceInterval <- apply(tmpBoot, 2,
quantile, probs=c(0.025,0.975) )
# compare to estimates used as the "true" parameters
obj$summary[2:5]
print( t(confidenceInterval) )
# compare to confidence interval using large sample theory
print( obj$CITable)
# }
# NOT RUN {
# }
# NOT RUN {
# conditional simulation with covariates
# colorado climate example
data(COmonthlyMet)
fit1E<- spatialProcess(CO.loc,CO.tmin.MAM.climate, Z=CO.elev )
# conditional simulation at missing data
good<- !is.na(CO.tmin.MAM.climate )
infill<- sim.spatialProcess( fit1E, xp=CO.loc[!good,],
Z= CO.elev[!good], M= 10)
# get an elevation grid ... NGRID<- 50 gives a nicer image but takes longer
NGRID <- 25
# get elevations on a grid
COGrid<- list( x=seq( -109.5, -101, ,NGRID), y= seq(39, 41.5,,NGRID) )
COGridPoints<- make.surface.grid( COGrid)
# elevations are a bilinear interpolation from the 4km
# Rocky Mountain elevation fields data set.
data( RMelevation)
COElevGrid<- interp.surface( RMelevation, COGridPoints )
# NOTE call to sim.spatialProcess treats the grid points as just a matrix
# of locations the plot has to "reshape" these into a grid
# to use with image.plot
SEout<- sim.spatialProcess( fit1E, xp=COGridPoints, Z= COElevGrid, M= 30)
# for just the smooth surface in lon/lat
# SEout<- sim.spatialProcess( fit1E, xp=COGridPoints, drop.Z=TRUE, M= 30)
# in practice M should be larger to reduce Monte Carlo error.
surSE<- apply( SEout, 2, sd )
image.plot( as.surface( COGridPoints, surSE))
points( fit1E$x, col="magenta", pch=16)
# }
# NOT RUN {
data( ozone2)
set.seed( 399)
# fit to day 16 from Midwest ozone data set.
out<- Krig( ozone2$lon.lat, ozone2$y[16,], Covariance="Matern",
aRange=1.0,smoothness=1.0, na.rm=TRUE)
# NOTE aRange =1.0 is not the best choice but
# allows the sim.rf circulant embedding algorithm to
# work without increasing the domain.
#six missing data locations
xp<- ozone2$lon.lat[ is.na(ozone2$y[16,]),]
# 5 draws from process at xp given the data
# this is an exact calculation
sim.Krig( out,xp, M=5)-> sim.out
# Compare: stats(sim.out)[3,] to Exact: predictSE( out, xp)
# simulations on a grid
# NOTE this is approximate due to the bilinear interpolation
# for simulating the unconditional random field.
# also more grids points ( nx and ny) should be used
sim.Krig.approx(out,M=5, nx=20,ny=20)-> sim.out
# take a look at the ensemble members.
predictSurface( out, grid= list( x=sim.out$x, y=sim.out$y))-> look
zr<- c( 40, 200)
set.panel( 3,2)
image.plot( look, zlim=zr)
title("mean surface")
for ( k in 1:5){
image( sim.out$x, sim.out$y, sim.out$z[,,k], col=tim.colors(), zlim =zr)
}
# }
# NOT RUN {
data( ozone2)
y<- ozone2$y[16,]
good<- !is.na( y)
y<-y[good]
x<- ozone2$lon.lat[good,]
O3.fit<- mKrig( x,y, Covariance="Matern", aRange=.5,smoothness=1.0, lambda= .01 )
set.seed(122)
O3.sim<- sim.mKrig.approx( O3.fit, nx=100, ny=100, gridRefinement=3, M=5 )
set.panel(3,2)
surface( O3.fit)
for ( k in 1:5){
image.plot( as.surface( O3.sim$predictionPoints, O3.sim$Ensemble[,k]) )
}
# conditional simulation at missing data
xMissing<- ozone2$lon.lat[!good,]
O3.sim2<- sim.mKrig.approx( O3.fit, xMissing, nx=80, ny=80,
gridRefinement=3, M=4 )
# }
# NOT RUN {
#An example for fastTps:
data(ozone2)
y<- ozone2$y[16,]
good<- !is.na( y)
y<-y[good]
x<- ozone2$lon.lat[good,]
O3Obj<- fastTps( x,y, aRange=1.5 )
# creating a quick grid list based on ranges of locations
grid.list<- fields.x.to.grid( O3Obj$x, nx=100, ny=100)
# controlling the grids
xR<- range( x[,1], na.rm=TRUE)
yR<- range( x[,2], na.rm=TRUE)
simulationGridList<- list( x= seq(xR[1],xR[2],,400),
y= seq( yR[1],yR[2], ,400))
# very fine localized prediction grid
O3GridList<- list( x= seq( -90.5,-88.5,,200), y= seq( 38,40,,200))
O3Sim<- sim.mKrig.approx( O3Obj, M=5, predictionPointsList=O3GridList,
simulationGridList = simulationGridList)
# }
Run the code above in your browser using DataLab