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
Consider a function of $p$ random variables
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
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
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
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)
x$X contains your DOE. Update
tell, and now
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;
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):