spup (version 1.4-0)

genSample: Methods for generating Monte Carlo realizations from uncertain inputs.

Description

Methods for classes: "MarginalNumericSpatial", "MarginalScalar", "MarginalCategoricalSpatial", "JointNumericSpatial", "JointScalar". Function that runs Monte Carlo simulations depending on the type of uncertain object. Facilitates unconditional Gaussian simulation of errors for spatially auto-correlated residuals, as well as random and stratified random sampling if no spatial auto-correlation is included.

Usage

genSample(
  UMobject,
  n,
  samplemethod,
  p = 0,
  asList = TRUE,
  debug.level = 1,
  ...
)

Value

A Monte Carlo sample of the variables of interest. If asList = TRUE returns list of all samples as lists.

Arguments

UMobject

an uncertain object to sample from, output of defineUM() or defineMUM().

n

integer, number of Monte Carlo realizations.

samplemethod

a string, "ugs", "randomSampling", "stratifiedSampling", "lhs" ("lhs" currently not in use).

p

A vector of quantiles. Optional. Only required if sample method is "stratifiedSampling" or "lhs".

asList

logical. If asList = TRUE returns list of all samples as a list. If asList = FALSE returns samples in a format of distribution parameters in UMobject.

debug.level

integer; set gstat internal debug level, see below for useful values. If set to -1 (or any negative value), a progress counter is printed.

...

Additional parameters that may be passed, e.g. in the "ugs" method. See examples.

Author

Kasia Sawicka, Stefan van Dam, Gerard Heuvelink

Details

Sampling methods:

"ugs" Unconditional Gaussian simulation of spatially auto-correlated and/or cross-correlated errors.

"randomSampling" Sampling multivariate distribution using eigenvalue decomposition (based on 'mvtnorm' package).

"stratifiedSampling" Number of samples (n) must be dividable by the number of quantiles to assure that each quantile is evenly represented.

"lhs" Not implemented yet. Sampling method for at least two uncertain inputs. The uncertain.object is then a list of two or more. It uses a stratified sampling method to generate inputs for the latin hypercube algorithm.

NOTE. Version 1.3-1 includes bug fixing related to derivation of cross-correlation matrix for multivariate uncertainty propagation analysis.

Examples

set.seed(12345)

### ------------------- "MarginalNumericSpatial" -------------------
# load data
data(dem30m, dem30m_sd)

# "ugs" method example
dem_crm <- makeCRM(acf0 = 0.78, range = 321, model = "Exp")
demUM <- defineUM(uncertain = TRUE, distribution = "norm", 
                   distr_param = c(dem30m, dem30m_sd), crm = dem_crm)
                   
# toy example
dem_sample <- genSample(UMobject = demUM, n = 2, samplemethod = "ugs", nmax = 4, asList = FALSE)
str(dem_sample)
# any meaningful Monte Carlo analysis should have normally much larger number of runs               
if (FALSE) {
dem_sample <- genSample(UMobject = demUM, n = 100, samplemethod = "ugs", nmax = 20, asList = FALSE)
str(dem_sample)
}

# "startifiedSampling" method example
demUM <- defineUM(uncertain = TRUE, distribution = "norm", distr_param = c(dem30m, dem30m_sd))
# toy example
dem_sample <- genSample(UMobject = demUM, n = 5, samplemethod = "stratifiedSampling", p = 0:5/5)
# any meaningful Monte Carlo analysis should have normally much larger number of runs
if (FALSE) {
dem_sample <- genSample(UMobject = demUM, n = 100, samplemethod = "stratifiedSampling", p = 0:5/5)
str(dem_sample)
}

# Examples with rasters
# (raster with auto-correlation)
data(OC, OC_sd)
OC_crm <- makeCRM(acf0 = 0.6, range = 1000, model = "Sph")
OC_UM <- defineUM(TRUE, distribution = "norm", distr_param = c(OC, OC_sd), crm = OC_crm, id = "OC")
class(OC_UM)
# toy example
some_sample <- genSample(OC_UM, n = 2, "ugs", nmax = 4)
some_sample
# any meaningful Monte Carlo analysis should have normally much larger number of runs
if (FALSE) {
some_sample <- genSample(OC_UM, n = 50, "ugs", nmax = 24)
some_sample
}


### ----------------------- "MarginalScalar" -----------------------
# example with normal distribution
scalarUM <- defineUM(uncertain = TRUE, distribution = "norm", distr_param = c(10, 1))
scalar_sample <- genSample(scalarUM, n = 10, samplemethod = "randomSampling")


### ----------------- "MarginalCategoricalSpatial" -----------------
# load data
data(woon)
woonUM <- defineUM(TRUE, categories = c(1,2,3), cat_prob = woon[, c(4:6)])
woon_sample <- genSample(woonUM, 10, asList = FALSE)
class(woon_sample)
str(woon_sample@data)


### -------------------- "JointNumericSpatial" ----------------------
# load data
data(OC, OC_sd, TN, TN_sd)

# define marginal UMs
OC_crm <- makeCRM(acf0 = 0.6, range = 1000, model = "Sph")
OC_UM <- defineUM(TRUE, distribution = "norm", distr_param = c(OC, OC_sd), crm = OC_crm, id = "OC")
TN_crm <- makeCRM(acf0 = 0.4, range = 1000, model = "Sph")
TN_UM <- defineUM(TRUE, distribution = "norm", distr_param = c(TN, TN_sd), crm = TN_crm, id = "TN")

# define joint UM
soil_prop <- list(OC_UM, TN_UM)
mySpatialMUM <- defineMUM(soil_prop, matrix(c(1,0.7,0.7,1), nrow=2, ncol=2))
class(mySpatialMUM)

# sample - "ugs" method
# toy example
my_cross_sample <- genSample(mySpatialMUM, n = 2, "ugs", nmax = 4, asList = TRUE)
class(my_cross_sample)
# any meaningful Monte Carlo analysis should have normally much larger number of runs
if (FALSE) {
my_cross_sample <- genSample(mySpatialMUM, n = 100, "ugs", nmax = 24, asList = TRUE)
class(my_cross_sample)
}


### ------------------------- "JointScalar" ------------------------
scalarUM <- defineUM(uncertain = TRUE, distribution = "norm",
                     distr_param = c(1, 2), id="Var1")                
scalarUM2 <- defineUM(uncertain = TRUE, distribution = "norm",
                      distr_param = c(3, 2), id="Var2")
scalarUM3 <- defineUM(uncertain = TRUE, distribution = "norm",
                      distr_param = c(10, 2.5), id="Var3")                
myMUM <- defineMUM(UMlist = list(scalarUM, scalarUM2, scalarUM3), 
               matrix(c(1,0.7,0.2,0.7,1,0.5,0.2,0.5,1), nrow = 3, ncol = 3))
my_sample <- genSample(myMUM, n = 5, samplemethod = "randomSampling", asList = FALSE)
my_sample