Gaussian Random Fields In Splus

This web page contains Splus functions for generating stationary [not necessarily isotropic] Gaussian random fields over regular grids, and for clipping and visualizing them. There is also a section on computing the Matérn covariance.

The related resources are

• The CGI program illustrating these concepts on-line;

• The source code of the above, including the C library of functions for generation, estimation, and visualization of the Gaussian random fields;

• Kriging in Splus;

• The ISR Technical Report with the underlying theory and algorithms.

Generating Gaussian Random Fields

To generate a realization of the stationary Gaussian random field, you need to provide a covariance function, in the form function (x, y, theta1, theta2, tau). This function depends on (x, y) -- the vector difference between the points where the covariance is computed, and on (theta1, theta2, tau) -- the covariance parameters. For example, consider the Spherical covariance function,

The corresponding Splus code is

sphercov <- function (x, y, theta1, theta2, tau) {         d <- sqrt (x^2 + y^2)         ifelse (d < theta1,              (1- 1.5*(d/theta1)+0.5*(d/theta1)^3) / tau,              0)       }

Note that the covariance function must be able to accept vectors for x and y.

Next, copy-and-paste the two following functions into your Splus session.

embed <- function (n1, n2, covfun, theta1, theta2, tau) {              xf <- rep (c(rep(1, n2), 0, n2:2), n1)              yf <- rep (c(1:n2, 0, rep(1, n2-1)), n1) +                    n2 * 0:(2*n1*n2 - 1) %/% (2*n2)              xb <- rep (c(rep(n1*n2 - n2 + 1, n2), 0,                    (n1*n2):(n1*n2 - n2 + 2)), n1 - 1) -                    n2 * 0:(2*(n1 - 1)*n2 - 1) %/% (2*n2)              yb <- rep (c(1:n2, 0, rep(1, n2 - 1)), n1 - 1)              indexX <- c(xf, rep(0, 2*n2), xb)              indexY <- c(yf, rep(0, 2*n2), yb)              z <- indexX==0 | indexY==0              x1 <- (indexX-1) %/% n2 + 1              y1 <- (indexX-1) %% n2 + 1              x2 <- (indexY-1) %/% n2 + 1              y2 <- (indexY-1) %% n2 + 1              ifelse (z, 0, covfun (x1-x2, y1-y2, theta1, theta2, tau))      }

genField <- function (n1, n2, covfun, theta1, theta2, tau) {         v <- embed(n1, n2, covfun, theta1, theta2, tau)         V <- matrix(v, ncol=2*n2, byrow=T)         LAMBDA <- fft(V, inverse=T)         if (sum(LAMBDA<0)) stop ("\nNegative embedding! \nUse larger grid!")          sqlambda <- sqrt(Re(c(t(LAMBDA))))         epsilon <- rnorm(4*n1*n2, 0, sqlambda) +                 1i*rnorm(4*n1*n2, 0, sqlambda)         EPSILON <- matrix(epsilon, ncol=2*n2, byrow=T)         W <- fft(EPSILON, inverse=T) / sqrt(4*n1*n2)         w <- Re(c(t(W)))         indexZ <- rep(1:n2, n1) + 2*n2 * 0:(n1*n2-1) %/% n2         matrix (w[indexZ], ncol=n2, byrow=T)      }

Once you have the functions cov, embed, and genField, you are ready to generate Gaussian random fields. For example, to generate a 128x128 sample with Spherical(100, 1, 1) correlation, type

z <- genField (128, 128, sphercov, 100, 1, 1)

This will put a 128x128 matrix with the desired field into z. Now you can view it by typing

persp (z)

or

image (z)

A Note On Negative Embedding

The algorithm carried out by the above functions works only if the covariance decreases sufficiently fast, so that the covariance between the points on the opposite sides of the grid is almost zero. Otherwise you get the "Negative embedding" error. For example, it is impossible to generate

z <- genField (13, 13, sphercov, 100, 1, 1)

Instead, you must generate a larger field, and take a portion of it of the desired size. This is OK to do since the field is stationary. Simply type

z <- genField (128, 128, sphercov, 100, 1, 1) z <- z[1:13, 1:13]

Also note that it is faster to generate n x k fields where both n and k are powers of two.

We suggest two functions for writing random fields into PPM images. The first function converts continuous (e.g. Gaussian) data into levels of gray.

writePPM <- function (field, filename){        if (length(dim(field)) != 2) stop ("\nField must be a matrix!")        fmin <- min(c(field))        fmax <- max(c(field))        seg <- (fmax-fmin)/255        greylev <- round((c(field)-fmin) / seg)        cat ("P3", dim(field), "255", sep="\n",             file=filename, append=F)        write (t(cbind(greylev, greylev, greylev)), ncol=15,        file=filename, append=T)     }

For example, if z is a realization of the Gaussian random field, as returned by the genField function, you can type

writePPM (z, "file.ppm")

at Splus prompt to create the file called file.ppm. You can then view it by typing

!xv file.ppm &

The second function is for clipping random fields and coloring the resulting discrete field. To use it, you must first download the file colors.txt and put it into the directory from which you start Splus. Then copy-and-paste the following code:

writeClippedPPM <- function (field, thresholds, colors, filename) {      if (length(dim(field)) != 2)          stop ("\nField must be a matrix!")      n <- length(colors)      if (n != length(thresholds)+1)          stop("\n#thresholds must be \n #colors less one!")      colorDef <- dget("colors.txt")      z <- c(field)      qz <- matrix (rep(colorDef[colors[n],], length(z)),                    ncol=3, byrow=T)      indz <- z < thresholds[1]      len <- sum(indz)      qz[indz,] <- c(rep(colorDef[colors[1], 1], len),      rep(colorDef[colors[1], 2], len),      rep(colorDef[colors[1], 3], len))      if (n>2)          for (i in 1:(n-2)) {               indz <- z < thresholds[i+1] & z > thresholds[i]               len <- sum(indz)               qz[indz,] <- c(rep(colorDef[colors[i+1], 1], len),                              rep(colorDef[colors[i+1], 2], len),                              rep(colorDef[colors[i+1], 3], len))          }      cat ("P3", dim(field), "255", sep="\n", file=filename, append=F)      write (t(qz), ncol=15, file=filename, append=T)    }

Now you can type

writeClippedPPM (z,                  c(-0.43, 0.43),                  c("dodger blue", "white", "grey"),                           "sky.ppm")

This means that you are clipping the random field stored in z at two thresholds, -0.43 and 0.43. Naturally, this results in a three color image. The particular colors to use are specified in the third parameter to the function. Note that the number of colors always exceeds the number of thresholds by one. You can browse the file colors.txt (near the end) to find out what color names are defined. Finally, the resulting image will be stored in the file sky.ppm.

To view the resulting image, type

!xv sky.ppm &

Computing the Matérn Covariance

The Matérn covariance is given by


We show how to evaluate it in Splus.

• Download the file matern.c and put it into the directory from which you start Splus.
• Type at the UNIX prompt (that is, outside Splus):

  Splus COMPILE matern.c

  If all is well, this creates a file matern.o.
• Run Splus and type

  dyn.load ("matern.o")

  This should be done each time you run Splus.
• Copy-and-paste the following function into your Splus session:

  matcov <- function (x, y, theta1, theta2, tau) {       d <- sqrt(x^2 + y^2)       answer <- numeric(length(d))       .C("matern_cov",            as.double(d),            as.integer(length(d)),            as.double(theta1),            as.double(theta2),            as.double(answer))[[5]]/tau   }