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
Solved – global sensitivity analysis input and output matrix
monte carlosensitivity analysis
Related Solutions
I think it is a reasonable method. The difficulty is usually how to select the number of trials (which affects simulation time).
If you conduct $n$ trials, of which $N \leq n$ give a "positive" result, you can estimate the probability $p$ of positive result as its relative frequency,
$$\hat p = N/n.$$
This estimator is unbiased and has mean-square error, or variance,
$$\mathrm E[(\hat p-p)^2] = \mathrm E[(N/n-p)^2] = \frac{p(1-p)}{n},$$
as follows from noting that $N$ is a binomial random variable with parameters $n$, $p$. The root-mean-square (RMS) error, or standard deviation, is just the square root of this.
In order to assess whether an RMS error value is acceptable or not, that value should be compared with the true $p$. For example, an RMS of $0.01$ may be acceptable for estimating a $p=0.1$, but not for $p=0.001$ (the error would be ten times the true value). The usual approach for this is to normalize the error by dividing it by $p$. So the normalized RMS error is
$$\frac{\sqrt{\mathrm E[(\hat p-p)^2]}}{p} = \sqrt\frac{1-p}{np}.$$
This can be approximated as $1/\sqrt{np}$ for $p$ small. As you can see, to maintain a certain normalized error you need to simulate a number of times $n$ inversely proportional to $p$. But $p$ is unknown, so it's difficult to know which $n$ you need.
A solution is to use sequential estimation, in which $n$ is not fixed in advance, but is adaptively selected according to a certain stopping rule to guarantee that the estimation error does no exceed a predefined level. A standard sequential method is inverse binomial sampling (also called negative-binomial Monte Carlo), which consists in the following: continue simulating until a target number $N$ of positive results is achieved. So now $N$ is fixed, and the number of trials, $n$, becomes the random variable.
The nice aspect of this approach is that by selecting the target $N$ you can control the normalized error level that will be achieved, irrespective of the unknown $p$. That is, to each $N$ corresponds a certain guaranteed value of normalized error. This works whether you consider error defined as mean-square-error, mean-absolute-error, or in terms of a confidence interval. A description of the procedure in each case is given in the following papers (please bear with my self-referencing):
- Normalized RMS error: http://arxiv.org/abs/0809.4047 (or see IEEE Communications Letters, November 2009)
- Normalized mean-absolute error: http://oa.upm.es/3222/ (or see IEEE Transactions on Communications, November 2008)
- Normalized confidence interval: http://arxiv.org/abs/0809.2402 (or see Bernoulli Journal, May 2010)
The sensitivity analysis you suggest corresponds to examining the partial derivatives of the outputs with respect to the inputs. Say the output vector $y \in \mathbb{R}^m$ is given by $y= f(x)$ , where $x \in \mathbb{R}^d$ is the input vector and $f$ is the function the network implements. The Jacobian of the outputs w.r.t. the inputs is:
$$J_{ij}(x) = \frac{\partial}{\partial x_j} f_i(x)$$
The Jacobian gives the local rate of change of each output w.r.t. each input, so it tells us how $f$ will behave in response to infinitesimal perturbations. If we start with input $x$ and add an infinitesimal value $\Delta$ to the $j$th input, we expect the $i$th output to increase by $\Delta J_{ij}(x)$.
If $J_{ij}(x)$ has large magnitude, it means that output $i$ is sensitive to input $j$ in the vicinity of $x$. Because $f$ is, in general, nonlinear, this notion of sensitivity depends on the input; it may be large in some regions and near zero in others. If you want some kind of summary measure of how strongly the outputs depend on the inputs, you'd have to aggregate over multiple input values. For example, you could take the absolute value of the Jacobian, averaged over all inputs in the training set (which acts as a surrogate for the expected value w.r.t. the underlying distribution of inputs). Of course, this kind of summary will end up discarding information, so could be misleading in some circumstances.
You can use the chain rule to derive an expression for the Jacobian, similarly to how you'd derive the gradient of the loss function w.r.t. the parameters for use with backprop. You can also compute it using automatic differentiation, using a library like Theano, TensorFlow, etc. There's not much reason to perform finite differencing (i.e. actually simulate the perturbation and measure the change in output), unless the function your network implements is nondifferentiable (in which case the Jacobian doesn't exist).
A couple caveats: If the inputs have different units/scales than each other, the sensitivities will also have different units/scales, and can't be directly compared. Standardizing/scaling the inputs is one possible solution. It's also important to keep in mind is that this type of analysis tells us about the model itself, but not necessarily the underlying distribution that generated the data. For example, if two inputs are correlated, the model might end up using the first but not the second. In this case, we'd find that the sensitivity is high for the first input and low for the second, but should not conclude that the first input is inherently more important for predicting the output in general.
This article should be of interest.
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=5\cdot10^6$ nonlinear FEM analyses and storing all the corresponding results (?), you can use
fast99from the R packagesensitivity:xis 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 withx$Xyou can consult the docs of
sensitivityfor additional information. Once you have the $M$ runs of your code, and the corresponding vector of responsesy, you can update the objectxwithNow,
S <- x$D1/x$VandT <- 1- x$Dt/x$Vare respectively the first-order indices (or main effect indices) and the total order indices.fast99is reasonably fast, as far as non-model based methods for the estimate of Sobol indices go. However, it has two drawbacks:An alternative to
fast99, which is more accurate in the estimation of the smaller indices, but require thrice the number of runs, issobolowen. 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:and then, for 95%-confidence interval Monte Carlo estimates of the main effect and total order indices, you call
As before,
x$Xcontains your DOE. Updatexwithtell, and nowx$Sandx$Tstore respectively the first-order indices and the total order indices. Consult the vignette ofsensitivityfor 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$Xby 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):
https://inldigitallibrary.inl.gov/sti/6468109.pdf
https://inldigitallibrary.inl.gov/sti/6799579.pdf