GaussRF             package:RandomFields             R Documentation

_G_a_u_s_s_i_a_n _R_a_n_d_o_m _F_i_e_l_d_s

_D_e_s_c_r_i_p_t_i_o_n:

     These functions simulate stationary spatial and spatio-temporal
     Gaussian random fields using turning bands/layers, circulant
     embedding, direct methods, and the random coin method.

_U_s_a_g_e:

     GaussRF(x, y=NULL, z=NULL, T=NULL, grid, model, param, trend,
             method=NULL, n=1, register=0, gridtriple=FALSE)

     InitGaussRF(x, y=NULL, z=NULL, T=NULL, grid, model, param, trend,
                 method=NULL, register=0, gridtriple=FALSE)

_A_r_g_u_m_e_n_t_s:

       x: matrix of coordinates, or vector of x coordinates

       y: vector of y coordinates

       z: vector of z coordinates

       T: vector of time coordinates, may only be given if the random
          field is defined as an anisotropic random field, i.e. if
          'model=list(list(model=,var=,k=,aniso=),...)'. 'T' must
          always be given in the 'gridtriple' format, independently how
          the spatial part is defined. 

    grid: logical; determines whether the vectors 'x', 'y', and 'z'
          should be interpreted as a grid definition, see Details. 
          'grid' does not apply for 'T'.

   model: string or list; covariance or variogram model, see
          'CovarianceFct', or type 'PrintModelList()' to get the list
          of all implemented models; see Details.

   param: vector or matrix of parameters or missing, see Details and
          'CovarianceFct';  The simplest form is that 'param' is vector
          of the form 'param=c(NA,variance,nugget,scale,...)', in this
          order;
           The dots '...' stand for additional parameters of the model.

   trend: Not programmed yet. trend surface: number (mean) or a vector
          of length d+1 (linear trend a_0 +a_1 x_1 + ... + a_d x_d), or
          function(x)

  method: 'NULL' or string; method used for simulating, see
          'RFMethods', or type 'PrintMethodList()' to get all options.
          If 'model' is given as list then 'method' may not be set if
          'model[[i]]$method', i=1,3,.. is given, and vice versa.

       n: number of realisations to generate

register: 0:9; place where intermediate calculations are stored; the
          numbers are aliases for 10 internal registers

gridtriple: logical. Only relevant if 'grid==TRUE'. If
          'gridtriple==TRUE' then 'x', 'y', and 'z' are of the form
          'c(start,end,step)'; if 'gridtriple==FALSE' then 'x', 'y',
          and 'z' must be vectors of ascending values 

_D_e_t_a_i_l_s:

     'GaussRF' can use different methods for the simulation, i.e.,
     circulant embedding, turning bands, direct methods, and random
     coin method.  If 'method==NULL' then 'GaussRF' searches for a
     valid method.  'GaussRF' may not find the fastest method neither
     the most precise one.  It just finds any method among the
     available methods. (However it guesses what is a good choice.) 
     Note that some of the methods do not work for all covariance  or
     variogram models.

        *  An isotropic random field is created by  'GaussRF' where
           'model' is the covariance or variogram model and the
           parameter is 'param=c(mean,variance,nugget,scale, ...)'.
           Alternatively the 'trend' can be given (not programmed yet);
           then 'param=c(variance,nugget,scale, ...)'.

        *  Nested models can be defined in the same way as a nested
           'CovarianceFct'.  If the 'trend' is not given it is set to
           0.

        *  An anisotropic random field (i.e. zonal anisotropy,
           geometrical anisotropy, separable models, non-separable
           space-time models) and a random field based on
           multiplicative or nested models is defined as in the case of
           an anisotropic 'CovarianceFct'.  If the 'trend' is not given
           it is set to 0. The 'method' may be specified by the global
           'method' or for each model separately, as additional
           parameter 'method' for each entry of the list; note that
           methods can not be mixed within a multiplicative part.

           If 'model=list(list(model=,var=,k=,aniso=),...)' then a time
           component might be given.  In case of 'model="nugget"',
           'aniso' must still be given as a matrix.  Namely if 'aniso'
           is a singular matrix then a zonal nugget effect is obtained.

     'GaussRF' calls initially 'InitGaussRF', which does some basic
     checks on the validity of the parameters.  Then, 'InitGaussRF'
     performs some first calculations, like the first Fourier transform
     in the circulant embedding method or the matrix decomposition for
     the direct methods.  Random numbers are not involved.  'GaussRF'
     then calls 'DoSimulateRF' which uses the intermediate results and
     random numbers to create a simulation.

     When 'InitGaussRF' checks the validity of the parameters, it also
     checks whether the previous simulation has had the same
     specification of the random field.  If so (and if
     'RFparameters()$STORING==TRUE'), the stored intermediate results
     are used instead of being recalculated. 

     Comments on specific parameters:

        *  'grid==FALSE' : the vectors 'x', 'y', and 'z' are
           interpreted as vectors of coordinates

        *  '(grid==TRUE) && (gridtriple==FALSE)' : the vectors 'x',
           'y', and 'z' are increasing sequences with identical lags
           for each sequence.  A corresponding grid is created (as
           given by 'expand.grid'). 

        *  '(grid==TRUE) && (gridtriple==FALSE)' : the vectors 'x',
           'y', and 'z' are triples of the form (start,end,step)
           defining a grid (as given by
           'expand.grid(seq(x$start,x$end,x$step),
           seq(y$start,y$end,y$step), seq(z$start,z$end,z$step))')

        *  'register' is a parameter which may never be used by most of
           the users (please let me know if you use it!).  In other
           words, the package will work fine if you ignore this
           parameter.  The parameter 'register' is of interest in the
           following situation.  Assume you wish to create sequentially
           several realisations of two random fields Z1 and Z2 that
           have different specifications of the covariance/variogram
           models, i.e. Z1, Z2, Z1, Z2,... Then, without using
           different registers, the algorithm will not be able to
           profit from already calculated intermediate results, as the
           specifications of the covariance/variogram model change
           every time.  However, using different registers allows for
           profiting from up to 10 stored intermediate results. 

        *  The strings for 'model' and 'method' may be abbreviated as
           long as the abbreviations match only one option.  See also
           'PrintModelList()' and 'PrintMethodList()'

        *  Further control parameters for the simulation are set by
           means of 'RFparameters(...)'.

_V_a_l_u_e:

     'InitGaussRF' returns 0 if no error has occurred and a positive
     value if failed.

     'GaussRF' and 'DoSimulateRF' return 'NULL' if an error has
     occurred; otherwise the returned object depends on the parameters
     'n' and 'grid':
      'n==1':
      * 'grid==FALSE'.  A vector of simulated values is returned
     (independent of the dimension of the random field)
      * 'grid==TRUE'.  An array of the dimension of the random field is
     returned.

     'n>1':
      * 'grid==FALSE'.  A matrix is returned.  The columns contain the
     repetitions.
      * 'grid==TRUE'.  An array of dimension d+1, where d is the
     dimension of the random field, is returned.  The last dimension
     contains the repetitions.

_N_o_t_e:

     The algorithms for all the simulation methods are controlled by
     additional parameters, see 'RFparameters()'.  These parameters
     have an influence on the speed of the algorithm and the precision
     of the result.  The default parameters are chosen such that the
     simulations are fine for many models and their parameters.  If in
     doubt modify the example in 'EmpiricalVariogram()' to check the
     precision.

_A_u_t_h_o_r(_s):

     Martin Schlather, martin.schlather@cu.lu <URL:
     http://www.cu.lu/~schlathe>

_R_e_f_e_r_e_n_c_e_s:

     See RFMethods for the references.

_S_e_e _A_l_s_o:

     'CovarianceFct', 'DeleteRegister', 'DoSimulateRF',
     'GetPracticalRange', 'EmpiricalVariogram', 'mleRF', 'MaxStableRF',
     'RFMethods', 'RandomFields', 'RFparameters', 'ShowModels',

_E_x_a_m_p_l_e_s:

      #############################################################
      ##                                                         ##
      ## Examples using the symmetric stable model, also called  ##
      ## "powered exponential model"                             ## 
      ##                                                         ##
      #############################################################
      PrintModelList()    ## the complete list of implemented models
      model <- "stable"   
      mean <- 0
      variance <- 4
      nugget <- 1
      scale <- 10
      alpha <- 1   ## see help("CovarianceFct") for additional
                   ## parameters of the covariance functions
      step <- 1    ## nicer, but also time consuming if step <- 0.1
      x <- seq(0, 20, step) 
      y <- seq(0, 20, step)     
      f <- GaussRF(x=x, y=y, model=model, grid=TRUE,
                   param=c(mean, variance, nugget, scale, alpha))
      image(x, y, f)

      #############################################################
      ## ... using gridtriple
      step <- 1    ## nicer, but also time consuming if step <- 0.1
      x <- c(0, 20, step)  ## note: vectors of three values, not a 
      y <- c(0, 20, step)  ##       sequence
      f <- GaussRF(grid=TRUE, gridtriple=TRUE,
                    x=x ,y=y, model=model,  
                    param=c(mean, variance, nugget, scale, alpha))
      image(seq(x[1],x[2],x[3]), seq(y[1],y[2],y[3]), f)

      #############################################################
      ## arbitrary points
      x <- runif(100, max=20) 
      y <- runif(100, max=20)
      z <- runif(100, max=20) # 100 points in 3 dimensional space
     (f <- GaussRF(grid=FALSE,
                   x=x, y=y, z=z, model=model, 
                   param=c(mean, variance, nugget, scale, alpha)))

      #############################################################
      ## usage of a specific method
      ## -- the complete list can be obtained by PrintMethodList()
      x <- runif(100, max=20) # arbitrary points
      y <- runif(100, max=20)
      (f <- GaussRF(method="dir",  # direct matrix decomposition
                   x=x, y=y, model=model, grid=FALSE, 
                   param=c(mean, variance, nugget, scale, alpha)))

      #############################################################
      ## simulating several random fields at once
      step <- 1    ## nicer, but also time consuming if step <- 0.1
      x <- seq(0, 20, step)  # grid
      y <- seq(0, 20, step)
      f <- GaussRF(n=3,  # three simulations at once
                   x=x, y=y, model=model, grid=TRUE,  
                   param=c(mean, variance, nugget, scale, alpha))
      image(x, y, f[,,1])
      image(x, y, f[,,2])
      image(x, y, f[,,3])
             

      #############################################################
      ##                                                         ##
      ##      Examples using the extended definition forms       ##
      ##                                                         ##
      ##                                                         ##  
      #############################################################

     ## note that the output seems plausible but not checked!!!!

     ## tbm may also be used for multiplicate models (if they have
     ## *exactly* the same anisotropy parameters)
     x <- (0:100)/10
     m <- matrix(c(1,2,3,4),ncol=2)/5
     z <- GaussRF(x=x, y=x, grid=TRUE,
                   model=list(
                     list(m="power",v=1,k=2,a=m),
                     "*", list(m="sph", v=1, a=m)
                     ),
                   me="TBM3", reg=0,n=1)
     print(c(mean(as.double(z)),var(as.double(z))))
     image(z,zlim=c(-3,3))

     ## non-separable space-time model applied for two space dimensions
     ## note that tbm method does not work nicely, but at least
     ## in some special cases.
     x <- y <- (1:32)/2     ## grid definition, but as a sequence
     T <- c(1,32,1)*10      ## note necessarily gridtriple definition
     aniso <- diag(c(0.5,8,1)) 
     k <- c(1,phi=1,1,1,psi=1,dim=2)
     model <- list(list(m="nsst", v=1, k=k, a=aniso))
     z <- GaussRF(x=x, y=y, T=T, grid=TRUE, model=model)
     rl <- function() if (interactive()) readline("Press return")
     for (i in 1:dim(z)[3]) { image(z[,,i]); rl();}
     for (i in 1:dim(z)[2]) { image(z[,i,]); rl();}
     for (i in 1:dim(z)[1]) { image(z[i,,]); rl();}


      #############################################################
      ##                                                         ##
      ##                    Brownian motion                      ##
      ##                  using Stein's method                   ##
      ##                                                         ##  
      #############################################################
      # 2d
      step <- 0.3  ## nicer, but also time consuming if step <- 0.1
      x <- seq(0, 10, step)
      kappa <- 1   # in [0,2)
      z <- GaussRF(x=x, y=x, grid=TRUE, model="2d", param=c(0,1,0,1,kappa))
      image(z,zlim=c(-3,3))

      # 3d
      x <- seq(0, 3, step)
      kappa <- 1   # in [0,2)
      z <- GaussRF(x=x, y=x, z=x, grid=TRUE, model="3d",
                   param=c(0,1,0,1,kappa))
     rl <- function() if (interactive()) readline("Press return")
     for (i in 1:dim(z)[1]) { image(z[i,,]); rl();}


      #############################################################
      ## This example shows the benefits from stored,            ##
      ## intermediate results: in case of the circulant          ##
      ## embedding method, the speed is doubled in the second    ##
      ## simulation.                                             ##  
      #############################################################
      DeleteAllRegisters()
      RFparameters(Storing=TRUE,PrintLevel=1)
      y <- x <- seq(0, 50, 0.2)
      (p <- c(runif(3), runif(1)+1))
      ut <- unix.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
                                   method="circ", param=p))
      image(x, y, f)
      hist(f)
      c( mean(as.vector(f)), var(as.vector(f)) )
      cat("unix time (first call)", format(ut,dig=3),"\n")

      # second call with the *same* parameters is much faster:
      ut <- unix.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen",
                                   method="circ",param=p)) 
      image(x, y, f)
      hist(f)
      c( mean(as.vector(f)), var(as.vector(f)) )
      cat("unix time (second call)", format(ut,dig=3),"\n")

