Papers {RandomFields} | R Documentation |
Code is provided that provides the results in the references below.
Martin Schlather, martin.schlather@math.uni-goettingen.de http://www.stochastik.math.uni-goettingen.de/~schlather;
Gneiting, T., Kleiber, W. and Schlather, M. (2011) Matern cross-covariance functions for multivariate random fields J. Amer. Statist. Assoc. 105, 1167-1177.
Gneiting, T., Sevcikova, H., Percival, D.B., Schlather, M., Jiang, Y. (2006) Fast and Exact Simulation of Large Gaussian Lattice Systems in R2: Exploring the Limits. J. Comput. Graph. Stat., 15, 483-501.
Schlather, M. (2002) Models for stationary max-stable random fields. Extremes 5, 33-44.
Scheuerer, M. and Schlather, M. (2011) Covariance Models for Random Vector Fields.
Schlather, M. (2010) On some covariance models based on normal scale mixtures. Bernoulli, 16, 780-797.
############################################################ ## ## ## M. Schlather, Extremes 5:1, 33-44, 2002 ## ## ## ############################################################ col=grey((100:0)/100) pts <- if (interactive()) 512 else 32 x <- (1:pts) / pts * 10 scalegauss <- 1.5 runif(1) root <- 1/2 RFparameters(mpp.plus=2*scalegauss, CE.force=TRUE) save.seed <- .Random.seed ms1 <- MaxStableRF(x, x, model="gauss", param=c(0, 1, 0, scalegauss), maxstable="Bool", grid=TRUE, Print=3) image(x, x, ms1^root, col=col) .Random.seed <- save.seed scalecirc<- 3.3 # 0.16 0.00 0.18 0.00 0.00 ms2 <- MaxStableRF(x, x, model="sph", param=c(0, 1, 0, scalecirc), maxstable="Bool", grid=TRUE, Print=3) image(x, x, ms2^root, col=col) .Random.seed <- save.seed ms3 <- MaxStableRF(x, x, model="exponen", param=c(0, 1, 0, 1), maxstable="extr", grid=TRUE, Print=3, method="ci") image(x, x, ms3^root, col=col) .Random.seed <- save.seed ms4 <- MaxStableRF(x, x, model="gauss", param=c(0, 1, 0, 1), maxstable="extr", grid=TRUE, Print=3, method="ci") image(x, x, ms4^root, col=col) ############################################################ ## ## ## M. Schlather, Bernoulli, 16, 780-797, 2010 ## ## ## ############################################################ # % library(RandomFields, lib="~/TMP") jpg <- jpeg Plot <- function(basicname, x, y, z, pixels=200, col=col, delay=1, basicn=10000, T, redo=TRUE, height=5, zi=1, speed=0.3, zlim = range(z)) { Print(try(dir.create(basicname))) if (length(T)==3) T <- seq(T[1], T[2], T[3]) cex = 2 par(mar=c(4,4.3, 0.8, 0.8), cex=0.7) for (file in list(TRUE, "jpg")) { h <- if (is.logical(file)) height else pixels args <- list(TRUE, file, ps=paste(basicname, "/X2", basicname, sep=""), quiet=FALSE, height=h, width=h) do.call("Dev", args) image(x, T, z[,zi,], zlim=zlim, col=col, cex.axis=cex, cex.lab=cex, ylab="time") Dev(FALSE) args <- list(TRUE, file, ps=paste(basicname, "/X3", basicname, sep=""), quiet=FALSE, height=h, width=h) do.call("Dev", args) image(x, y, z[,,1], zlim=zlim, col=col, cex.axis=cex, cex.lab=cex) Dev(FALSE) } k <- 0 time <- length(T) if (redo) { system(paste("rm ", basicname, "/", basicname, "*.jpg", sep="")) par(mar=rep(0, 4)) Dev(TRUE, "jpg", ps=paste(basicname, "/", basicname, sep=""), quiet=FALSE, height=pixels * height, width=pixels * height) plot(Inf, Inf, axes=FALSE, xlim=c(0,1), ylim=c(0,1)) Dev(FALSE) for (i in 1:(time)) { for (j in 1:delay) { k <- k + 1 name <- paste(basicname, "/", basicname, k + basicn, sep="") Dev(TRUE, "jpg", ps=name, quiet=FALSE, height=pixels, width=pixels) image(x, y, z[,,i], zlim=zlim, col=col, axes=FALSE) Dev(FALSE) } } } system(paste("cd ", basicname, ";mencoder -mf fps=30 -ffourcc DX50 -ovc lavc ", " -speed ", speed, " -lavcopts vcodec=mpeg4:vbitrate=9800:aspect=4/3:vhq:keyint=15", " -vf scale=720:576 -o ", basicname, ".avi mf://", basicname, "*.jpg &", sep="")) } ### Example 10 in Schlather (2010) basicname <- "coxisham" y <- x <- seq(0, 10, len=if (interactive()) 256 else 5) T <- c(0, 25.5, if (interactive()) 0.1 else 5) col <- c(topo.colors(300)[1:100], cm.colors(300)[c((1:50) * 2, 101:150)]) model <- list("coxisham", mu=c(1, 1), D=matrix(nr=2, c(1, 0.5, 0.5, 1)), list("whittle", nu=1) ) z <- GaussRF(x, y, T=T, grid =TRUE, spectral.lines=1500, every=10, model = model, Print=2) zlim <- range(z) time <- dim(z)[3] for (i in 1:time) { cat(i,"\n") image(x, y, z[, , i], add=i>1, col=col, zlim=zlim) } # load(paste(basicname, ".dat", sep="")) Plot(basicname, x, y, z, col=col, T=T) save(file=paste(basicname, "/", basicname, ".dat", sep=""), z) ### Example 13 in Schlather (2010) basicname <- "moving" y <- x <- seq(0, 10,len = if (interactive()) 256 else 10) T <- c(0, 25.5, if (interactive()) 0.1 else 10) col <- c(topo.colors(300)[1:100], cm.colors(300)[c((1:50) * 2, 101:150)]) model <- list("ave1", A = matrix(nc=2,c(0.5, 0, 0, 1)), z= c(2,0), list("whittle", nu=1) ) z <- GaussRF(x, x, T=T, model=model, grid=TRUE, Print=2, me="av", every = 10, CE.trial=2, CE.force=TRUE, CE.maxmem=16777216*8) zlim <- range(z) time <- dim(z)[3] for (i in 1:time) { Print(i) image(x, y, z[, , i], add=i>1, col=col, zlim=zlim) } # load(paste(basicname, ".dat", sep="")) Plot(basicname, x, y, z, col=col, T=T) save(file=paste(basicname, "/", basicname, ".dat", sep=""), z) ### Example 16 in Schlather (2010) intens <- if (interactive()) 100 else 3 ## 1000 takes a huge amount ## ## of time; take smaller values basicname <- "cyclone" len <- if (interactive()) 81 else 10 y <- x <- seq(-3, 3, len=len) T <- seq(0, 6, len=len) col <- c(topo.colors(300)[1:100], cm.colors(300)[c((1:50) * 2, 101:150)]) model <- list("stp", M=matrix(nc=3, rep(0, 9)), Sx=list("EtAxxA", E=c(1, 1, 1), alpha = -2 * pi, A=t(matrix(nc=3, c(2, 0, 0, 1, 1 , 0, 0, 0, 0)))), H = list("Rotat", phi= -2 * pi), phi = list("nonstWM", nu = 1) ) z <- GaussRF(x, x, z=T, model=model, grid=TRUE, me="coin", every = 10, mpp.intens=intens, mpp.p = 0.1, mpp.beta=3.5, mpp.plus = 8, Print=3) zlim <- c(-3.5, 3.5) time <- dim(z)[3] for (i in 1:time) {Print(i);image(x, y, z[,,i], add=i>1, col=col, zlim=zlim)} # load(paste(basicname, ".dat", sep="")) Plot(basicname, x, y, z, col=col, T=T, pixels=256, zi=0.5 + dim(z)[2]/2, speed=0.2, zlim=zlim) save(file=paste(basicname, "/", basicname, ".dat", sep=""), z) ############################################################ ## ## ## T. Gneiting, W. Kleiber, M. Schlather, JASA 2011 ## ## ## ############################################################ ## see help("weather") ############################################################ ## ## ## Figure 1 in ## ## Gneiting, T., Sevcikova, H., ## ## Percival, D.B., Schlather, M. and Jiang, Y. (2005) ## ## ## ############################################################ stabletest <- function(alpha, theta, size=if (interactive()) 512 else 20) { RFparameters(CE.trials=1, CE.tolIm = 1e-8, CE.tolRe=0, CE.force = FALSE, CE.useprimes=TRUE, CE.strategy=0, skipchecks=!FALSE, Print=2, Storing=TRUE) model <- list("cutoff", diameter=theta, a=1, list(model="stable", alpha=alpha)) CovarianceFct(0, model=model, dim=2) r <- GetModelInfo(-1)$sub[[1]]$internalq[5] # theor R x <- c(0, r, r / (size - 1)) * theta err <- try(InitGaussRF(x, x, grid=TRUE, model=model, meth="circulant")) return(if (!is.numeric(err) || err != 0) NA else r) } alphas <- seq(1.52, 2.0, if (interactive()) 0.02 else 0.1) thetas <- if (interactive()) seq(0.05, 3.5, 0.05) else seq(0.1, 3.5, 0.2) m <- matrix(NA, nrow=length(thetas), ncol=length(alphas)) for (it in 1:length(thetas)) { theta <- thetas[it] for (ia in 1:length(alphas)) { alpha <- alphas[ia] cat("alpha=", alpha, "theta=", theta,"\n") print(unix.time(m[it, ia] <- stabletest(alpha=alpha, theta=theta))) if (is.na(m[it, ia])) break } if (any(is.finite(m))) image(thetas, alphas, m, col=rainbow(100)) } ############################################################ ## ## ## M. Scheuerer, M. Schlather, 2011 ## ## ## ############################################################ # % library(RandomFields, lib="~/TMP") my.legend <- function(lu.x, lu.y, zlim, col, cex=1) { ## uses already the legend code of R-1.3.0 cn <- length(col) filler <- vector("character", length=(cn-3)/2) legend(lu.x, lu.y, y.i=0.03, x.i=0.1, legend=c(format(zlim[2], dig=2), filler, format(mean(zlim), dig=2), filler, format(zlim[1], dig=2)), lty=1, col=rev(col),cex=cex) } my.arrows <- function(xy, z, r, thinning) { startx <- as.vector(xy[,1] - r/2*z[as.integer(dim(z)[1]/3) + 1,,]) starty <- as.vector(xy[,2] - r/2*z[as.integer(dim(z)[1]/3) + 2,,]) endx <- as.vector(xy[,1] + r/2*z[as.integer(dim(z)[1]/3) + 1,,]) endy <- as.vector(xy[,2] + r/2*z[as.integer(dim(z)[1]/3) + 2,,]) startx[c(rep(TRUE, thinning), FALSE)] <- NA starty[c(rep(TRUE, thinning), FALSE)] <- NA endx[c(rep(TRUE, thinning), FALSE)] <- NA endy[c(rep(TRUE, thinning), FALSE)] <- NA arrows(x0=startx, y0=starty, x1=endx, y1=endy, length=0.03) } x <- c(-3, 3, if (interactive()) 0.049 else 0.5) nu <- 3 col <- grey(seq(0, 1, 0.01)) thinning <- 21 length.arrow <- 1.5 / thinning runif(1) seed <- .Random.seed eps <- interactive() # true falls eps/pdf drucken for (modelname in c("divfree", "curlfree")) { cat(modelname, "\n") model <- list(modelname, list("matern", nu=nu)) xx <- seq(x[1], x[2], x[3]) RFparameters(Print=2) if (!eps) { cf <- Covariance(x=cbind(xx,0), model=model) do.call(getOption("device"), list(height=5, width=5)) par(mfcol=c(1,1)) j <- 3 plot(xx, cf[(j-1) * length(xx) + (1 : length(xx)), j]) } .Random.seed <- seed z <- GaussRF(x, x, grid=TRUE, gridtr=TRUE, model=model, n=1, CE.trial=2, Stor=TRUE, me="ci", Print=3, CE.force=!TRUE) ## z[1,,] : Potentialfeld ## z[2:3,,] : vectorfeld ## z[4,,] : div bzw. rot if (eps) { do.call(getOption("device"), list(height=5, width=5)) par(mfcol=c(1,1)) } else { do.call(getOption("device"), list(height=5, width=10)) par(mfcol=c(1,2)) } par(mar=c(2.2,2.5,0.5,0.5), cex.axis=2, bg="white") for (no.vectors in c(TRUE,FALSE)) for (i in c(1,2,4)) { ## image i=1: Potential feld + Vektorfeld ## image i=4: div/rot feld + Vektorfeld if (i==1 || i==4) { image(xx, xx, col=col, z[i,,] ) if (no.vectors) my.legend(max(xx) - 0.3 * diff(range(xx)), min(xx) + 0.3 * diff(range(xx)), zlim=range(z[i,,]), col=col, cex=1.5) } else plot(Inf, Inf, xlim=range(xx), ylim=range(xx)) if (!no.vectors || i!=1 && i!=4) { xy <- as.matrix(expand.grid(xx, xx)) my.arrows(xy, z, length.arrow, thinning) } if (all(par()$mfcol==1)) { name <- paste(modelname, "_", nu, "_", !no.vectors, "_", i, sep="") cat(name,"\n") dev.copy2eps(file=paste(name, ".eps", sep="")) dev.copy2pdf(file=paste(name, ".pdf", sep="")) } } }