sobolrec: Recursive estimation of Sobol' indices

Description

sobolrec implements a recursive version of the procedure introduced by Tissot & Prieur (2015) using two replicated nested designs. This function estimates either all first-order indices or all closed second-order indices at a total cost of $2 \times N$ model evaluations where $N$ is the size of each replicated nested design.

Usage

sobolrec(model=NULL, factors, layers, order, precision, method=NULL, tail=TRUE, 
          ...)
# S3 method for sobolrec
ask(x, index, ...)
# S3 method for sobolrec
tell(x, y = NULL, index, ...)
# S3 method for sobolrec
print(x, ...)
# S3 method for sobolrec
plot(x, ylim = c(0,1), ...)

Value

sobolrec returns a list of class "sobolrec", containing all the input arguments detailed before, plus the following components:

call: the matched call.
X: a data.frame containing the design of experiments (row concatenation of the two replicated designs).
y: a list of the response used at each step.
V: a list of the model variance estimated at each step.
S: a list of the Sobol' indices estimated at each step.
steps: the number of steps performed.
N: the size of each replicated nested design.

Arguments

model

a function, or a model with a predict method, defining the model to analyze.

factors

an integer giving the number of factors, or a vector of character strings giving their names.

layers

If order=1, a vector specifying the respective sizes of each layer (see "Details"). If order=2, an integer specifying the size of all layers.

order

an integer specifying which indices to estimate: 1 for first-order indices, 2 for closed second-order indices.

precision

a vector containing:

the target precision for the stopping criterion.
the number of steps for the stopping criterion (must be greater than 1).

tail

a boolean specifying the method used to choose the number of levels of the orthogonal array (see "Warning messages").

method

If order=2, a character specifying the method to construct the orthogonal arrays (see "Details"):

"al" for the algebraic method
"ar" for the accept-reject method

Set to NULL if order=1.

x

a list of class "sobolrec" storing the state of the sensitivity study (parameters, data, estimates).

index

an integer specifying the step of the recursion

y

the model response.

ylim

y-coordinate plotting limits.

...

any other arguments for model which are passed unchanged each time it is called.

Warning messages

"The value entered for layers is not the square of a prime number. It has been replaced by: ": When order=2, the value of layers must be the square of a prime power number. This warning message indicates that it was not the case and the value has been replaced depending on tail. If tail=TRUE (resp. tail=FALSE) the new value of layers is equal to the square of the prime number preceding (resp. following) the square root of layers.
"The value entered for layers is not satisfying the constraint. It has been replaced by: ": the value $N$ for layers must satisfied the constraint $N \geq (d-1)^{2}$ where $d$ is the number of factors. This warning message indicates that N was replaced by the square of the prime number following (or equals to) $d-1$.

Details

For first-order indices, layers is a vector: $$\left(s_1, ...,s_m \right)$$ specifying the number $m$ of layers of the nested design whose respective size are given by: $$\prod_{i=1}^{k-1} s_i, \ k=2, ...,m+1$$ For closed second-order indices, layers directly specifies the size of all layers.

For each Sobol' index $S$ the stopping criterion writes: $$\mid S_{l-1}-S_{l} \mid < \epsilon$$ This criterion is tested for the last $l_0$ steps (including the current one). $\epsilon$ and $l_0$ are respectively the target precision and the number of steps of the stopping criterion specified in precision.

sobolrec uses either an algebraic or an accept-rejet method to construct the orthogonal arrays for the estimation of closed second-order indices. The algebraic method is less precise than the accept-reject method but offers more steps when the number of factors is small.

sobolrec automatically assigns a uniform distribution on [0,1] to each input. Transformations of distributions (between U[0,1] and the wanted distribution) have to be performed before the call to tell().

References

A.S. Hedayat, N.J.A. Sloane and J. Stufken, 1999, Orthogonal Arrays: Theory and Applications, Springer Series in Statistics.

L. Gilquin, E. Arnaud, H. Monod and C. Prieur, 2021, Recursive estimation procedure of Sobol' indices based on replicated designs, Computational and Applied Mathematics, 40:1--23.

Examples

Run this code

	# \donttest{
# Test case: the non-monotonic Sobol g-function

# The method of sobol requires 2 samples
# (there are 8 factors, all following the uniform distribution on [0,1])

# first-order indices estimation
x <- sobolrec(model = sobol.fun, factors = 8, layers=rep(2,each=15), order=1,
              precision = c(5*10^(-2),2), method=NULL, tail=TRUE)
print(x)

# closed second-order indices estimation
x <- sobolrec(model = sobol.fun, factors = 8, layers=11^2, order=2,
              precision = c(10^(-2),3), method="al", tail=TRUE)
print(x)


# Test case: dealing with external model 
# put in comment because of bug with ask use !

#x <- sobolrec(model = NULL, factors = 8, layers=rep(2,each=15), order=1,
#              precision = c(5*10^(-2),2), method=NULL, tail=TRUE)
#toy <- sobol.fun
#k <- 1
#stop_crit <- FALSE
#while(!(stop_crit) & (k

Run the code above in your browser using DataLab