************************** Code for Problem 1 ****************************** library(gstat) run = 1000 delta = 1 # 1.5 or 2 # step 1 # 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# res.mat1 <- matrix(0, nrow = 14, ncol = run) for(i in 1:run) { res.mat1[, i] <- variogram(simdata[[i]] ~ 1, data = simdata)$gamma } # step 3 & 4# res.mat2 <- matrix(0, nrow = 15, ncol = run) for(j in 1:run) { xy.jittered <- xy xy.jittered[, 1] <- xy.jittered[, 1] + runif(400, min = -delta, max = delta) xy.jittered[, 2] <- xy.jittered[, 2] + runif(400, min = -delta, max = delta) names(xy.jittered) <- c("x.jittered", "y.jittered") simdata2 <- data.frame(xy.jittered, simdata[[j]]) colnames(simdata2) <- c("x.jittered", "y.jittered", "z") coordinates(simdata2) = ~ x.jittered + y.jittered res.mat2[, j] <- variogram(simdata2$z ~ 1, data = simdata2)$gamma } # step 5 # par(mfrow = c(1, 2)) boxplot(t(res.mat1), ylim = c(0, max(res.mat1, res.mat2) + 0.1), main = "Original Locations", xlab = "Lag", ylab = "Sample Semivariogram") boxplot(t(res.mat2), ylim = c(0, max(res.mat1, res.mat2) + 0.1), main = "Jittered Locations", xlab = "Lag", ylab = "Sample SemiVariogram")