### R code from vignette source 'SpatialExtremesGuide.Rnw' ################################################### ### code chunk number 1: SpatialExtremesGuide.Rnw:33-43 ################################################### library(SpatialExtremes) set.seed(12) n.site <- 50 coord <- matrix(runif(2*n.site, 0, 10), ncol = 2) ##Simulate a max-stable process - with unit Frechet margins data <- rmaxstab(40, coord, "whitmat", nugget = 0, range = 1, smooth = 2) ##Compute the lambda-madogram lmadogram(data, coord, n.bins = 25, border = "grey", box = FALSE) ################################################### ### code chunk number 2: SpatialExtremesGuide.Rnw:57-59 ################################################### options(SweaveHooks=list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1)))) ################################################### ### code chunk number 3: SpatialExtremesGuide.Rnw:100-101 (eval = FALSE) ################################################### ## vignette("SpatialExtremesGuide") ################################################### ### code chunk number 4: Smith2Sim (eval = FALSE) ################################################### ## x <- seq(0, 10, length = 100) ## seed <- 19 ## set.seed(seed) ## data0 <- rmaxstab(1, cbind(x,x), "gauss", cov11 = 9/8, cov12 = 0, cov22 = 9/8, ## grid = TRUE) ## set.seed(seed) ## data1 <- rmaxstab(1, cbind(x,x), "gauss", cov11 = 9/8, cov12 = 1, cov22 = 9/8, ## grid = TRUE) ## ## png("Figures/Smith2Sim.png", width = 1400, height = 700) ## par(mfrow=c(1,2)) ## image(x, x, log(data0), col = terrain.colors(64)) ## image(x, x, log(data1), col = terrain.colors(64)) ## dev.off() ################################################### ### code chunk number 5: covariances ################################################### par(mfrow=c(1,4), mar = c(4, 4, 0.2, 0.2)) covariance(nugget = 0, sill = 1, range = 1, smooth = 4, cov.mod = "whitmat", xlim = c(0,9), ylim = c(0, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "whitmat", add = TRUE, col = 2, xlim = c(0,9)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", add = TRUE, col = 3, xlim = c(0,9)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.5, cov.mod = "whitmat", col = 4, add = TRUE, xlim = c(0,9)) legend("topright", c(expression(nu == 4), expression(nu == 2), expression(nu == 1), expression(nu == 0.5)), col = 1:4, lty = 1, inset = 0.05) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "powexp", xlim = c(0, 4), ylim = c(0, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1.5, cov.mod = "powexp", add = TRUE, col = 2, xlim = c(0, 4)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "powexp", add = TRUE, col = 3, xlim = c(0, 4)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.75, cov.mod = "powexp", add = TRUE, col = 4, xlim = c(0, 4)) legend("topright", c(expression(nu == 2), expression(nu == 1.5), expression(nu == 1), expression(nu == 0.75)), col = 1:4, lty = 1, inset = 0.05) covariance(nugget = 0, sill = 1, range = 1, smooth = 4, cov.mod = "cauchy", xlim = c(0, 5), ylim = c(0, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "cauchy", add = TRUE, col = 2, xlim = c(0, 5)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "cauchy", add = TRUE, col = 3, xlim = c(0, 5)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.75, cov.mod = "cauchy", add = TRUE, col = 4, xlim = c(0, 5)) legend("topright", c(expression(nu == 4), expression(nu == 2), expression(nu == 1), expression(nu == 0.75)), col = 1:4, lty = 1, inset = 0.05) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "bessel", xlim = c(0, 20), ylim = c(-0.3, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "bessel", add = TRUE, col = 2, xlim = c(0, 20)) covariance(nugget = 0, sill = 1, range = 1, smooth = .5, cov.mod = "bessel", add = TRUE, col = 3, xlim = c(0, 20)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.25, cov.mod = "bessel", add = TRUE, col = 4, xlim = c(0, 20)) legend("topright", c(expression(nu == 2), expression(nu == 1), expression(nu == .5), expression(nu == 0.25)), col = 1:4, lty = 1, inset = 0.05) ################################################### ### code chunk number 6: SpatialExtremesGuide.Rnw:556-557 ################################################### par(mfrow=c(1,4), mar = c(4, 4, 0.2, 0.2)) covariance(nugget = 0, sill = 1, range = 1, smooth = 4, cov.mod = "whitmat", xlim = c(0,9), ylim = c(0, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "whitmat", add = TRUE, col = 2, xlim = c(0,9)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", add = TRUE, col = 3, xlim = c(0,9)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.5, cov.mod = "whitmat", col = 4, add = TRUE, xlim = c(0,9)) legend("topright", c(expression(nu == 4), expression(nu == 2), expression(nu == 1), expression(nu == 0.5)), col = 1:4, lty = 1, inset = 0.05) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "powexp", xlim = c(0, 4), ylim = c(0, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1.5, cov.mod = "powexp", add = TRUE, col = 2, xlim = c(0, 4)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "powexp", add = TRUE, col = 3, xlim = c(0, 4)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.75, cov.mod = "powexp", add = TRUE, col = 4, xlim = c(0, 4)) legend("topright", c(expression(nu == 2), expression(nu == 1.5), expression(nu == 1), expression(nu == 0.75)), col = 1:4, lty = 1, inset = 0.05) covariance(nugget = 0, sill = 1, range = 1, smooth = 4, cov.mod = "cauchy", xlim = c(0, 5), ylim = c(0, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "cauchy", add = TRUE, col = 2, xlim = c(0, 5)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "cauchy", add = TRUE, col = 3, xlim = c(0, 5)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.75, cov.mod = "cauchy", add = TRUE, col = 4, xlim = c(0, 5)) legend("topright", c(expression(nu == 4), expression(nu == 2), expression(nu == 1), expression(nu == 0.75)), col = 1:4, lty = 1, inset = 0.05) covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "bessel", xlim = c(0, 20), ylim = c(-0.3, 1)) covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "bessel", add = TRUE, col = 2, xlim = c(0, 20)) covariance(nugget = 0, sill = 1, range = 1, smooth = .5, cov.mod = "bessel", add = TRUE, col = 3, xlim = c(0, 20)) covariance(nugget = 0, sill = 1, range = 1, smooth = 0.25, cov.mod = "bessel", add = TRUE, col = 4, xlim = c(0, 20)) legend("topright", c(expression(nu == 2), expression(nu == 1), expression(nu == .5), expression(nu == 0.25)), col = 1:4, lty = 1, inset = 0.05) ################################################### ### code chunk number 7: SpatialExtremesGuide.Rnw:610-621 (eval = FALSE) ################################################### ## covariance(nugget = 0, sill = 1, range = 1, smooth = 4, cov.mod = "whitmat", ## xlim = c(0,9), ylim = c(0, 1)) ## covariance(nugget = 0, sill = 1, range = 1, smooth = 2, cov.mod = "whitmat", ## add = TRUE, col = 2, xlim = c(0,9)) ## covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", ## add = TRUE, col = 3, xlim = c(0,9)) ## covariance(nugget = 0, sill = 1, range = 1, smooth = 0.5, cov.mod = "whitmat", ## col = 4, add = TRUE, xlim = c(0,9)) ## legend("topright", c(expression(nu == 4), expression(nu == 2), ## expression(nu == 1), expression(nu == 0.5)), ## col = 1:4, lty = 1, inset = 0.05) ################################################### ### code chunk number 8: Schlather2Sim (eval = FALSE) ################################################### ## x <- y <- seq(0, 10, length = 100) ## set.seed(12) ## ms0 <- rmaxstab(1, cbind(x, y), "whitmat", nugget = 0, range = 1, smooth = 1, grid=TRUE) ## set.seed(12) ## ms1 <- rmaxstab(1, cbind(x, y), "powexp", nugget = 0, range = 1.5, smooth = 1, grid=TRUE) ## ## png("Figures/Schlather2Sim.png", width = 1400, height = 700) ## par(mfrow=c(1,2)) ## image(x, y, log(ms0), col = terrain.colors(64)) ## image(x, y, log(ms1), col = terrain.colors(64)) ## dev.off() ################################################### ### code chunk number 9: anisoCovFun ################################################### cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "powexp", plot = FALSE) phi <- pi / 4 r <- sqrt(0.5) A <- matrix(c(cos(phi), r^2 * sin(phi), r^2 * sin(phi), cos(phi)), 2) cov.fun.aniso <- function(vec.pos) cov.fun(sqrt(vec.pos %*% A %*% vec.pos)) rho1 <- rho2 <- matrix(NA, 100, 100) xs <- ys <- seq(-4, 4, length = 100) for (i in 1:100){ x <- xs[i] for (j in 1:100){ rho1[i,j] <- cov.fun(sqrt(x^2+ys[j]^2)) rho2[i,j] <- cov.fun.aniso(c(x, ys[j])) } } par(mfrow=c(1,2)) contour(xs, ys, rho1, xlab = expression(paste(Delta, x)), ylab = expression(paste(Delta, y))) contour(xs, ys, rho2, xlab = expression(paste(Delta, x)), ylab = expression(paste(Delta, y))) ################################################### ### code chunk number 10: SpatialExtremesGuide.Rnw:686-687 ################################################### cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "powexp", plot = FALSE) phi <- pi / 4 r <- sqrt(0.5) A <- matrix(c(cos(phi), r^2 * sin(phi), r^2 * sin(phi), cos(phi)), 2) cov.fun.aniso <- function(vec.pos) cov.fun(sqrt(vec.pos %*% A %*% vec.pos)) rho1 <- rho2 <- matrix(NA, 100, 100) xs <- ys <- seq(-4, 4, length = 100) for (i in 1:100){ x <- xs[i] for (j in 1:100){ rho1[i,j] <- cov.fun(sqrt(x^2+ys[j]^2)) rho2[i,j] <- cov.fun.aniso(c(x, ys[j])) } } par(mfrow=c(1,2)) contour(xs, ys, rho1, xlab = expression(paste(Delta, x)), ylab = expression(paste(Delta, y))) contour(xs, ys, rho2, xlab = expression(paste(Delta, x)), ylab = expression(paste(Delta, y))) ################################################### ### code chunk number 11: tbm (eval = FALSE) ################################################### ## x <- y <- seq(0, 10, length = 100) ## coord <- cbind(x, y) ## seed <- 3 ## set.seed(seed) ## data1 <- rgp(1, coord, cov.mod = "whitmat", sill = 1, range = 1, smooth = 1, ## grid = TRUE, control = list(nlines = 1)) ## set.seed(seed) ## data2 <- rgp(1, coord, cov.mod = "whitmat", sill = 1, range = 1, smooth = 1, ## grid = TRUE, control = list(nlines = 25)) ## set.seed(seed) ## data3 <- rgp(1, coord, cov.mod = "whitmat", sill = 1, range = 1, smooth = 1, ## grid = TRUE, control = list(nlines = 625)) ## ## png("Figures/tbm.png", width = 1400, height = 700) ## par(mfrow=c(1,3)) ## image(x, y, data1, col = terrain.colors(64)) ## image(x, y, data2, col = terrain.colors(64)) ## image(x, y, data3, col = terrain.colors(64)) ## dev.off() ################################################### ### code chunk number 12: SpatialExtremesGuide.Rnw:1009-1010 (eval = FALSE) ################################################### ## x <- y <- seq(0, 10, length = 100) ## coord <- cbind(x, y) ## seed <- 3 ## set.seed(seed) ## data1 <- rgp(1, coord, cov.mod = "whitmat", sill = 1, range = 1, smooth = 1, ## grid = TRUE, control = list(nlines = 1)) ## set.seed(seed) ## data2 <- rgp(1, coord, cov.mod = "whitmat", sill = 1, range = 1, smooth = 1, ## grid = TRUE, control = list(nlines = 25)) ## set.seed(seed) ## data3 <- rgp(1, coord, cov.mod = "whitmat", sill = 1, range = 1, smooth = 1, ## grid = TRUE, control = list(nlines = 625)) ## ## png("Figures/tbm.png", width = 1400, height = 700) ## par(mfrow=c(1,3)) ## image(x, y, data1, col = terrain.colors(64)) ## image(x, y, data2, col = terrain.colors(64)) ## image(x, y, data3, col = terrain.colors(64)) ## dev.off() ################################################### ### code chunk number 13: rmaxstabSmith (eval = FALSE) ################################################### ## x <- y <- seq(0, 10, length = 100) ## coord <- cbind(x, y) ## set.seed(8) ## M0 <- rmaxstab(1, coord, "gauss", cov11 = 9/8, cov12 = 0, cov22 = 9/8, ## grid = TRUE) ## set.seed(8) ## M1 <- rmaxstab(1, coord, "gauss", cov11 = 9/8, cov12 = 3/4, cov22 = 9/8, ## grid = TRUE) ## ## png("Figures/rmaxstabSmith.png", width = 1400, height = 700) ## par(mfrow = c(1,2)) ## image(x, y, log(M0), col = terrain.colors(64)) ## image(x, y, log(M1), col = terrain.colors(64)) ## dev.off() ################################################### ### code chunk number 14: SpatialExtremesGuide.Rnw:1164-1165 (eval = FALSE) ################################################### ## x <- y <- seq(0, 10, length = 100) ## coord <- cbind(x, y) ## set.seed(8) ## M0 <- rmaxstab(1, coord, "gauss", cov11 = 9/8, cov12 = 0, cov22 = 9/8, ## grid = TRUE) ## set.seed(8) ## M1 <- rmaxstab(1, coord, "gauss", cov11 = 9/8, cov12 = 3/4, cov22 = 9/8, ## grid = TRUE) ## ## png("Figures/rmaxstabSmith.png", width = 1400, height = 700) ## par(mfrow = c(1,2)) ## image(x, y, log(M0), col = terrain.colors(64)) ## image(x, y, log(M1), col = terrain.colors(64)) ## dev.off() ################################################### ### code chunk number 15: SpatialExtremesGuide.Rnw:1191-1193 (eval = FALSE) ################################################### ## coord <- matrix(runif(100, 0, 10), ncol = 2) ## data1 <- rmaxstab(100, coord, "whitmat", nugget = 0, range = 1, smooth = 1) ################################################### ### code chunk number 16: SpatialExtremesGuide.Rnw:1196-1200 (eval = FALSE) ################################################### ## x <- seq(0, 10, length = 100) ## coord <- cbind(x, x) ## data2 <- rmaxstab(1, coord, "powexp", nugget = 0, range = 1, smooth = 2, ## grid = TRUE) ################################################### ### code chunk number 17: extCoeffModels ################################################### smith <- function(h) 2 * pnorm(h/2) cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, plot = FALSE) schlather <- function(h) 1 + sqrt((1-cov.fun(h))/2) alpha <- 0.1 plot(smith, from = 0, to = 5, xlab = "h", ylab = expression(theta(h))) plot(schlather, add = TRUE, col = 2, from = 0, to = 5) legend("bottomright", c("Smith", "Schlather"), col = 1:2, lty = 1, inset = 0.05) ################################################### ### code chunk number 18: SpatialExtremesGuide.Rnw:1305-1306 ################################################### smith <- function(h) 2 * pnorm(h/2) cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, plot = FALSE) schlather <- function(h) 1 + sqrt((1-cov.fun(h))/2) alpha <- 0.1 plot(smith, from = 0, to = 5, xlab = "h", ylab = expression(theta(h))) plot(schlather, add = TRUE, col = 2, from = 0, to = 5) legend("bottomright", c("Smith", "Schlather"), col = 1:2, lty = 1, inset = 0.05) ################################################### ### code chunk number 19: extCoeffST-Smith ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", plot = FALSE) ext.coeff <- function(h) 1 + sqrt((1 - cov.fun(h))/2) par(mfrow=c(1,2), ps=14) ST <- fitextcoeff(data, coord, loess = FALSE, ylim = c(1, 2)) plot(ext.coeff, from = 0, to = 12, col = 2, add = TRUE, lwd = 1.5) Smith <- fitextcoeff(data, coord, estim = "Smith", loess = FALSE, ylim = c(1, 2)) plot(ext.coeff, from = 0, to = 12, col = 2, add = TRUE, lwd = 1.5) ################################################### ### code chunk number 20: SpatialExtremesGuide.Rnw:1409-1410 ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", plot = FALSE) ext.coeff <- function(h) 1 + sqrt((1 - cov.fun(h))/2) par(mfrow=c(1,2), ps=14) ST <- fitextcoeff(data, coord, loess = FALSE, ylim = c(1, 2)) plot(ext.coeff, from = 0, to = 12, col = 2, add = TRUE, lwd = 1.5) Smith <- fitextcoeff(data, coord, estim = "Smith", loess = FALSE, ylim = c(1, 2)) plot(ext.coeff, from = 0, to = 12, col = 2, add = TRUE, lwd = 1.5) ################################################### ### code chunk number 21: SpatialExtremesGuide.Rnw:1426-1435 (eval = FALSE) ################################################### ## n.site <- 40 ## n.obs <- 100 ## coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) ## data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, ## smooth = 1) ## ## par(mfrow=c(1,2)) ## fitextcoeff(data, coord, loess = FALSE) ## fitextcoeff(data, coord, estim = "Smith", loess = FALSE) ################################################### ### code chunk number 22: madogram ################################################### cov.fun1 <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", plot = FALSE) ext.coeff1 <- function(h) 1 + sqrt((1 - cov.fun1(h))/2) mado1 <- function(h) log(ext.coeff1(h)) cov.fun2 <- covariance(nugget = 0, sill = 1, range = 1.5, smooth = 1, cov.mod = "powexp", plot = FALSE) ext.coeff2 <- function(h) 1 + sqrt((1 - cov.fun2(h))/2) mado2 <- function(h) log(ext.coeff2(h)) n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) par(mfrow=c(1,2)) madogram(data, coord, which = "mado", ylim = c(0, log(2))) plot(mado1, from = 0, to = 12, add = TRUE, col = 2, lwd = 1.5) madogram(data, coord, which = "mado", ylim = c(0, log(2)), n.bins = 100) plot(mado1, from = 0, to = 12, add = TRUE, col = 2, lwd = 1.5) ################################################### ### code chunk number 23: SpatialExtremesGuide.Rnw:1521-1522 ################################################### cov.fun1 <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", plot = FALSE) ext.coeff1 <- function(h) 1 + sqrt((1 - cov.fun1(h))/2) mado1 <- function(h) log(ext.coeff1(h)) cov.fun2 <- covariance(nugget = 0, sill = 1, range = 1.5, smooth = 1, cov.mod = "powexp", plot = FALSE) ext.coeff2 <- function(h) 1 + sqrt((1 - cov.fun2(h))/2) mado2 <- function(h) log(ext.coeff2(h)) n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) par(mfrow=c(1,2)) madogram(data, coord, which = "mado", ylim = c(0, log(2))) plot(mado1, from = 0, to = 12, add = TRUE, col = 2, lwd = 1.5) madogram(data, coord, which = "mado", ylim = c(0, log(2)), n.bins = 100) plot(mado1, from = 0, to = 12, add = TRUE, col = 2, lwd = 1.5) ################################################### ### code chunk number 24: SpatialExtremesGuide.Rnw:1554-1562 (eval = FALSE) ################################################### ## n.site <- 50 ## n.obs <- 100 ## coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) ## data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) ## ## par(mfrow=c(1,2)) ## madogram(data, coord, which = "mado") ## madogram(data, coord, which = "mado", n.bins = 100) ################################################### ### code chunk number 25: madogramExtCoeff ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", plot = FALSE) ext.coeff <- function(h) 1 + sqrt((1 - cov.fun(h))/2) mado <- function(h) log(ext.coeff(h)) par(mfrow=c(1,2)) madogram(data, coord, which = "mado", ylim = c(0, log(2))) plot(mado, from = 0, to = 12, col = 2, lwd = 1.5, add = TRUE) madogram(data, coord, which = "ext") plot(ext.coeff, from = 0, to = 12, col = 2, lwd = 1.5, add = TRUE) ################################################### ### code chunk number 26: SpatialExtremesGuide.Rnw:1599-1600 ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) cov.fun <- covariance(nugget = 0, sill = 1, range = 1, smooth = 1, cov.mod = "whitmat", plot = FALSE) ext.coeff <- function(h) 1 + sqrt((1 - cov.fun(h))/2) mado <- function(h) log(ext.coeff(h)) par(mfrow=c(1,2)) madogram(data, coord, which = "mado", ylim = c(0, log(2))) plot(mado, from = 0, to = 12, col = 2, lwd = 1.5, add = TRUE) madogram(data, coord, which = "ext") plot(ext.coeff, from = 0, to = 12, col = 2, lwd = 1.5, add = TRUE) ################################################### ### code chunk number 27: SpatialExtremesGuide.Rnw:1617-1623 (eval = FALSE) ################################################### ## n.site <- 40 ## n.obs <- 100 ## coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) ## data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0, range = 1, smooth = 1) ## ## madogram(data, coord) ################################################### ### code chunk number 28: F-madogram ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "gauss", cov11 = 1, cov12 = 0, cov22 = 1) ext.coeff <- function(h) 2 * pnorm(h/2) Fmado <- function(h) 0.5 * (ext.coeff(h) - 1) / (ext.coeff(h) + 1) par(mfrow=c(1,2)) fmadogram(data, coord, which = "mado", ylim = c(0, 1/5)) plot(Fmado, from = 0, to = 13, add = TRUE, col = 2, lwd = 1.5) fmadogram(data, coord, which = "ext") plot(ext.coeff, from = 0, to = 13, add = TRUE, col = 2, lwd = 1.5) ################################################### ### code chunk number 29: SpatialExtremesGuide.Rnw:1662-1663 ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "gauss", cov11 = 1, cov12 = 0, cov22 = 1) ext.coeff <- function(h) 2 * pnorm(h/2) Fmado <- function(h) 0.5 * (ext.coeff(h) - 1) / (ext.coeff(h) + 1) par(mfrow=c(1,2)) fmadogram(data, coord, which = "mado", ylim = c(0, 1/5)) plot(Fmado, from = 0, to = 13, add = TRUE, col = 2, lwd = 1.5) fmadogram(data, coord, which = "ext") plot(ext.coeff, from = 0, to = 13, add = TRUE, col = 2, lwd = 1.5) ################################################### ### code chunk number 30: SpatialExtremesGuide.Rnw:1746-1753 (eval = FALSE) ################################################### ## n.site <- 40 ## n.obs <- 100 ## coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) ## data <- rmaxstab(n.obs, coord, "gauss", cov11 = 1, cov12 = 0, cov22 = 1) ## ## par(mfrow=c(1,2)) ## fmadogram(data, coord) ################################################### ### code chunk number 31: lambda-madogram ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data1 <- rmaxstab(n.obs, coord, "powexp", nugget = 0, range = 1, smooth = 1) data2 <- rmaxstab(n.obs, coord, "cauchy", nugget = 0, range = 1, smooth = 1) par(mfrow=c(1,2), pty = "s") lmadogram(data1, coord, n.bins = 60) lmadogram(data2, coord, n.bins = 60) ################################################### ### code chunk number 32: SpatialExtremesGuide.Rnw:1901-1902 ################################################### n.site <- 40 n.obs <- 100 coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) data1 <- rmaxstab(n.obs, coord, "powexp", nugget = 0, range = 1, smooth = 1) data2 <- rmaxstab(n.obs, coord, "cauchy", nugget = 0, range = 1, smooth = 1) par(mfrow=c(1,2), pty = "s") lmadogram(data1, coord, n.bins = 60) lmadogram(data2, coord, n.bins = 60) ################################################### ### code chunk number 33: SpatialExtremesGuide.Rnw:1916-1925 (eval = FALSE) ################################################### ## n.site <- 40 ## n.obs <- 100 ## coord <- matrix(runif(2 * n.site, 0, 10), ncol = 2) ## data1 <- rmaxstab(n.obs, coord, "powexp", nugget = 0, range = 1, smooth = 1) ## data2 <- rmaxstab(n.obs, coord, "cauchy", nugget = 0, range = 1, smooth = 1) ## ## par(mfrow=c(1,2), pty = "s") ## lmadogram(data1, coord, n.bins = 60) ## lmadogram(data2, coord, n.bins = 60) ################################################### ### code chunk number 34: SpatialExtremesGuide.Rnw:1979-1985 ################################################### n.site <- 40 n.obs <- 80 coord <- matrix(runif(2*n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "gauss", cov11 = 9/8, cov12 = 1, cov22 = 9/8) fitcovmat(data, coord, marge = "emp") ################################################### ### code chunk number 35: SpatialExtremesGuide.Rnw:2120-2127 ################################################### n.obs <- 80 n.site <- 40 set.seed(12) coord <- matrix(runif(2*n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "whitmat", nugget = 0.8, range = 3, smooth = 1.2) fitmaxstab(data, coord, cov.mod = "whitmat") ################################################### ### code chunk number 36: SpatialExtremesGuide.Rnw:2148-2151 (eval = FALSE) ################################################### ## fitmaxstab(data, coord, cov.mod = "whitmat", smooth = 1.2) ## fitmaxstab(data, coord, cov.mod = "whitmat", nugget = 0) ## fitmaxstab(data, coord, cov.mod = "whitmat", range = 3) ################################################### ### code chunk number 37: SpatialExtremesGuide.Rnw:2157-2159 (eval = FALSE) ################################################### ## fitmaxstab(data, coord, cov.mod = "cauchy") ## fitmaxstab(data, coord, cov.mod = "powexp") ################################################### ### code chunk number 38: SpatialExtremesGuide.Rnw:2163-2164 (eval = FALSE) ################################################### ## fitmaxstab(data, coord, cov.mod = "gauss") ################################################### ### code chunk number 39: SpatialExtremesGuide.Rnw:2170-2171 (eval = FALSE) ################################################### ## fitmaxstab(data, coord, cov.mod = "gauss", cov12 = 0, method = "BFGS") ################################################### ### code chunk number 40: SpatialExtremesGuide.Rnw:2343-2344 ################################################### set.seed(1) ################################################### ### code chunk number 41: SpatialExtremesGuide.Rnw:2347-2355 ################################################### n.obs <- 80 n.site <- 40 coord <- matrix(runif(2*n.site, 0, 10), ncol = 2) data <- rmaxstab(n.obs, coord, "cauchy", nugget = 0.2, range = 3, smooth = 1.2) M0 <- fitmaxstab(data, coord, cov.mod = "powexp") M1 <- fitmaxstab(data, coord, cov.mod = "cauchy") TIC(M0, M1) ################################################### ### code chunk number 42: SpatialExtremesGuide.Rnw:2421-2422 ################################################### set.seed(7) ################################################### ### code chunk number 43: SpatialExtremesGuide.Rnw:2425-2433 ################################################### n.obs <- 50 n.site <- 30 coord <- matrix(rnorm(2*n.site, sd = sqrt(.2)), ncol = 2) data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220) M0 <- fitmaxstab(data, coord, "gauss", cov11 = 100) M1 <- fitmaxstab(data, coord, "gauss") anova(M0, M1) ################################################### ### code chunk number 44: SpatialExtremesGuide.Rnw:2511-2512 ################################################### anova(M0, M1, method = "CB") ################################################### ### code chunk number 45: SpatialExtremesGuide.Rnw:2547-2550 ################################################### x <- c(2.2975896, 1.6448808, 1.3323833, -0.4464904, 2.2737603, -0.2581876, 9.5184398, -0.5899699, 0.4974283, -0.8152157) z <- gev2frech(x, 1, 2, .2) ################################################### ### code chunk number 46: SpatialExtremesGuide.Rnw:2561-2562 ################################################### frech2gev(z, 1, 2, .2) ################################################### ### code chunk number 47: SpatialExtremesGuide.Rnw:2619-2620 (eval = FALSE) ################################################### ## fitmaxstab(data, coord, "gauss", fit.marge = TRUE) ################################################### ### code chunk number 48: knotsEffect ################################################### set.seed(12) x <- runif(100) fun <- function(x) sin(3 * pi * x) y <- fun(x) + rnorm(100, 0, 0.15) knots1 <- quantile(x, prob = 1:2 / 3) knots2 <- quantile(x, prob = 1:10 / 11) knots3 <- quantile(x, prob = 1:50 / 51) M0 <- rbpspline(y, x, knots = knots1, degree = 3, penalty = 0) M1 <- rbpspline(y, x, knots = knots2, degree = 3, penalty = 0) M2 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0) par(mfrow=c(1,3)) plot(x, y, col = "lightgrey") rug(knots1) lines(M0) plot(x, y, col = "lightgrey") rug(knots2) lines(M1, col = 2) plot(x, y, col = "lightgrey") rug(knots3) lines(M2, col = 3) ################################################### ### code chunk number 49: SpatialExtremesGuide.Rnw:2854-2855 ################################################### set.seed(12) x <- runif(100) fun <- function(x) sin(3 * pi * x) y <- fun(x) + rnorm(100, 0, 0.15) knots1 <- quantile(x, prob = 1:2 / 3) knots2 <- quantile(x, prob = 1:10 / 11) knots3 <- quantile(x, prob = 1:50 / 51) M0 <- rbpspline(y, x, knots = knots1, degree = 3, penalty = 0) M1 <- rbpspline(y, x, knots = knots2, degree = 3, penalty = 0) M2 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0) par(mfrow=c(1,3)) plot(x, y, col = "lightgrey") rug(knots1) lines(M0) plot(x, y, col = "lightgrey") rug(knots2) lines(M1, col = 2) plot(x, y, col = "lightgrey") rug(knots3) lines(M2, col = 3) ################################################### ### code chunk number 50: SpatialExtremesGuide.Rnw:2871-2872 (eval = FALSE) ################################################### ## set.seed(12) ## x <- runif(100) ## fun <- function(x) sin(3 * pi * x) ## y <- fun(x) + rnorm(100, 0, 0.15) ## knots1 <- quantile(x, prob = 1:2 / 3) ## knots2 <- quantile(x, prob = 1:10 / 11) ## knots3 <- quantile(x, prob = 1:50 / 51) ## M0 <- rbpspline(y, x, knots = knots1, degree = 3, penalty = 0) ## M1 <- rbpspline(y, x, knots = knots2, degree = 3, penalty = 0) ## M2 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0) ## ## par(mfrow=c(1,3)) ## plot(x, y, col = "lightgrey") ## rug(knots1) ## lines(M0) ## plot(x, y, col = "lightgrey") ## rug(knots2) ## lines(M1, col = 2) ## plot(x, y, col = "lightgrey") ## rug(knots3) ## lines(M2, col = 3) ################################################### ### code chunk number 51: smoothingParamEffect ################################################### M0 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0) M1 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0.1) M2 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 10) par(mfrow=c(1,3)) plot(x, y, col = "lightgrey") lines(M0) plot(x, y, col = "lightgrey") lines(M1, col = 2) plot(x, y, col = "lightgrey") lines(M2, col = 3) ################################################### ### code chunk number 52: SpatialExtremesGuide.Rnw:2897-2898 ################################################### M0 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0) M1 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0.1) M2 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 10) par(mfrow=c(1,3)) plot(x, y, col = "lightgrey") lines(M0) plot(x, y, col = "lightgrey") lines(M1, col = 2) plot(x, y, col = "lightgrey") lines(M2, col = 3) ################################################### ### code chunk number 53: SpatialExtremesGuide.Rnw:2930-2931 (eval = FALSE) ################################################### ## M0 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0) ## M1 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 0.1) ## M2 <- rbpspline(y, x, knots = knots3, degree = 3, penalty = 10) ## ## par(mfrow=c(1,3)) ## plot(x, y, col = "lightgrey") ## lines(M0) ## plot(x, y, col = "lightgrey") ## lines(M1, col = 2) ## plot(x, y, col = "lightgrey") ## lines(M2, col = 3) ################################################### ### code chunk number 54: CVandGCV ################################################### par(mfrow=c(1,3)) lambda.cv <- cv(y, x, knots = knots3, degree = 3)$penalty abline(v = lambda.cv, lty = 2) lambda.gcv <- gcv(y, x, knots = knots3, degree = 3)$penalty abline(v = lambda.gcv, lty = 2) cv.fit <- rbpspline(y, x, knots3, degree = 3, penalty = "cv") gcv.fit <- rbpspline(y, x, knots3, degree = 3, penalty = "gcv") plot(x, y, col = "lightgrey") lines(cv.fit, col = 2) lines(gcv.fit, col = 3) ################################################### ### code chunk number 55: SpatialExtremesGuide.Rnw:3022-3023 ################################################### par(mfrow=c(1,3)) lambda.cv <- cv(y, x, knots = knots3, degree = 3)$penalty abline(v = lambda.cv, lty = 2) lambda.gcv <- gcv(y, x, knots = knots3, degree = 3)$penalty abline(v = lambda.gcv, lty = 2) cv.fit <- rbpspline(y, x, knots3, degree = 3, penalty = "cv") gcv.fit <- rbpspline(y, x, knots3, degree = 3, penalty = "gcv") plot(x, y, col = "lightgrey") lines(cv.fit, col = 2) lines(gcv.fit, col = 3) ################################################### ### code chunk number 56: SpatialExtremesGuide.Rnw:3038-3039 (eval = FALSE) ################################################### ## par(mfrow=c(1,3)) ## lambda.cv <- cv(y, x, knots = knots3, degree = 3)$penalty ## abline(v = lambda.cv, lty = 2) ## lambda.gcv <- gcv(y, x, knots = knots3, degree = 3)$penalty ## abline(v = lambda.gcv, lty = 2) ## cv.fit <- rbpspline(y, x, knots3, degree = 3, penalty = "cv") ## gcv.fit <- rbpspline(y, x, knots3, degree = 3, penalty = "gcv") ## plot(x, y, col = "lightgrey") ## lines(cv.fit, col = 2) ## lines(gcv.fit, col = 3) ################################################### ### code chunk number 57: SpatialExtremesGuide.Rnw:3115-3116 ################################################### set.seed(15) ################################################### ### code chunk number 58: SpatialExtremesGuide.Rnw:3119-3132 ################################################### n.site <- 20 n.obs <- 50 coord <- matrix(runif(2*n.site, 0, 10), ncol = 2) colnames(coord) <- c("lon", "lat") data <- rmaxstab(n.obs, coord, "gauss", cov11 = 100, cov12 = 25, cov22 = 220) param.loc <- -10 + 2 * coord[,2] param.scale <- 5 + 2 * coord[,1] + coord[,2]^2 param.shape <- rep(0.2, n.site) for (i in 1:n.site) data[,i] <- frech2gev(data[,i], param.loc[i], param.scale[i], param.shape[i]) ################################################### ### code chunk number 59: SpatialExtremesGuide.Rnw:3137-3144 ################################################### loc.form <- y ~ lat scale.form <- y ~ lon + I(lat^2) shape.form <- y ~ 1 shape.form2 <- y ~ lon M1 <- fitspatgev(data, coord, loc.form, scale.form, shape.form) M2 <- fitspatgev(data, coord, loc.form, scale.form, shape.form2) M1 ################################################### ### code chunk number 60: SpatialExtremesGuide.Rnw:3151-3153 ################################################### anova(M1, M2) TIC(M1, M2)