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 isotropic Gaussian random fields using
     turning bands, circulant embedding, direct methods, and the random
     coin method.

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

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

     InitGaussRF(x, y=NULL, z=NULL, grid, model, param,
                    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

    grid: logical; determines whether the vectors `x', `y', and `z'
          should be interpreted as a grid definition, see Details.

   model: string; covariance or variogram model, see `CovarianceFct',
          or type `PrintModelList()' to get all options

   param: parameter vector: `param=c(mean, variance, nugget,
          scale,...)'; the parameters must be given in this order;
          further parameters are to be added in case of a parametrised
          class of models, see `CovarianceFct'

  method: `NULL' or string; Method used for simulating, see
          `RFMethods', or type `PrintMethodList()' to get all options

       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' creates an isotropic Gaussian random field with `model'
     as covariance function/variogram model and parameters
     `param=c(mean,variance,nugget,scale,...)'.  The sill of the
     variogram equals `variance + nugget'.

     `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.

     `GaussRF' is split up in an initial `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. 

     Using `InitGaussRF' and `DoSimulateRF' in sequence might be
     slightly faster than `GaussRF' (but less convenient). 

     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))')

        *  `model="nugget"' is one possibility to create independent
           Gaussian random variables.  Without loss of efficiency, any
           covariance function with parameter vector `c(mean, 0,
           nugget, scale, ...)' can also be used.  If `model="nugget"'
           is used the second component of `param' must be zero.

        *  The sum of the components variance and nugget in the
           argument `param' equals the sill of the variogram. 

        *  `register' is a parameter which may never be used by most
           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 occured and a positive
     value if failed.

     `GaussRF' and `DoSimulateRF' return `NULL' if an error has
     occured; 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@uni-bayreuth.de <URL:
     http://www.geo.uni-bayreuth.de/~martin>

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

     Gneiting, T. and Schlather, M. (2001) Statistical modeling with
     covariance functions. In preparation.

     Schlather, M. (1999) An introduction to positive definite
     functions and to unconditional simulation of random fields.
     Technical report ST 99-10, Dept. of Maths and Statistics,
     Lancaster University.

_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
      x <- seq(0, 20, 0.1) 
      y <- seq(0, 20, 0.1)     
      f <- GaussRF(x=x, y=y, model=model, grid=TRUE,
                   param=c(mean, variance, nugget, scale, alpha))
      image(x, y, f)

      #############################################################
      ## ... using gridtriple
      x <- c(0, 20, 0.1)  ## note: vectors of three values, not a 
      y <- c(0, 20, 0.1)  ##       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))
      f

      #############################################################
      ## 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))
      f

      #############################################################
      ## simulating several random fields at once
      x <- seq(0, 20, 0.1)  # grid
      y <- seq(0, 20, 0.1)
      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])


      #############################################################
      ## 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")

