Last chance! 50% off unlimited learning
Sale ends in
sim.Krig.standard(object, xp, M = 1, verbose = FALSE, sigma2 = NA, rho = NA)sim.Krig.grid(object, grid.list = NA, M = 1, nx = 40, ny = 40, xy=c(1,2), verbose =
FALSE,
sigma2 = NA, rho = NA, extrap = FALSE)
sim.Krig.standard
a matrix with columns indexed by the locations in xp
and
M
rows.For sim.Krig.grid
a list with arguments x
and y
defining the grid
locations in the usual manner and z
contains the values of the simulated conditional
field(s). z
is a three dimesional array where the first two indices are "x" and "y" and
the third index is between 1 and M and indexes the simulated fields.
Given a spatial process Z(x)= P(x) + h(x) observed at
Y.k = P(x.k) + h(x.k) + e.k
where P(x) is a low order, fixed polynomial and h(x) a Gaussian spatial process. With Y= Y.1, ..., Y.N, the goal is to sample the conditional distribution of the process. [Z(x) | Y ]
For fixed a covariance this is just a multivariate normal sampling problem.
sim.Krig.standard
samples this conditional process at the points
xp
and is exact for fixed covariance parameters.
sim.Krig.grid
also assumes fixed covariance parameters and does
approxiamte sampling on a grid.
The outline of the algorithm is
0) Find the spatial prediction at the unobserved locations based on the actual data. Call this Z.hat(x).
1) Generate an unconditional spatial process and from this process simluate synthetic observations.
2) Use the spatial prediction model ( using the true covariance) to estimate the spatial process at unobserved locations.
3) Find the difference between the simulated process and its prediction based on synthetic observations. Call this e(x).
4) Z.hat(x) + e(x) is a draw from [Z(x) | Y ].
sim.Krig.standard
follows this algorithm exactly.
sim.Krig.grid
evaluates the
conditional surface on grid and simulates the values of h(x) off the grid using bilinear
interpolation of the four nearest grid points. Because of this approximation it is important to
choose the grid to be fine relative to the spacing of the observations. The advantage of this
approximation is that one can consider conditional simulation for large grids -- beyond the
size possible with exact methods. Here the method for simulation is circulant embedding and
so is restricted to correlation stationary fields.
data( ozone2)
set.seed( 399)
# fit to day 16 from Midwest ozone data set.
Krig( ozone2$lon.lat, ozone2$y[16,], Covariance="Matern",
theta=1.0,smoothness=1.0, na.rm=TRUE)-> out
# NOTE theta =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,]),]
# 50 draws from process at xp given the data
# this is an exact calculation
sim.Krig.standard( out,xp, M=50)-> sim.out
# Compare: stats(sim.out)[3,] to Exact: predict.se( out, xp)
# simulations on a grid
# NOTE this is approximate due to the bilinear interpolation
# for simulating the unconditional random field.
sim.Krig.grid(out,M=5)-> sim.out
# take a look at the ensemble members.
predict.surface( 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)
}
Run the code above in your browser using DataLab