I am performing nonlinear analysis in ABAQUS software, with 5 input variables. For each realization of the input vector, I run the code and I get a scalar response value. I am using the Monte Carlo method to generate a sample for the input vector. I generated a Monte Carlo sample of size $N=500000$ for the input vector, and correspondigly I obtained a sample of the same size for the response by running the code. How can I perform a global sensitivity analysis (GSA) without a predefined function for it? Thanks

**Contents**hide

#### Best Answer

Consider a function of $p$ random variables

$$Y=f(X_1,dots,X_p)$$

When the output is neither a linear nor a monotone function of the outputs, **but** $X_1,dots,X_p$ are independent, **and** f is square-integrable, then a good way to estimate the global sensitivity of $Y$ to $X_1,dots,X_p$ is to use the Sobol indices. The definition of Sobol indices is related to the ANOVA decomposition(note this is a permalink). Assuming that you are only interested in the main effect indices and total order indices, and assuming that you are really going to run $N=5cdot10^6$ nonlinear FEM analyses and storing all the corresponding results (?), you can use `fast99`

from the R package `sensitivity`

:

`N <- 10^6 p <- 5 library(sensitivity) x <- fast99(model = NULL, factors = p, n = N, q = "qnorm", q.arg = list(mean = 0, sd = 1)) `

`x`

is an object which, among the other things, stores the design you will need to run, i.e., the $M=pN$ realizations of the input vector, corresponding to the $M$ runs of your code. You can access them with

`x$X`

you can consult the docs of `sensitivity`

for additional information. Once you have the $M$ runs of your code, and the corresponding vector of responses `y`

, you can update the object `x`

with

`tell(x, y) `

Now, `S <- x$D1/x$V`

and `T <- 1- x$Dt/x$V`

are respectively the first-order indices (or main effect indices) and the total order indices.

`fast99`

is reasonably fast, as far as non-model based methods for the estimate of Sobol indices go. However, it has two drawbacks:

- it only gives you a point estimate, but not a confidence interval
- it can be inaccurate when the Sobol indices are small in magnitude w.r.t. 1, even though with millions of runs this won't be a problem.

An alternative to `fast99`

, which is more accurate in the estimation of the smaller indices, but require thrice the number of runs, is `sobolowen`

. This returns a Monte-Carlo based estimator of the Sobol indices, thus you get a confidence interval along with the point estimates. In this case you generate 3 design matrices this way:

`X1 <- data.frame(matrix(runif(p * N), nrow = N)) X2 <- data.frame(matrix(runif(p * N), nrow = N)) X3 <- data.frame(matrix(runif(p * N), nrow = N)) `

and then, for 95%-confidence interval Monte Carlo estimates of the main effect and total order indices, you call

`x <- sobolowen(model = NULL, X1, X2, X3, nboot = 100) `

As before, `x$X`

contains your DOE. Update `x`

with `tell`

, and now `x$S`

and `x$T`

store respectively the first-order indices and the total order indices. Consult the vignette of `sensitivity`

for many other estimators of the Sobol indices, with varying computational efforts.

PS note that here I assumed all your variables to be standard normal. I guess $X_1,dots,X_p$ are design parameters, whose distributions you assumed normal, but which are likely not standard normal. This means that to obtain your actual DOE, you need to transform each column of `x$X`

by multiplying it for the standard deviation of variable $X_i$, and then add the mean of $X_i$.

**EDIT**: I just got the abstract of this paper in my feeds;

http://www.ans.org/pubs/journals/nt/a_38989

It discusses exactly your problem (though I assume their ABAQUS model will be more complex than yours), so you should give it a look. If you cannot access the paper from the above link, these two pdfs seem quite similar, even if a little older (1 year):