swh:1:snp:33a53053e50f7abe7d281cc0c803be827debf4a3
Tip revision: a4ba0967b3c7f05b70d7e3f625228786cc4290e5 authored by Edzer Pebesma on 06 April 2023, 09:32:40 UTC
version 2.1-1
version 2.1-1
Tip revision: a4ba096
lhs.R
library(sp)
library(gstat)
# roughly follows the case presented in:
# E.J. Pebesma and G.B.M. Heuvelink, 1999. Latin hypercube sampling
# of Gaussian random fields. Technometrics 41 (4), pp. 303-312.
data(meuse)
data(meuse.grid)
coordinates(meuse) = ~x+y
gridded(meuse.grid) = ~x+y
x <- variogram(log(zinc) ~ 1, meuse)
v <- vgm(.5, "Sph", 800, nug = .1)
v.fit = fit.variogram(x, model = v)
plot(x, model = v.fit)
n = 100
ok = krige(log(zinc)~1, meuse, meuse.grid, v.fit, nmax=40)
sim = krige(log(zinc)~1, meuse, meuse.grid, v.fit, nsim = n, nmax=40)
simo = t(apply(as.data.frame(sim)[1:n], 1, order)) # rank order
nr = nrow(simo) # number of prediction locations
simo = (simo - 1.0 + matrix(runif(n * nr), nr, n))/n
summary(simo) # LHS on uniform [0,1] distribution; back to Gaussian:
lhs = t(apply(cbind(ok$var1.pred, sqrt(ok$var1.var), simo), 1,
function(x) qnorm(x[-(1:2)], x[1], x[2])))
sim2 = sim
sim2@data = data.frame(lhs)
spplot(sim2[1:10], main = 'lhs', col.regions=bpy.colors())
# verify that simulated and true mean/var are close:
m = apply(lhs, 1, mean)
v = apply(lhs, 1, var)
summary(m - ok$var1.pred)
summary(v - ok$var1.var)