library(gstat) library(sp) #### problem 5 # step 1# run = 1000 set.seed(22) xy <- expand.grid(1:20,1:20) names(xy) <- c("x", "y") g.dummy <- gstat(formula = z ~ 1, location = ~ x + y, dummy = T, beta = 0, model = vgm(psill = 1, model = "Exp", range = 3)) simdata <- predict(g.dummy, newdata = xy, nsim = run) coordinates(simdata) = ~ x + y # step 2 - 4# vgm.mat <- matrix(0, nrow = 20, ncol = run) for(i in 1:run) { vgm.mat[, i] <- variogram(simdata[[i]] ~ 1, data = simdata, width = 1, cutoff = 20)$gamma } cor(t(vgm.mat)) # step 5 # covgm.mat <- matrix(0, nrow = 20, ncol = run) for(i in 1:run) { covgm.mat[, i] <- variogram(simdata[[i]] ~ 1, data = simdata, width = 1, cutoff = 20, covariogram = TRUE)$gamma[-21] } cor(t(covgm.mat)) # boxplots # boxplot(t(vgm.mat), ylim = c(0, max(vgm.mat) + 0.1), xlab = "Lag", ylab = "Sample Semivariogram", main = "Semivariogram") boxplot(t(covgm.mat), ylim = c(min(covgm.mat) - 0.1, max(covgm.mat) + 0.1), xlab = "Lag", ylab = "Sample Covariogram", main = "Covariogram")