GeoSim: Simulation of Gaussian and non-Gaussian random fields.

In GeoModels: Procedures for Gaussian and Non Gaussian Geostatistical (Large) Data Analysis

Simulation of Gaussian and non-Gaussian random fields.

Description

Simulates a realization of a Gaussian or non-Gaussian spatial, spatio-temporal, or spherical random field, and Gaussian bivariate random fields, for a specified covariance or correlation model. The covariance parameters are supplied through param; available correlation models are documented in GeoCovmatrix.

Usage

GeoSim(coordx=NULL, coordy=NULL,coordz=NULL, coordt=NULL, coordx_dyn=NULL, corrmodel, 
 distance="Eucl", grid=FALSE, method="cholesky", 
 model='Gaussian', n=1, param,anisopars=NULL,radius=1, 
 sparse=FALSE,X=NULL,spobj=NULL,nrep=1,progress=TRUE,check.duplicates=FALSE)

Arguments

coordx	May be NULL when coordx_dyn or spobj is supplied. Otherwise, a numeric (d \times 2)-matrix or (d \times 3)-matrix. Coordinates on a sphere for a fixed radius radius are passed in lon/lat format expressed in decimal degrees.
coordy	A numeric vector giving 1-dimension of spatial coordinates; Optional argument, the default is NULL.
coordz	A numeric vector giving 1-dimension of spatial coordinates; Optional argument, the default is NULL.
coordt	A numeric vector giving the temporal coordinates at which the field is simulated. Optional argument; the default is NULL, in which case a spatial random field is expected. Temporal coordinates may be irregularly spaced; temporal lags are computed from the supplied coordinate values.
coordx_dyn	For dynamic simulation sites, a list with one two- or three-column coordinate matrix per element of coordt. The output at time t follows the row order of coordx_dyn[[t]]. See GeoModels-spacetime-ordering.
corrmodel	String; the name of a correlation model, for the see GeoCovmatrix for the list of implemented correlation models.
distance	String; the name of the spatial distance. The default is Eucl, the Euclidean distance. See GeoFit for details.
grid	Logical; if FALSE (the default) the data are interpreted as spatial or spatio-temporal realisations on a set of non-equispaced spatial sites (irregular grid).
method	String; the type of matrix decomposition used in the simulation. Default is cholesky. The other possible choices is svd.
model	String; the type of RF and therefore the densities associated to the likelihood objects. Gaussian is the default, see the Section Details.
n	Positive integer size parameter. For Binomial it may be scalar or site-specific; for direct Negative Binomial it is the single common number r of successes, with Geometric corresponding to r=1.
param	A list of parameter values required in the simulation procedure of random fields, see Examples.
anisopars	A list of two elements "angle" and "ratio" i.e. the anisotropy angle and the anisotropy ratio, respectively.
radius	Numeric; a value indicating the radius of the sphere when using the great circle distance. Default value is 1.
sparse	Logical; if TRUE then cholesky decomposition is performed using sparse matrices algorithms (spam packake). It should be used with compactly supported covariance models.FALSE is the default.
X	Numeric design matrix for the mean. For fixed locations, rows are ordered by time blocks, with all sites at the first time followed by all sites at the second time. For dynamic locations, supply either a stacked matrix in temporal-block order or a list with X[[t]] aligned with coordx_dyn[[t]].
spobj	An object of class sp or spacetime. For space-time objects, the current sp2Geo() conversion uses sequential temporal indices and does not preserve irregular original time spacing; to retain irregular temporal distances, use the explicit-coordinate interface with numeric coordt.
nrep	Numeric; Numbers of indipendent replicates.
progress	Logic; If TRUE then a progress bar is shown.
check.duplicates	Logical. If TRUE, perform a fast exact scan for duplicated observation locations at this user entry point. The default is FALSE, so no duplicate-location scan is imposed. Internal bootstrap, cross-validation, and refitting calls do not repeat the scan.

Details

Bivariate simulation is currently implemented only for model = "Gaussian"; other bivariate marginal models are rejected explicitly. Binary and Bernoulli are aliases of a binomial field with n = 1; Geom and Geometric are aliases of a negative-binomial field with n = 1. Direct simulation of model = "Beta2" is not implemented by GeoSim; use GeoSimCopula with an explicit copula instead. The function stops explicitly rather than returning a latent Gaussian draw. Models whose names contain Gaussian_misp_ are inferential working models and are rejected as data-generating mechanisms. Other model names that are valid elsewhere in GeoModels but do not have a direct simulator are also rejected explicitly.

Several direct non-Gaussian constructions use an integer number of independent latent Gaussian fields. Consequently, Gamma requires integer shape; Beta requires integer shape1 and shape2; and for StudentT, SkewStudentT, and TwoPieceStudentT, 1/df must be an integer at least 3. For PoissonGamma and PoissonGammaZIP, 2*shape must be a positive integer. These parameters are never rounded silently. These restrictions apply to the direct latent-field constructions in GeoSim; they do not apply to quantile-transformed copula margins in GeoSimCopula.

nrep must be a positive integer. For binomial and negative-binomial fields, n must contain positive integers and may be scalar or have one value per simulated observation. The returned param component preserves the parameter list supplied by the user rather than the standardized latent-Gaussian working parameters.

For univariate simulation, the location parameter is \mu=X\beta, with coefficients named mean, mean1, and so on in the order of the columns of X. If X=NULL, the simulation is intercept-only and uses scalar mean. Alternatively, param$mean may be a vector with one value per simulated observation; this external mean is mutually exclusive with X.

For the Tukey transformed-Gaussian models, Tukeyh requires 0 <= tail < 0.5, while Tukeyh2 requires both 0 <= tail1 < 0.5 and 0 <= tail2 < 0.5. In Tukeyh2, tail1 is the right-tail parameter and tail2 is the left-tail parameter. Zero is an allowed boundary and recovers the Gaussian transformation on the corresponding side. SinhAsinh requires a strictly positive tail.

Value

Returns an object of class GeoSim. An object of class GeoSim is a list containing at most the following components:

bivariate	Logical:TRUE if the Gaussian RF is bivariate, otherwise FALSE;
coordx	A d-dimensional vector of spatial coordinates;
coordy	A d-dimensional vector of spatial coordinates;
coordz	A d-dimensional vector of spatial coordinates;
coordt	A t-dimensional vector of temporal coordinates;
coordx_dyn	A list of dynamical (in time) spatial coordinates;
corrmodel	The correlation model; see GeoCovmatrix.
data	The simulated data. For fixed-location space-time simulation this is a matrix with times in rows and sites in columns; for dynamic locations it is a list with one vector per time, aligned with coordx_dyn. See GeoModels-spacetime-ordering.
distance	The type of spatial distance;
method	The method of simulation
model	The type of RF, see GeoFit.
n	The Binomial number of trials; for direct Negative Binomial, the common number r of successes.
numcoord	The number of spatial coordinates;
numtime	The number the temporal realisations of the RF;
param	The parameter list supplied to the simulation call;
radius	The radius of the sphere if coordinates are passed in lon/lat format;
spacetime	TRUE if spatio-temporal and FALSE if spatial RF;
nrep	The number of indipendent replicates;

Spatio-temporal ordering

With fixed locations, a simulated space-time realization is returned as a T \times N matrix: rows correspond to coordt and columns to rows of coordx. The corresponding internal and X row order is c(t(data)), i.e. time then site.

With dynamic locations, the simulated realization is a list of length T. Element data[[t]] has one value per row of coordx_dyn[[t]], in the same row order. Replicates, when requested, contain objects with this same fixed or dynamic layout. See GeoModels-spacetime-ordering.

Author(s)

Moreno Bevilacqua, [e-mail],https://sites.google.com/view/moreno-bevilacqua/home, Víctor Morales Oñate, [e-mail], https://sites.google.com/site/moralesonatevictor/, Christian Caamaño-Carrillo, [e-mail],https://www.researchgate.net/profile/Christian-Caamano

See Also

GeoCovmatrix for covariance matrix construction, GeoFit for parameter estimation, GeoSimcond and GeoSimapprox for conditional and approximate simulation.

Examples

library(GeoModels)


################################################################
###
### Example 1. Simulation of a spatial Gaussian RF 
### with Matern and Generalized Wendland correlations
###############################################################

# Define the spatial-coordinates of the points:
x <- runif(500);y <- runif(500)
coords=cbind(x,y)
set.seed(261)
# Simulation of a spatial Gaussian RF with Matern correlation function
data1 <- GeoSim(coordx=coords, corrmodel="Matern", param=list(smooth=0.5,
 mean=0,sill=1,scale=0.4/3,nugget=0))$data

set.seed(261)
data2 <- GeoSim(coordx=coords, corrmodel="GenWend", param=list(smooth=0,
 power2=4,mean=0,sill=1,scale=0.4,nugget=0))$data
opar=par(no.readonly = TRUE)
par(mfrow=c(1,2))
if (requireNamespace("fields", quietly = TRUE)) {
 fields::quilt.plot(coords, data1, main = "Matern", xlab = "", ylab = "")
}
if (requireNamespace("fields", quietly = TRUE)) {
 fields::quilt.plot(coords, data2, main = "Wendland", xlab = "", ylab = "")
}
par(opar)
 

################################################################
###
### Example 2. Simulation of a spatial geometric RF 
### with underlying Wend0 correlation
###
################################################################

# Define the spatial-coordinates of the points:
x <- runif(800);y <- runif(800)
coords <- cbind(x,y)
set.seed(251)
# Simulation of a spatial Binomial RF:
sim <- GeoSim(coordx=coords, corrmodel="Wend0",
 model="BinomialNeg",n=1,sparse=TRUE,
 param=list(nugget=0,mean=0,scale=.2,power2=4))

if (requireNamespace("fields", quietly = TRUE)) {
 fields::quilt.plot(
  coords, sim$data, nlevel = max(sim$data),
  col = terrain.colors(max(sim$data + 1))
 )
}

################################################################
###
### Example 3. Simulation of a spatial Weibull RF
### with underlying Matern correlation on a regular grid
###
###############################################################
# Define the spatial-coordinates of the points:
x <- seq(0,1,0.032)
y <- seq(0,1,0.032)
set.seed(261)
# Simulation of a spatial Gaussian RF with Matern correlation function
data1 <- GeoSim(x,y,grid=TRUE, corrmodel="Matern",model="Weibull", 
 param=list(shape=1.2,mean=0,scale=0.3/3,nugget=0,smooth=0.5))$data
if (requireNamespace("fields", quietly = TRUE)) {
 fields::image.plot(x, y, data1, main = "Weibull RF", xlab = "", ylab = "")
}

################################################################
###
### Example 4. Simulation of a spatial t RF
### with with underlying Generalized Wendland correlation 
###
###############################################################
# Define the spatial-coordinates of the points:
x <- seq(0,1,0.03)
y <- seq(0,1,0.03)
set.seed(268)
# Simulation of a spatial Gaussian RF with Matern correlation function
data1 <- GeoSim(x,y,grid=TRUE, corrmodel="GenWend",model="StudentT", sparse=TRUE,
 param=list(df=1/4,mean=0,sill=1,scale=0.3,nugget=0,smooth=1,power2=5))$data
if (requireNamespace("fields", quietly = TRUE)) {
 fields::image.plot(
  x, y, data1, col = terrain.colors(100), main = "Student-t RF",
  xlab = "", ylab = ""
 )
}


################################################################
###
### Example 5. Simulation of a sinhasinh RF
### with underlying Wend0 correlation.
###
###############################################################

# Define the spatial-coordinates of the points:
x <- runif(500, 0, 2)
y <- runif(500, 0, 2)
coords <- cbind(x,y)
set.seed(261)
corrmodel="Wend0"
# Simulation of a spatial Gaussian RF:
param=list(power2=4,skew=0,tail=1,
 mean=0,sill=1,scale=0.2,nugget=0) ## gaussian case
data0 <- GeoSim(coordx=coords, corrmodel=corrmodel,
 model="SinhAsinh", param=param,sparse=TRUE)$data
plot(density(data0),xlim=c(-7,7))

param=list(power2=4,skew=0,tail=0.7,
 mean=0,sill=1,scale=0.2,nugget=0) ## heavy tails
data1 <- GeoSim(coordx=coords, corrmodel=corrmodel,
 model="SinhAsinh", param=param,sparse=TRUE)$data
lines(density(data1),lty=2)

param=list(power2=4,skew=0.5,tail=1,
 mean=0,sill=1,scale=0.2,nugget=0) ## asymmetry
data2 <- GeoSim(coordx=coords, corrmodel=corrmodel,
 model="SinhAsinh", param=param,sparse=TRUE)$data
lines(density(data2),lty=3)

################################################################
###
### Example 6. Simulation of a bivariate Gaussian RF
### with bivariate Matern correlation model
###
###############################################################

# Define the spatial-coordinates of the points:
x <- runif(500, 0, 2)
y <- runif(500, 0, 2)
coords <- cbind(x,y)

# Simulation of a bivariate spatial Gaussian RF:
# with a separable Bivariate Matern
param=list(mean_1=4,mean_2=2,smooth_1=0.5,smooth_2=0.5,smooth_12=0.5,
 scale_1=0.12,scale_2=0.1,scale_12=0.15,
 sill_1=1,sill_2=1,nugget_1=0,nugget_2=0,pcol=0.5)
data <- GeoSim(coordx=coords,corrmodel="Bi_matern",
 param=param)$data
opar=par(no.readonly = TRUE)
par(mfrow=c(1,2))
if (requireNamespace("fields", quietly = TRUE)) {
 fields::quilt.plot(
  coords, data[1, ], col = terrain.colors(100), main = "1",
  xlab = "", ylab = ""
 )
}
if (requireNamespace("fields", quietly = TRUE)) {
 fields::quilt.plot(
  coords, data[2, ], col = terrain.colors(100), main = "2",
  xlab = "", ylab = ""
 )
}
par(opar)


################################################################
###
### Example 7. Simulation of a spatio temporal Gaussian random field.
### observed on fixed location sites with double Matern correlation 
###
###############################################################



coordt=1:5

# Define the spatial-coordinates of the points:
x <- runif(50, 0, 2)
y <- runif(50, 0, 2)
coords <- cbind(x,y)

param<-list(nugget=0,mean=0,scale_s=0.2/3,scale_t=2/3,sill=1,smooth_s=0.5,smooth_t=0.5)
data <- GeoSim(coordx=coords, coordt=coordt, corrmodel="Matern_Matern",
 param=param)$data
dim(data)

################################################################
###
### Example 8. Simulation of a spatio temporal Gaussian random field.
### observed on dynamical location sites with double Matern correlation 
###
###############################################################

# Define the dynamical spatial-coordinates of the points:

coordt=1:5
coordx_dyn=list()
maxN=30
set.seed(8)
for(k in 1:length(coordt))
{
NN=sample(1:maxN,size=1)
x <- runif(NN, 0, 1)
y <- runif(NN, 0, 1)
coordx_dyn[[k]]=cbind(x,y)
}
coordx_dyn

param<-list(nugget=0,mean=0,scale_s=0.2/3,scale_t=2/3,sill=1,smooth_s=0.5,smooth_t=0.5)
data <- GeoSim(coordx_dyn=coordx_dyn, coordt=coordt, corrmodel="Matern_Matern",
 param=param)$data
## spatial realization at first temporal instants
data[[1]]
## spatial realization at third temporal instants
data[[3]]




################################################################
###
### Example 9. Simulation of a Gaussian RF 
### with a Wend0 correlation in the north emisphere of the planet earth
### using geodesic distance
###############################################################
distance="Geod";radius=6371

NN=3000 ## total point on the sphere on lon/lat format
set.seed(80)
coords=cbind(runif(NN,-180,180),runif(NN,0,90))
## Set the wendland parameters
corrmodel <- "Wend0"
param<-list(mean=0,sill=1,nugget=0,scale=1000,power2=3)
# Simulation of a spatial Gaussian RF on the sphere
#set.seed(2)
data <- GeoSim(coordx=coords,corrmodel=corrmodel,sparse=TRUE,
 distance=distance,radius=radius,param=param)$data
#require(globe)
#globe::globeearth(eye=place("newyorkcity"))
#globe::globepoints(loc=coords,pch=20,col = cm.colors(length(data),alpha=0.4)[rank(data)])

GeoModels documentation built on Sept. 23, 2026, 5:07 p.m.