Package 'RandomFields'

Description

The package RandomFields offers various tools for

1. model estimation (ML) and inference (tests) for regionalized variables and data analysis,

2. simulation of different kinds of random fields, including

   • multivariate, spatial, spatio-temporal, and non-stationary Gaussian random fields,

   • Poisson fields, binary fields, Chi2 fields, t fields and

   • max-stable fields.

   It can also deal with non-stationarity and anisotropy of these processes and conditional simulation (for Gaussian random fields, currently).

See https://www.wim.uni-mannheim.de/schlather/publications/software for intermediate updates.

Details

The following features are provided by the package:

1. Bayesian Modelling

   • See Bayesian Modelling for an introduction to hierarchical modelling.

2. Coordinate systems

   • Cartesian, earth and spherical coordinates are recognized, see coordinate systems for details.

   • A list of valid models is given by spherical models.

3. Data and example studies: Some data sets and published code are provided to illustrate the syntax and structure of the package functions.

   • soil : soil physical data

   • weather : UWME weather data

   • papers : code used in the papers published by the author(s)

4. Estimation of parameters (for second-order random fields)

   • RFfit : general function for estimating parameters; (for Gaussian random fields)

   • RFhurst : estimation of the Hurst parameter

   • RFfractaldim : estimation of the fractal dimension

   • RFvariogram : calculates the empirical variogram

   • RFcov : calculates the empirical (auto-)covariance function

5. Graphics

   • Fitting a covariance function manually RFgui

   • the generic function plot

   • global graphical parameters with RFpar

6. Inference (for Gaussian random fields)

   • RFcrossvalidate : cross validation

   • RFlikelihood : likelihood

   • RFratiotest : likelihood ratio test

   • AIC, AICc, BIC, anova, logLik

7. Models

   • For an introduction and general properties, see RMmodels.

   • For an overview over classes of covariance and variogram models –e.g. for geostatistical purposes– see RM. More sophisticated models and covariance function operators are included.

   • RFformula reports a new style of passing a model since version 3.3.

   • definite models are evaluated by RFcov, RFvariogram and RFcovmatrix. For a quick impression use plot(model).

   • non-definite models are evaluated by RFfctn and RFcalc

   • RFlinearpart returns the linear part of a model

   • RFboxcox deals explicitly with Box-Cox transformations. In many cases it is performed implicitly.

8. Prediction (for second-order random fields)

   • RFinterpolate : kriging, including imputing

9. Simulation

   • RFsimulate: Simulation of random fields, including conditional simulation. For a list of all covariance functions and variogram models see RM. Use plot for visualisation of the result.

10. S3 and S4 objects

    • The functions return S4 objects based on the package sp, if spConform=TRUE. This is the default.

      If spConform=FALSE, simple objects as in version 2 are returned. These simple objects are frequently provided with an S3 class. This option makes the returning procedure much faster, but currently does not allow for the comfortable use of plot.

    • plot, print, summary, sometimes also str recognise these S3 and S4 objects

    • use sp2RF for an explicit transformation of sp objects to S4 objects of RandomFields.

    • Further generic functions are available for fitted models, see ‘Inference’ above.

11. Xtended features, especially for package programmers

    • might decide on a large variety of arguments of the simulation and estimation procedures using the function RFoptions

    • may use ‘./configure –with-tcl-config=/usr/lib/tcl8.5/tclConfig.sh –with-tk-config=/usr/lib/tk8.5/tkConfig.sh’ to configure R

Changings

A list of major changings from Version 2 to Version 3 can be found in MajorRevisions.

Changings lists some further changings, in particular of argument and argument names.

RandomFields should be rather stable when running it with parallel. However RandomFields might crash severely if an error occurs when running in parallel. When used with parallel, you might set RFoptions(cores = 1). Note that RFoptions(cores = ...) with more than 1 core uses another level of parallelism which will be in competetions with parallel during runtime.

Update

Current updates are available through https://www.wim.uni-mannheim.de/schlather/publications/software.

Contributions

• Contributions to version 3.0 and following:
  Felix Ballani (TU Bergakademie Freiberg; Poisson Polygons, 2014)
  Daphne Boecker (Univ. Goettingen; RFgui, 2011)
  Katharina Burmeister (Univ. Goettingen; testing, 2012)
  Sebastian Engelke (Univ. Goettingen; RFvariogram, 2011-12)
  Sebastian Gross (Univ. Goettingen; tilde formulae, 2011)
  Alexander Malinowski (Univ. Mannheim; S3, S4 classes 2011-13)
  Juliane Manitz (Univ. Goettingen; testing, 2012)
  Johannes Martini (Univ. Goettingen; RFvariogram, 2011-12)
  Ulrike Ober (Univ. Goettingen; help pages, testing, 2011-12)
  Marco Oesting (Univ. Mannheim; Brown-Resnick processes, Kriging, Trend, 2011-13)
  Paulo Ribeiro (Unversidade Federal do Parana; code adopted from geoR, 2014)
  Kirstin Strokorb (Univ. Mannheim; help pages, 2011-13)
  

• Contributions to version 2.0 and following:
  Peter Menck (Univ. Goettingen; multivariate circulant embedding)
  R Core Team, Richard Singleton (fft.c and advice)

• Contributions to version 1 and following:
  Ben Pfaff, 12167 Airport Rd, DeWitt MI 48820, USA making available an algorithm for AVL trees (avltr*)

Thanks

Patrick Brown : comments on Version 3
Paulo Ribeiro : comments on Version 1
Martin Maechler : advice for Version 1

Financial support

• V3.0 has been financially supported by the German Science Foundation (DFG) through the Research Training Group 1953 ‘Statistical Modeling of Complex Systems and Processes — Advanced Nonparametric Approaches’ (2013-2018).

• V3.0 has been financially supported by Volkswagen Stiftung within the project ‘WEX-MOP’ (2011-2014).

• Alpha versions for V3.0 have been financially supported by the German Science Foundation (DFG) through the Research Training Groups 1644 ‘Scaling problems in Statistics’ and 1023 ‘Identification in Mathematical Models’ (2008-13).

• V1.0 has been financially supported by the German Federal Ministry of Research and Technology (BMFT) grant PT BEO 51-0339476C during 2000-03.

• V1.0 has been financially supported by the EU TMR network ERB-FMRX-CT96-0095 on “Computational and statistical methods for the analysis of spatial data” in 1999.

Note

The following packages enable further choices for the optimizer (instead of optim) in RandomFields: optimx, soma, GenSA, minqa, pso, DEoptim, nloptr, RColorBrewer, colorspace

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Singleton, R.C. (1979). In Programs for Digital Signal Processing Ed.: Digital Signal Processing Committee and IEEE Acoustics, Speech, and Signal Processing Committe (1979) IEEE press.

• Schlather, M., Malinowski, A., Menck, P.J., Oesting, M. and Strokorb, K. (2015) Analysis, simulation and prediction of multivariate random fields with package RandomFields. Journal of Statistical Software, 63 (8), 1-25, url = ‘http://www.jstatsoft.org/v63/i08/’

• see also the corresponding vignette.

See Also

See also RF, RM, RP, RR, RC, R.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

# simulate some data first (Gaussian random field with exponential
# covariance; 6 realisations)
model <- RMexp()
x <- seq(0, 10, 0.1)
z <- RFsimulate(model, x, x, n=6)

## select some data from the simulated data
xy <- coordinates(z)
pts <- sample(nrow(xy), min(100, nrow(xy) / 2))
dta <- matrix(nrow=nrow(xy), as.vector(z))[pts, ]
dta <- cbind(xy[pts, ], dta)
plot(z, dta)


## re-estimate the parameter (true values are 1)
estmodel <- RMexp(var=NA, scale=NA)
(fit <- RFfit(estmodel, data=dta))

## show a kriged field based on the estimated parameters
kriged <- RFinterpolate(fit, x, x, data=dta)
plot(kriged, dta)

Description

These models define particular ways to simulate Brown-Resnick processes.

Usage

RPbrmixed(phi, tcf, xi, mu, s, meshsize, vertnumber, optim_mixed,
          optim_mixed_tol,lambda, areamat, variobound) 

RPbrorig(phi, tcf, xi, mu, s)

RPbrshifted(phi, tcf, xi, mu, s)

RPloggaussnormed(variogram, prob, optimize_p, nth, burn.in, rejection)

Arguments

Details

The argument xi is always a number, i.e. $\xi$ is constant in space. In contrast, $\mu$ and $s$ might be constant numerical values or given an RMmodel, in particular by an RMtrend model.

The functions RPbrorig, RPbrshifted and RPbrmixed simulate a Brown-Resnick process, which is defined by

$Z(x) = \max_{i=1}^\infty X_i \exp(W_i(x) - \gamma),$

where the $X_i$ are the points of a Poisson point process on the positive real half-axis with intensity $x^{-2} dx$, $W_i \sim W$ are iid centered Gaussian processes with stationary increments and variogram $\gamma$ given by model. The functions correspond to the following ways of simulation:

RPbrorig

simulation using the original definition (method 0 in Oesting et al., 2012)

RPbrshifted

simulation using a random shift (similar to method 1 and 2)

RPbrmixed

simulation using M3 representation (method 4)

Value

The functions return an object of class RMmodel.

Note

Advanced options for RPbroriginal and RPbrshifted are maxpoints and max_gauss, see RFoptions.

Author(s)

Marco Oesting, [email protected], https://www.isa.uni-stuttgart.de/institut/team/Oesting/; Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Oesting, M., Kabluchko, Z. and Schlather M. (2012) Simulation of Brown-Resnick Processes, Extremes, 15, 89-107.

See Also

RPbrownresnick, RMmodel, RPgauss, maxstable, maxstableAdvanced.

Examples

#
RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

## currently does not work

Description

RPbrownresnick defines a Brown-Resnick process.

Usage

RPbrownresnick(phi, tcf, xi, mu, s)

Arguments

Details

The argument xi is always a number, i.e. $\xi$ is constant in space. In contrast, $\mu$ and $s$ might be constant numerical values or (in future!) be given by an RMmodel, in particular by an RMtrend model.
For $xi=0$, the default values of $mu$ and $s$ are $0$ and $1$, respectively. For $xi\not=0$, the default values of $mu$ and $s$ are $1$ and $|\xi|$, respectively, so that it defaults to the standard Frechet case if $\xi > 0$.

The functions RPbrorig, RPbrshifted and RPbrmixed perform the simulation of a Brown-Resnick process, which is defined by

$Z(x) = \max_{i=1}^\infty X_i \exp(W_i(x) - \gamma^2),$

where the $X_i$ are the points of a Poisson point process on the positive real half-axis with intensity $x^{-2} dx$, $W_i \sim W$ are iid centered Gaussian processes with stationary increments and variogram $\gamma$ given by phi.

For simulation, internally, one of the methods RPbrorig, RPbrshifted and RPbrmixed is chosen automatically.

Note

Advanced options are maxpoints and max_gauss, see RFoptions.

Further advanced options related to the simulation methods RPbrorig, RPbrshifted and RPbrmixed can be found in the paragraph ‘Specific method options for Brown-Resnick Fields’ in RFoptions.

Author(s)

Marco Oesting, [email protected], https://www.isa.uni-stuttgart.de/institut/team/Oesting/; Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Brown, B.M. and Resnick, S.I. (1977). Extreme values of independent stochastic processes. J. Appl. Probab. 14, 732-739.

• Buishand, T., de Haan , L. and Zhou, C. (2008). On spatial extremes: With application to a rainfall problem. Ann. Appl. Stat. 2, 624-642.

• Kabluchko, Z., Schlather, M. and de Haan, L (2009) Stationary max-stable random fields associated to negative definite functions Ann. Probab. 37, 2042-2065.

• Oesting, M., Kabluchko, Z. and Schlather M. (2012) Simulation of Brown-Resnick Processes, Extremes, 15, 89-107.

See Also

RPbrorig, RPbrshifted, RPbrmixed, RMmodel, RPgauss, maxstable, maxstableAdvanced.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again


## for some more sophisticated models see 'maxstableAdvanced'

Description

This data set contains the calcium content measured in soil samples taken from the 0-20cm layer at 178 locations within a certain study area divided in three sub-areas. The elevation at each location was also recorded.

The first region is typically flooded during the rain season and not used as an experimental area. The calcium levels would represent the natural content in the region. The second region has received fertilisers a while ago and is typically occupied by rice fields. The third region has received fertilisers recently and is frequently used as an experimental area.

Usage

data(ca20)

Format

The object ca20 belongs to the class geodata and is a list with the following elements:

coords

a matrix with the coordinates of the soil samples.

data

calcium content measured in $mmol_c/dm^3$.

covariate

a data-frame with the covariates

altitude

a vector with the elevation of each sampling location, in meters ($m$).

area

a factor indicating the sub area to which the locations belongs.

borders

a matrix with the coordinates defining the borders of the area.

reg1

a matrix with the coordinates of the limits of the sub-area 1.

reg1

a matrix with the coordinates of the limits of the sub-area 2.

reg1

a matrix with the coordinates of the limits of the sub-area 3.

Source

The data was collected by researchers from PESAGRO and EMBRAPA-Solos, Rio de Janeiro, Brasil and provided by Dra. Maria Cristina Neves de Oliveira.

Capeche, C.L.; Macedo, J.R.; Manzatto, H.R.H.; Silva, E.F. (1997) Caracterização pedológica da fazenda Angra - PESAGRO/RIO - Estação experimental de Campos (RJ). (compact disc). In: Congresso BRASILEIRO de Ciência do Solo. 26., Informação, globalização, uso do solo; Rio de Janeiro, 1997. trabalhos. Rio de Janeiro: Embrapa/SBCS.

References

Oliveira, M.C.N. (2003) Métodos de estimação de parâmetros em modelos geoestatísticos com diferentes estruturas de covariâncias: uma aplicação ao teor de cálcio no solo. Tese de Doutorado, ESALQ/USP/Brasil.

Further information on the package geoR can be found at:
http://www.leg.ufpr.br/geoR/.

Description

• Version 3.3

  • RFempiricalvariogram, RFempiricalcovariance and RFempiricalmadogram became obsolete. Use RFvariogram, RFcov, RFmadogram instead.

  • RFoptions(grDefault=FALSE) returns to the old style of graphical device handling. Otherwise there is no handling.

  • C code is started to be parallelized.

  • Some new Multivariate RMmodels

  • New way of passing models, see RFformula, which allows connections (formulae) between parameters, e.g. one parameter value might be twice as large as another parameter value. Also dummy variables can be RMdeclared.

• Options getting obsolete (Version 3 and older)

  • oldstyle is becoming warn_oldstyle

  • newstyle is becoming warn_newstyle

  • newAniso is becoming warn_newAniso

  • ambiguous is becoming warn_ambiguous

  • normal_mode is becoming warn_normal_mode

  • colour_palette is becoming warn_colour_palette

• Changings in option names

  • several changes in RFoptions()$graphics in version 3.1.11

  • pdfnumber became in version 3.0.42 filenumber

  • pdfonefile became in version 3.0.42 onefile

  • pdffile became in version 3.0.42 file

  • tbmdim became in version 3.0.41 reduceddim

  • coord_units became in version 3.0.39 coordunits

  • new_coord_units became in version 3.0.39 new_coordunits

  • variab_units became in version 3.0.39 varunits

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

See Also

MajorRevisions, RandomFields.

Examples

## no examples given

Description

Circulant embedding is a fast simulation method for stationary (possibly anisotropic) Gaussian random fields on regular grids based on Fourier transformations. It is guaranteed to be an exact method for covariance functions with finite support, e.g. the spherical model. The method is admissible for any dimension apart from memory restrictions.
The simulation is performed on a torus which represents the bended grid. To remove wrong dependencies occuring at different borders of the grid which would be close on the torus, the simulation area is multiplied by a natural number. There is also a multivariate version of circulant embedding.

Cut-off embedding is a fast simulation method for stationary, isotropic Gaussian random fields on square lattices based on the standard RPcirculant method, so that exact simulation is garantueed for further covariance models, e.g. the RMwhittle model.

In fact, the circulant embedding is called with the cutoff hypermodel, see RMcutoff. Cutoff halfens the maximum number of elements models used to define the covariance function of interest (from 10 to 5).

Here, multiplicative models are not allowed (yet).
For details see RMcutoff.

Intrinsic embedding is a fast simulation method for intrinsically stationary, isotropic Gaussian random fields on square lattices based on the standard RPcirculant method, for further variogram models, e.g. RMfbm.

Note that the simulated random field is always non-stationary. In fact, the circulant embedding is called with the Intrinsic hypermodel, see RMintrinsic.

Here, multiplicative models are not allowed (yet).
For details see RMintrinsic.

Usage

RPcirculant(phi, boxcox, force, mmin, strategy,
 maxGB, maxmem,  tolIm, tolRe, trials, useprimes, dependent,
 approx_step, approx_maxgrid)

RPcutoff(phi, boxcox, force, mmin, strategy,
 maxGB, maxmem, tolIm, tolRe, trials, useprimes, dependent,
 approx_step, approx_maxgrid, diameter, a) 
 

RPintrinsic(phi, boxcox, force, mmin, strategy,
 maxGB, maxmem, tolIm, tolRe, trials, useprimes, dependent,
 approx_step, approx_maxgrid, diameter, rawR)

Arguments

Details

Here, the algorithms by Dietrich and Newsam are implemented.

Value

An object of class RMmodel.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

Circulant Embedding

• Chan, G. and Wood, A.T.A. (1997) An Algorithm for Simulating Stationary Gaussian Random Fields. Journal of the Royal Statistical Society: Series C (Applied Statistics), 46, 171–181.

• Dietrich, C. R. and G. N. Newsam (1993) A fast and exact method for multidimensional gaussian stochastic simulations. Water Resour. Res. 29(8), 2861–2869.

• Dietrich, C. R. and G. N. Newsam (1996) A fast and exact method for multidimensional Gaussian stochastic simulations: Extension to realizations conditioned on direct and indirect measurements Water Resour. Res. 32(6), 1643–1652.

• Dietrich, C. R. and Newsam, G. N. (1997) Fast and Exact Simulation of Stationary Gaussian Processes through Circulant Embedding of the Covariance Matrix. SIAM J. Sci. Comput. 18, 1088–1107.

• Wood, A. T. A. and Chan, G. (1994) Simulation of Stationary Gaussian Processes in $[0, 1]^d$. Journal of Computational and Graphical Statistics 3, 409–432.

Cutoff and Intrinsic

• Gneiting, T., Sevecikova, H, Percival, D.B., Schlather M., Jiang Y. (2006) Fast and Exact Simulation of Large Gaussian Lattice Systems in $R^2$: Exploring the Limits. J. Comput. Graph. Stat. 15, 483–501.

• Stein, M.L. (2002) Fast and exact simulation of fractional Brownian surfaces. J. Comput. Graph. Statist. 11, 587–599

See Also

Gaussian, RP

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

model <- RMstable(s=1, alpha=1.8)
x <- seq(-3,3,0.1)

z <- RFsimulate(model=RPcirculant(model), x=x, y=x, n=1)
plot(z)

model <- RMexp(var=10, s=2)
z <- RFsimulate(model=RPcirculant(model), 1:10)
plot(z)

model <- RMfbm(Aniso=diag(c(1,2)), alpha=1.5)
z <- RFsimulate(model, x=1:10, y=1:10)
plot(z)

Description

The random coin method (or dilution method) is a simulation method for stationary Gaussian random fields. It is based on the following procedure: For a stationary Poisson point process on ${\bf R}^d$ consider the random field

$Y(y) = \sum_{x\in X} f(y-x)$

for a function $f$. The covariance of $Y$ is proportional to the convolution

$C(h) = \int f(x)f(x+h) dx$

If the intensity of the Poisson point process increases, the random field $Y$ approaches a Gaussian random field with covariance function $C$.

Usage

RPcoins(phi, shape, boxcox, intensity, method) 

RPaverage(phi, shape, boxcox, intensity, method)

Arguments

Value

RPcoins returns an object of class RMmodel.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Lantuejoul, C. (2002) Geostatistical Simulation: Models and Algorithms. Springer.

See Also

Gaussian, RP, RPhyperplane, RPspectral, RPtbm.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

Description

Several constants are provided that might make the use of some functions easier, e.g. RFgetModelNames.

Value

RC_TYPE_NAMES = c( "tail correlation", "positive definite", "variogram", "negative definite", "point-shape function", "shape function", "trend", "distribution or shape", "of manifold type", "process", "method for Gauss process", "normed process (non-negative values with maximum value being 0 or 1)", "method for Brown-Resnick process", "Smith", "Schlather", "Poisson", "PoissonGauss", "distribution family", "interface", "mathematical operator", "other type")

RC_DOMAIN_NAMES = c("single variable", "kernel", "framework dependent", "submodel dependent", "parameter dependent", "<keep copy>", "mismatch")

RC_ISO_NAMES = c("isotropic", "space-isotropic", "vector-isotropic", "symmetric", "cartesian system", "gnomonic projection", "orthographic projection", "spherical isotropic", "spherical symmetric", "spherical system", "earth isotropic", "earth symmetric", "earth system", "cylinder system", "non-dimension-reducing", "framework dependent", "submodel dependent", "parameter dependent", "<internal keep copy>", "<mismatch>")

RC_MONOTONE_NAMES = c( "not set", "mismatch in monotonicity", "submodel dependent monotonicity", "previous model dependent monotonicity", "parameter dependent monotonicity", "not monotone", "monotone", "Gneiting-Schaback class", "normal mixture", "completely monotone", "Bernstein")

RC_ISOTROPIC gives the numerical code for option "isotropic"

RC_DOUBLEISOTROPIC gives the numerical code for option "space-isotropic"

RC_CARTESIAN_COORD gives the numerical code for option "cartesian system"

RC_GNOMONIC_PROJ gives the numerical code for the gnomonic projection, see also zenit in RFoptions.

RC_ORTHOGRAPHIC_PROJ gives the numerical code for the orthographic projection, see also zenit in RFoptions.

RC_EARTH_COORDS gives the numerical code for option "earth coordinates"

RC_EARTH_ISOTROPIC gives the numerical code for option "earth isotropic"

RC_SPHERICAL_COORDS gives the numerical code for option "earth coordinates"

RC_OPTIMISER_NAMES and RC_NLOPTR_NAMES give the names for the options optimiser and algorithm, respectively, RFfitoptimiser.

RC_LIKELIHOOD_NAMES = c("auto", "full", "composite", "tesselation") gives the names of the ML variants: (i) internal choice according to the number of data, (ii) full likelihood, (iii) (pairwise) composite likelihood, and (iv) composite likelihood based on a tessellation of the space.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

See Also

RF, RM, RP, RR, R., RFgetModelNames, RMmodelgenerator-class, RMtrafo.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

RC_ISO_NAMES
RC_ISO_NAMES[RC_ISOTROPIC:RC_CARTESIAN_COORD + 1]
## Not run: 
RFgetModelNames(isotropy=RC_ISO_NAMES[RC_ISOTROPIC:RC_CARTESIAN_COORD +
1])

## End(Not run)

Description

Generate an object of class RFsp from conventional objects.

Usage

conventional2RFspDataFrame(data, coords=NULL, gridTopology=NULL, n=1,
                           vdim=1, T=NULL, vdim_close_together)

Arguments

Value

Object of class RFspatialGridDataFrame, RFspatialPointsDataFrame, RFgridDataFrame or RFpointsDataFrame.

Author(s)

Alexander Malinowski, Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again
x <- 1:20
z <- RFsimulate(RMexp(), x, spConform=FALSE)
z2 <- conventional2RFspDataFrame(z, coord=x)
Print(z, z2)

Description

Implemented Coordinate Systems.

Implemented coordinate systems

• Cartesian coordinate system

• Earth coordinate systems
  The earth is considered as an ellipsoid; The first angle takes values in $[0, 360)$, the second angle takes values in $[-90, 90]$.

• Spherical coordinate systems
  The earth is considered as an ellipsoid; The first angle takes values in $[0, 2\pi)$, the second angle takes values in $[-\pi/2, \pi/2]$.

Transformations between the system

• Earth to cartesian
  The 3-dimensional resulting coordinates are either given in ‘km’ or in ‘miles’.

• Gnomonic and orthographic projections
  The 2-dimensional resulting coordinates are either given in ‘km’ or in ‘miles’. The projection direction is given by the zenit.

• Earth to spherical
  In this case the Earth is considered as a ball.

Cartesian systems cannot be transformed to earth or spherical coordinate systems, nor a spherical system to earth coordinates.

Options

coord_system

character. One of the values "auto", "cartesian", "earth"

If "auto", then the coordinates are considered as "cartesian" except the names of the given coordinates indicate a different system. Currently, only "longitude" and "latidute" (or abbreviations of them) are excepted as names for given coordinates and indicate an earth coordinate systems. See the examples below.

Default: "auto"

coordidx

integer vector of column numbers of the variables in the data frame. varidx can be set alternatively to coordnames. This parameter gives the coordinate columns in a data frame by starting column and ending column or the sequence. An NA in the second component means ‘until the end’.

coordnames

vector of characters that can be set alternatively to coordidx. This parameter gives the coordinate columns in a data frame by names. If it is "NA", then, depending on the context, either an error message is returned or it is assumed that the first columns give the coordinates.

coordunits

any string. If coordinate_system = "earth" and longitude and latitude are transformed to 3d cartesian coordinates, coordunits determines whether the radius is given in kilometers ("km") or miles ("miles"). If empty, then "km" is chosen.

Default: ""

new_coord_system

One of the values "keep", "cartesian", "earth", "plane".

1. "keep"
   The coord_system is kept (except an explicit transformation is given, see RMtrafo.

   Note that some classes of models, e.g. completely monotone functions and compactly supported covariance models with range less than $\pi$ are valid models on a sphere. In this case the models are considered as models on the sphere. See spherical models for lists.

2. "cartesian"
   If coord_system is "earth" the coordinates are transformed to cartesian coordinates before any model is considered.

3. "orthographic", "genomic"
   If coord_system is "earth" the locations are projected to a plane before any model is considered.

Default: "keep"

new_coordunits

internal and should not be set by the user.

Default: ""

polar_coord

logical. If FALSE the spherical coordinates agree with the earth coordinate parametrization, except that radians are used for spherical coordinates instead of degrees for the earth coordinates.

If TRUE the spherical coordinates signify polar coordinates.

Default : FALSE

varidx

integer vector of length 2. varidx can be set alternatively to varnames. This parameter gives the data columns in a data frame, either by starting column and ending column. An NA in the second component means ‘until the end’.

varnames

vector of characters that can be set alternatively to varidx. This parameter gives the data columns in a data frame by names.

if varnames equals "NA" then for keywords ‘data’, ‘value’ and ‘variable’ will be searched for keywords. If none of them are found, depending on the context, either an error message is returned or it is assumed that the last columns give the data.

varunits

vector of characters. The default units of the variables.

Default: ""

xyz_notation

logical or NA. Used by RMuser only.

NA : automatic choice (if possible)

FALSE : notation (x, y) should not be understood as kernel definition, not as xyz notation

TRUE: xyz notation used

zenit

two angles of the central projection direction for the gnomonic projection (https://en.wikipedia.org/wiki/Gnomonic_projection, https://de.wikipedia.org/wiki/Gnomonische_Projektion) and the orthographic projection, (https://en.wikipedia.org/wiki/Orthographic_projection_in_cartography, https://de.wikipedia.org/wiki/Orthografische_Azimutalprojektion).

If any(is.na(zenit)) then either the value of either of the components may not be NA, whose value will be denoted by $p$.

If $p=1$ then the mean of the locations is calculated; if $p=Inf$ then the mean of the range is calculated.

Default: c(1, NA)

References

Covariance models in a cartesian system

• Schlather, M. (2011) Construction of covariance functions and unconditional simulation of random fields. In Porcu, E., Montero, J.M. and Schlather, M., Space-Time Processes and Challenges Related to Environmental Problems. New York: Springer.

Covariance models on a sphere

• Gneiting, T. (2013) Strictly and non-strictly positive definite functions on spheres. Bernoulli, 19, 1327-1349.

Tail correlation function

• Strokorb, K., Ballani, F., and Schlather, M. (2014) Tail correlation functions of max-stable processes: Construction principles, recovery and diversity of some mixing max-stable processes with identical TCF. Extremes, Submitted.

See Also

RMtrafo, RFearth2cartesian, RPdirect, models valid on a sphere, RFoptions.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

z <- 1:4
x <- cbind(z, 0)
y <- cbind(0, z)
model <- RMwhittle(nu=0.5)
RFcov(model, x, y, grid=FALSE) ##  standard is (cartesian) model


## same as above, but explicit:
RFcov(model, x, y, grid=FALSE, coord_sys="cartesian") 


## model is not valid on a sphere; x,y coordinates are
## transformed from earth coordinates to spherical coordinates
RFcov(model, x, y, grid=FALSE, coord_sys="earth")


## now the scale is chosen such that the covariance
## values are comparable to those in the cartesian case
RFcov(RMS(model, s= 1 / 180 * pi), x, y, grid=FALSE,
      coord_sys="earth")


## projection onto a plane first. Then the scale is interpreted
## in the usual, i.e. cartesian, sense, i.e. the model does not
## really make sense
RFoptions(zenit = c(2.5, 2.5))
RFcov(model, x, y, grid=FALSE,
      coord_sys="earth", new_coord_sys="orthographic")


## again, here the scale is chosen to be comparable to the cartesian case
## here the (standard) units are [km]
(z1 <- RFcov(RMS(model, s= 6350 / 180 * pi), x, y, grid=FALSE,
             coord_sys="earth", new_coord_sys="orthographic"))


## as above, but in miles
(z2 <- RFcov(RMS(model, s= 6350 / 1.609344 / 180 * pi), x, y, grid=FALSE,
             coord_sys="earth", new_coord_sys="orthographic",
             new_coordunits="miles"))
stopifnot(all.equal(z1, z2))


## again, projection onto a plane first, but now using the
## gnomonic projection
## here the (standard) units are [km]
(z1 <- RFcov(RMS(model, s= 6350 / 180 * pi), x, y, grid=FALSE,
             coord_sys="earth", new_coord_sys="gnomonic"))

## as above, but in miles
(z2 <- RFcov(RMS(model, s= 6350 / 1.609344 / 180 * pi), x, y, grid=FALSE,
             coord_sys="earth", new_coord_sys="gnomonic",
             new_coordunits="miles"))
stopifnot(all.equal(z1, z2, tol=1e-5))

Description

Distribution families to specify random parameters in the model definition.

Details

See Bayesian Modelling for a less technical introduction to hierarchical modelling.

When simulating Gaussian random fields, the random parameters are drawn only once at the very beginning. So, if the argument n in RFsimulate is greater than 1 then n simulations conditional on a single realization of the random parameters are performed. See the examples below.

There are (simple) multivariate versions and additional versions to the distributions families implemented. Further, any distribution family defined in R can be used, see the examples below.

These functions will allow for Bayesian modelling. (Future project).

Implemented models

Note

The allowance of random parameters is a very recent, developing feature of RandomFields.

Future changings of the behaviour are not unlikely.

Note

A further random element is RMsign, which is an operator on shape functions. As an exception its name starts with RM and not with RR.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

See Also

RC, RF, RM, RP, Other models, RFdistr, RMmodelgenerator, R.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

## here, the scale is given by an exponential variable:
model <- RMgauss(scale=exp())
for (i in 1:4) {
  RFoptions(seed = i)
  # each leads to a simulation with a different scale parameter
  plot(model) ## random
  plot(RFsimulate(model, x=seq(0,10,0.1)))
  readline("press return")
}

# but here, all 4 simulations have the same (but random) scale:
plot(RFsimulate(model, x=seq(0,10,0.1), n=4)) 


## hierarchical models are also possible:
## here, the scale is given by an exponential variable whose
## rate is given by a uniform variable
model <- RMgauss(scale=exp(rate=unif()))
plot(model)
plot(RFsimulate(model, x=seq(0,10,0.1)))


## HOWEVER, the next model is deterministic with scale \code{e=2.718282}.
model <- RMgauss(scale=exp(1))
plot(model)
plot(RFsimulate(model, x=seq(0,10,0.1)))

Description

RPopitz defines an extremal t process.

Usage

RPopitz(phi, xi, mu, s, alpha)

Arguments

Details

The argument xi is always a number, i.e. $\xi$ is constant in space. In contrast, $\mu$ and $s$ might be constant numerical values or (in future!) be given by an RMmodel, in particular by an RMtrend model.
For $xi=0$, the default values of $mu$ and $s$ are $0$ and $1$, respectively. For $xi\not=0$, the default values of $mu$ and $s$ are $1$ and $|\xi|$, respectively, so that it defaults to the standard Frechet case if $\xi > 0$.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Davison, A.C., Padoan, S., Ribatet, M. (2012). Statistical modelling of spatial extremes. Stat. Science 27, 161-186.

• Opitz, T. (2012) A spectral construction of the extremal t process. arxiv 1207.2296.

See Also

RMmodel, RPgauss, maxstable, maxstableAdvanced.

Examples

## sorry, does not work savely yet

Description

RPschlather defines an extremal Gaussian process.

Usage

RPschlather(phi, tcf, xi, mu, s)

Arguments

Details

The argument xi is always a number, i.e. $\xi$ is constant in space. In contrast, $\mu$ and $s$ might be constant numerical values or (in future!) be given by an RMmodel, in particular by an RMtrend model.
For $xi=0$, the default values of $mu$ and $s$ are $0$ and $1$, respectively. For $xi\not=0$, the default values of $mu$ and $s$ are $1$ and $|\xi|$, respectively, so that it defaults to the standard Frechet case if $\xi > 0$.

The argument phi can be any random field for which the expectation of the positive part is known at the origin.

It simulates an Extremal Gaussian process $Z$ (also called “Schlather model”), which is defined by

$Z(x) = \max_{i=1}^\infty X_i \max(0, Y_i(x)),$

where the $X_i$ are the points of a Poisson point process on the positive real half-axis with intensity $c x^{-2} dx$, $Y_i \sim Y$ are iid stationary Gaussian processes with a covariance function given by phi, and $c$ is chosen such that $Z$ has standard Frechet margins. phi must represent a stationary covariance model.

Note

Advanced options are maxpoints and max_gauss, see RFoptions.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

See Also

RMmodel, RPgauss, maxstable, maxstableAdvanced.

Examples

RFoptions(seed=0, xi=0)
## seed=0: *ANY* simulation will have the random seed 0; set
##         RFoptions(seed=NA) to make them all random again

## xi=0: any simulated max-stable random field has extreme value index 0
x <- seq(0, 2,0.01)

## standard use of RPschlather (i.e. a standardized Gaussian field)
model <- RMgauss()
z1 <- RFsimulate(RPschlather(model), x)
plot(z1, type="l")

## the following refers to the generalized use of RPschlather, where
## any random field can be used. Note that 'z1' and 'z2' have the same
## margins and the same .Random.seed (and the same simulation method),
## hence the same values
model <- RPgauss(RMgauss(var=2))
z2 <- RFsimulate(RPschlather(model), x)
plot(z2, type="l")
all.equal(z1, z2) # true

## Note that the following definition is incorrect
try(RFsimulate(model=RPschlather(RMgauss(var=2)), x=x))


## check whether the marginal distribution (Gumbel) is indeed correct:
model <- RMgauss()
z <- RFsimulate(RPschlather(model, xi=0), x, n=100)
plot(z)
hist(unlist(z@data), 50, freq=FALSE)
curve(exp(-x) * exp(-exp(-x)), from=-3, to=8, add=TRUE)

Description

Here, some details of RFfit are given concerning the fitting of models for Gaussian random fields.

This documentation is far from being complete.

Maximum likelihood

The application of the usual maximum likelihood method and reporting the result is the default.

Least squares

The weighted least squares methods minimize

$\sum_{i} w_i (\hat \gamma(h_i) - \gamma(h_i))^2$

over all parametrized models of $\gamma$. Here, $i$ runs over all $N$ bins of the binned variogram $\hat \gamma$ and $h_i$ is the centre of bin $i$.

The following variants of the least squares methods, passed as sub.methods in RFfit are implemented:

'self'

$w_i = (\gamma(h_i))^{-2}$

'plain'

$w_i = 1$ for all $i$.

'sqrt.nr'

$w_i^2$ equals the number of points $n_i$ in bin $i$.

'sd.inv'

$1 / w_i$ equals the standard deviation of the variogram cloud within bin $i$.

'internal'

Three subvariants are implemented:

'internal1'

$w_i^2 = (N-i+1) n_i$

'internal2'

$w_i = N-i+1$

'internal3'

$w_i^2 = N-i+1$

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

See Also

RFfit, RFfit-class.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again
## See 'RFfit'.

Description

Here, all the methods (models) for simulating Gaussian random fields are listed.

Implemented models

Computing demand for simulations

Assume at $n$ locations in $d$ dimensions a $v$-variate field has to be simulated. Let

$f(n, d) = 2^d n \log(n)$

The following table gives in particular the time and memory needed for the specific simulation method.

Computing demand for interpolation

Assume $v$-variate data are given at $n$ locations in $d$ dimensions. To interpolate at $k$ locations RandomFields needs

Computing demand for conditional simulation

Assume $v$-variate data are given at $n$ locations $x_1,\ldots, x_n$ in $d$ dimensions. To conditionally simulate at $k$ locations $y_1,\ldots, y_k$, the computing demand equals the sum of the demand for interpolating and the demand for simulating on the $k+n$ locations. (Grid algorithms for simulating will apply if the $k$ locations $y_1,\ldots, y_k$ are defined by a grid and the $n$ locations $x_1,\ldots, x_n$ are a subset of $y_1,\ldots, y_k$, a situation typical in image analysis.)

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Chiles, J.-P. and Delfiner, P. (1999) Geostatistics. Modeling Spatial Uncertainty. New York: Wiley.

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

• Schlather, M. (2010) On some covariance models based on normal scale mixtures. Bernoulli, 16, 780-797.

• Schlather, M. (2011) Construction of covariance functions and unconditional simulation of random fields. In Porcu, E., Montero, J.M. and Schlather, M., Space-Time Processes and Challenges Related to Environmental Problems. New York: Springer.

• Yaglom, A.M. (1987) Correlation Theory of Stationary and Related Random Functions I, Basic Results. New York: Springer.

• Wackernagel, H. (2003) Multivariate Geostatistics. Berlin: Springer, 3nd edition.

See Also

RP, Other models, RMmodel, RFgetMethodNames, RFsimulateAdvanced.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

set.seed(1)
x <- runif(90, 0, 500)
z <- RFsimulate(RMspheric(), x)
z <- RFsimulate(RMspheric(), x, max_variab=10000)

Description

Here, the code of the paper on ‘Fast and Exact Simulation of Large Gaussian Lattice Systems in R2’ is given.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

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

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

## Figure 1 (pretty time consuming)
stabletest <- function(alpha, theta, size=512) {
  RFoptions(trials=1, tolIm = 1e-8, tolRe=0, force = FALSE,
            useprimes=TRUE, strategy=0, skipchecks=!FALSE,
             storing=TRUE)
  model <- RMcutoff(diameter=theta, a=1, RMstable(alpha=alpha))
  RFcov(dist=0, model=model, dim=2, seed=0)
  r <- RFgetModelInfo(modelname="RMcutoff", level=3)$storage$R_theor
  x <- seq(0, r, by= r / (size - 1)) * theta
  err <- try(RFsimulate(x, x, model=RPcirculant(model), n=0))
  return(if (class(err) == "try-error") NA else r)
}

alphas <- seq(1.52, 2.0, 0.02) 
thetas <- seq(0.05, 3.5, 0.05)

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

Description

RandomFields provides Bayesian modelling to some extend: (i) simulation of hierarchical models at arbitrary depth; (ii) estimation of the parameters of a hierarchical model of depth 1 by means of maximizing the likelihood.

Details

A Bayesian approach can be taken for scalar, real valued model parameters, e.g. the shape parameter nu in the RMmatern model. A random parameter can be passed through a distribution of an existing family, e.g. (dnorm, pnorm, qnorm, rnorm) or self-defined. It is passed without the leading letter d, p, q, r, but as a function call e.g norm(). This function call may contain arguments that must be named, e.g. norm(mean=3, sd=5).

Usage:

• exp() denotes the exponential distribution family with rate 1,

• exp(3) is just the scalar $e^3$ and

• exp(rate=3) is the exponential distribution family with rate $3$.

The family can be passed in three ways:

• implicitly, e.g. RMwhittle(nu=exp()) or

• explicitly through RRdistr, e.g. RMwhittle(nu=RRdistr(exp())).

• by use of RRmodels of the package.

The first is more convenient, the second more flexible and slightly safer.

Note

• While simulating any depth of hierarchical modelling is possible, estimation is currently restricted to one level of hierarchy.

• The effect of the distribution family varies between the different processes:

  • in max-stable fields and RPpoisson, a new realization of the prior distribution(s) is drawn for each shape function

  • in all other cases: a realization of the prior(s) is only drawn once. This effects, in particular, Gaussian fields with argument n>1, where all realizations are based on the same realization out of the prior distribution(s).

  Note that checking the validity of the arguments is rather limited for such complicated models, in general.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

See Also

RMmodelsAdvanced. For hierarchical modelling see RR.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again
## See 'RRmodels' for hierarchical models

## the following model defines the argument nu of the Whittle-Matern
## model to be an exponential random variable with rate 5.
model <- ~ 1 + RMwhittle(scale=NA, var=NA, nu=exp(rate=5)) + RMnugget(var=NA)

data(soil)
fit <- RFfit(model, x=soil$x, y=soil$y, data=soil$moisture, modus="careless")
print(fit)

Description

The Hyperplane method is a simulation method for stationary, isotropic random fields with exponential covariance function. It is based on a tessellation of the space by hyperplanes. Each cell takes a spatially constant value of an i.i.d. random variable. The superposition of several such random fields yields approximatively a Gaussian random field.

Usage

RPhyperplane(phi, boxcox, superpos, maxlines, mar_distr, mar_param ,additive)

Arguments

Value

RPhyperplane returns an object of class RMmodel.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

• Lantuejoul, C. (2002) Geostatistical Simulation: Models and Algorithms. Springer.

See Also

Gaussian, RP.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again
model <- RPhyperplane(RMexp(s=2), superpos=1)
x <- seq(0, 3, 0.04)
z <- RFsimulate(x=x, x, model=model, n=1)
plot(z)

Description

Method to simulate the Nugget effect. (Only for advanced users)

Usage

RPnugget(phi, boxcox, tol, vdim)

Arguments

Details

General

This method only allows RMnugget as a submodel.

Anisotropy

The method also allows for zonal nugget effects. Only there the argument tol becomes important. For the zonal nugget effect, the anisotropy matrix Aniso should be given in RMnugget. There, only the kernel of the matrix is important.

Points close together

The locations at a distance less than or equal to the RFoptions nugget.tol are considered as being identical. This strategy applies to the simulation method and the covariance function itself. Hence, the covariance function is only positive definite if nugget.tol=0.0. However, if the anisotropy matrix does not have full rank and nugget.tol=0.0, then the simulations are likely to be odd. The value of nugget.tol should be of order $10^{-15}$.

Repeated measurements

Measurement errors are mathematically not distinguishable from spatial nugget effects as long as measurements are not repeated at the very same space-time location. So there is no need to distinguish the spatial nugget effect from a measurement error. This is the default, see allow_duplicated_locations in RFoptions.

In case several measurement have been taken in single space-time locations, measurement errors can be separated from spatial noise. In this case RMnugget() models the measurement error (which corresponds to a non-stationary model in an abstract space) by default and the measurement error model cannot be extended beyond the given locations. On the other hand RMnugget(Ansio=something) and RMnugget(proj=something) model the spatial nugget effect (with and without zonal anisotropy in case Aniso has low and full rank respectively).

Role of RPnugget

Even for advanced users, there is no need to call RPnugget directly, as this is done internally when the RMnugget is involved in the covariance model.

Value

RPnugget returns an object of class RMmodel.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

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

See Also

Gaussian, RP, RPcoins, RPhyperplane, RPspectral, RPtbm.

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again


x <- y <- 1:2
xy <- as.matrix(expand.grid(x, y)) ## we get 4 locations
## Standard use of the nugget effect
model <- RMnugget(var = 100)
RFcovmatrix(model, x=xy)
as.vector(RFsimulate(model, x=x, y=x, tol=1e-10))

## zonal nugget effect, which is not along the axes
model <- RMnugget(Aniso=matrix(1, nr=2, nc=2))
RFcovmatrix(model, x=xy)
as.vector(RFsimulate(model, x=x, y=x, tol=1e-10))


## All the following examples refer to repeated measurements
RFoptions(allow_duplicated_locations = TRUE) 
(xy <- rbind(xy, xy)) ## now, the 4 locations are repeated twice 

## standard situation: the nugget is interpreted as measurement error:
model <- RMnugget()
RFcovmatrix(model, x=xy)
as.matrix(RFsimulate(model, x=xy)) 

## any anisotropy matrix with full rank: spatial nugget effect
model <- RMnugget(Aniso=diag(2))
RFcovmatrix(model, x=xy)
as.matrix(RFsimulate(model, x=xy))

## anisotropy matrix with lower rank: zonal nugget effect
model <- RMnugget(Aniso=matrix(c(1, 0, 0, 0), nc=2))
RFcovmatrix(model, x=xy)
as.matrix(RFsimulate(model, x=xy))

## same as before: zonal nugget effect
model <- RMnugget(Aniso=t(c(1,0)))
RFcovmatrix(model, x=xy)
as.matrix(RFsimulate(model, x=xy))

Description

These functions are internal and should not be used.

Usage

rfGenerateModels(package = "RandomFields", assigning,
                 RFpath ="~/svn/RandomFields/RandomFields",
                 RMmodels.file=paste(RFpath, "R/RMmodels.R", sep="/"),
                 PL = RFoptions()$basic$printlevel)

rfGenerateConstants(package="RandomFields", aux.package = "RandomFieldsUtils",
	   RFpath = paste0("~/svn/",package, "/", package),
           RCauto.file = paste(RFpath, "R/aaa_auto.R", sep="/"),
	   header.source =
	   c(if (length(aux.package) > 0) paste0("../../", aux.package,"/",
				 aux.package, "/src/Auto", aux.package, ".h"),
	     paste0("src/Auto",package,".h")),
	   c.source = paste0("src/Auto", package, ".cc"))

rfGenerateTest(package = "RandomFields", files = NULL,
               RFpath = paste0("~/svn/", package, "/", package))

rfGenerateMaths(package = "RandomFields",
			    files = "/usr/include/tgmath.h",
			    do.cfile = FALSE,
                            ## copy also in ../private/lit
                            Cfile = "QMath",
                            Rfile = "RQmodels",
                            RFpath = paste0("~/svn/",package,"/", package))


checkExamples(exclude = NULL, include=1:length(.fct.list),
               ask=FALSE, echo=TRUE, halt=FALSE, ignore.all = FALSE,
               path=package, package = "RandomFields", 
               read.rd.files=TRUE, local = TRUE, libpath = NULL,
               single.runs = FALSE)

ScreenDevice(height, width)

FinalizeExample()
StartExample(reduced = TRUE, save.seed = TRUE)
showManpages(path="/home/schlather/svn/RandomFields/RandomFields/man")

plotWithCircles(data, factor=1.0,
                xlim=range(data[,1])+c(-maxr,maxr),
                ylim=range(data[,2])+c(-maxr,maxr),
                col=1, fill=0, ...)

maintainers.machine()

Arguments

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

Examples

RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
##                   RFoptions(seed=NA) to make them all random again

## internal functions: no examples given


# for (i in dep.packages) cat(maintainer(i), "\n")

Description

Here, the code of the paper on ‘Analysis, simulation and prediction of multivariate random fields with package RandomFields’ is given.

Author(s)

Martin Schlather, [email protected], https://www.wim.uni-mannheim.de/schlather/

References

Schlather, M., Malinowski, A., Menck, P.J., Oesting, M. and Strokorb, K. (2015) Analysis, simulation and prediction of multivariate random fields with package RandomFields. Journal of Statistical Software, 63 (8), 1-25, url = ‘http://www.jstatsoft.org/v63/i08/’

See Also

weather, SS12, S10.

Examples

## Not run: 
               ###########################################
               ##  SECTION 4: UNCONDITIONAL SIMULATION  ##
               ###########################################

RFoptions(seed = 0, height = 4)
## seed=0:  *ANY* simulation will have the random seed 0; set
##          RFoptions(seed=NA) to make them all random again
## height : height of X11
## always_close_device=FALSE: the pictures are kept on the screen


## Fig. 1: linear model of coregionalization
M1 <- c(0.9, 0.6)
M2 <- c(sqrt(0.19), 0.8)
model <- RMmatrix(M = M1, RMwhittle(nu = 0.3)) + 
         RMmatrix(M = M2, RMwhittle(nu = 2))
x <- y <- seq(-10, 10,  0.2)
simu <- RFsimulate(model, x, y)
plot(simu)


## Fig. 2: Wackernagel's delay model
model <- RMdelay(RMstable(alpha = 1.9, scale = 2), s = c(4, 4))
simu <- RFsimulate(model, x, y)
plot(simu, zlim = 'joint')

## Fig. 3: extended Wackernagel's delay model
model <- RMdelay(RMstable(alpha = 1.9, scale = 2), s = c(0, 4)) + 
         RMdelay(RMstable(alpha = 1.9, scale = 2), s = c(4, 0))
simu <- RFsimulate(model, x, y)
plot(simu, zlim = 'joint')
plot(model, dim = 2, xlim = c(-5, 5), main = "Covariance function", 
     cex = 1.5, col = "brown")

## Fig. 4: latent dimension model
## MARGIN.slices has the effect of choosing the third dimension
##               as latent dimension
## n.slices has the effect of choosing a bivariate model
model <- RMgencauchy(alpha = 1.5, beta = 3)
simu <- RFsimulate(model, x, y, z = c(0,1))
plot(simu, MARGIN.slices = 3, n.slices = 2)


## Fig. 5: Gneiting's bivariate Whittle-Matern model
model <- RMbiwm(nudiag = c(1, 2), nured = 1, rhored = 1, cdiag = c(1, 5), 
                s = c(1, 1, 2))
simu <- RFsimulate(model, x, y)
plot(simu)


## Fig. 6: anisotropic linear model of coregionalization
M1 <- c(0.9, 0.6)
M2 <- c(sqrt(0.19), 0.8)
A1 <- RMangle(angle = pi/4, diag = c(0.1, 0.5))
A2 <- RMangle(angle = 0, diag = c(0.1, 0.5))
model <- RMmatrix(M = M1, RMgengneiting(kappa = 0, mu = 2, Aniso = A1)) +
         RMmatrix(M = M2, RMgengneiting(kappa = 3, mu = 2, Aniso = A2))
simu <- RFsimulate(model, x, y)
plot(simu)


## Fig. 7: random vector field whose paths are curl free
## A 4-variate field is simulated, where the first component
## refers to the potential field, the second and third component
## to the curl free vector field and the forth component to the
## field of sinks and sources
model <- RMcurlfree(RMmatern(nu = 5), scale = 4)
simu <- RFsimulate(model, x, y)
plot(simu, select.variables = list(1, 2 : 3, 4))
plot(model, dim = 2, xlim = c(-3, 3), main = "", cex = 2.3, col="brown") 


## Fig. 8: Kolmogorov's model of turbulence
## See the physical literature for its derivation and details
x <- y <- seq(-2, 2, len = 20)
model <- RMkolmogorov()
simu <- RFsimulate(model, x, y, z = 1)
plot(simu, select.variables = list(1 : 2), col = c("red"))
plot(model, dim = 3, xlim = c(-3, 3), MARGIN = 1 : 2, cex = 2.3,
     fixed.MARGIN = 1.0, main = "", col = "brown")



               ###########################################
               ## SECTION 5: DATA ANALYSIS              ##
               ###########################################

## get the data     
data("weather")
PT <- weather[ , 1 : 2]  ## full data set takes more than 
##                                     half an hour to be analysed
## transformation of earth coordinates to Euclidean coordinates
Dist.mat <- as.vector(RFearth2dist(weather[ , 3 : 4]))
All <- TRUE
\dontshow{if(RFoptions()$internal$examples_reduced){warning("reduced data set")
All <- 1:10
PT <- weather[All , 1 : 2] 
Dist.mat <- as.vector(RFearth2dist(weather[All , 3 : 4]))
}}



#################################
## model definition            ##
#################################
## bivariate pure nugget effect:
nug <- RMmatrix(M = matrix(nc = 2, c(NA, 0, 0, NA)), RMnugget())
## parsimonious bivariate Matern model
pars.model <- nug + RMbiwm(nudiag = c(NA, NA), scale = NA, cdiag = c(NA, NA),
                           rhored = NA)
## whole bivariate Matern model
whole.model <- nug + RMbiwm(nudiag = c(NA, NA), nured = NA, s = rep(NA, 3),
                            cdiag = c(NA, NA), rhored = NA)



#################################
## model fitting and testing   ## 
#################################
## 'parsimonious model'
## fitting takes a while ( > 10 min)
pars <- RFfit(pars.model, distances = Dist.mat, dim = 3, data = PT)
print(pars)
print(pars, full=TRUE)
RFratiotest(pars)
#RFcrossvalidate(pars, x = weather[All , 3 : 4], data = PT, full = TRUE)

## 'whole model'
## fitting takes a while ( > 10 min)
whole <- RFfit(whole.model, distances = Dist.mat, dim = 3, data = PT)
print(whole, full=TRUE)
RFratiotest(whole)
#RFcrossvalidate(whole, x = weather[All , 3 : 4], data = PT, full = TRUE)

## compare parsimonious and whole
RFratiotest(nullmodel = pars, alternative = whole)


#################################
## kriging                     ##
#################################
## First, the coordinates are projected on a plane
a <- colMeans(weather[All , 3 : 4]) * pi / 180
P <- cbind(c(-sin(a[1]), cos(a[1]), 0),
           c(-cos(a[1]) * sin(a[2]), -sin(a[1]) * sin(a[2]), cos(a[2])),
           c(cos(a[1]) * cos(a[2]), sin(a[1]) * cos(a[2]), sin(a[2])))
x <- RFearth2cartesian(weather[All , 3 : 4])
plane <- (t(x) %*%P)[ , 1 : 2]

## here, kriging is performed on a rectangular that covers
## the projected points above. The rectangular is approximated
## by a grid of length 200 in each direction.
n <- 200 
r <- apply(plane, 2, range)
dta <- cbind(plane, weather[All , 1 : 2])
z <- RFinterpolate(pars, data = dta, dim = 2,
                   x = seq(r[1, 1], r[2, 1], length = n),
                   y = seq(r[1, 2], r[2, 2], length = n),
                   varunits = c("Pa", "K"), spConform = TRUE)
plot(z)

## End(Not run)

Description

This man page documents some major changings in RandomFields.

Changes done in 3.1.0 (Summer 2015)

• full (univariate) trend modelling

• error in particular in RFfit corrected

• RFfit runs much faster now

• effects of modus operandi changed for estimating

Corrections done in 3.0.56 (Jan 2015)

• log Gauss field corrected (has been a log log Gauss field)

• RMconstant is now called RMfixcov

Corrections done in 3.0.55 (Jan 2015)

• Conditional simulation: several severe typos corrected.

Major Revision: changings from Version 2 to Version 3 (Jan 2014)

• S4 objects

  • RandomFields is now based on S4 objects using the package sp. The functions accept both sp objects and simple objects as used in version 2. See also above.

• Documentation

  • each model has now its own man page;

  • classes of models and functions are bundled in several pages: Covariance models start with RM, distribution families with RR, processes with RP, user functions with RF

  • the man pages of several functions are split into two parts:

    (i) a beginners man page which includes a link to

    (ii) man pages for advanced users

• Interfaces

  • The interfaces become simpler, at the same time more powerful than the functions in version 2. E.g., RFsimulate can perform unconditional simulation, conditional simulation and random imputing.

  • Only those arguments are kept in the functions that are considered as being absolutely necessary. All the other arguments can be included as options.

  • RFgui is an instructive interface based on tcl/tk, replacing the former ShowModels

• Inference for Gaussian random fields

  • RFfit has undergone a major revision. E.g.:

    (i) estimation of random effect models with spatial covariance structure

    (ii) automatic estimation of 10 and more arguments in multivariate and/or space-time models

  • RFvariogram is now based on an fft algorithm if the data are on a grid, even allowing for missing values.

  • RFratiotest has been added.

• Processes

  • Modelling of maxstable processes has been enhanced, including

    (i) the simulation of Brown-Resnick processes

    (ii) initial support of tail correlation functions;

  • Further processes chi2 processes, compound Poisson processes, binary processes added.

• Models

  • the formula notation for linear models may now be defined

  • Novel, user friendly definition of the covariance models

  • Multivariate and vector-valued random fields are now fully included

  • The user may now define his own functions, to some extend.

  • The trend allows for much more flexibility

  • Distributions may now be included which will be extended to Bayesian modelling in future.