Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models

Variance-based sensitivity analysis has been extensively studied and successfully developed in the context of deterministic simulators. For a deterministic model ${\mathcal{M}}_{d}$, Sobol’ indices quantify the contribution of each input variable $\left\{X_{i},i=1,\ldots,M\right\}$, or combination thereof, to the variance of the model output $Y={\mathcal{M}}_{d}(\bm{X})$.

In this paper, we assume that $X_{i}$’s are mutually independent. Let us split the input vector into two subsets $\bm{X}=(\bm{X}_{{\bm{\mathsf{u}}}},\bm{X}_{\sim{\bm{\mathsf{u}}}})$, where ${\bm{\mathsf{u}}}\subset\left\{1,\ldots,M\right\}$ and $\sim{\bm{\mathsf{u}}}$ is the complement of ${\bm{\mathsf{u}}}$, i.e., $\sim{\bm{\mathsf{u}}}=\left\{1,\ldots,M\right\}\setminus{\bm{\mathsf{u}}}$. From the total variance theorem, the variance of the output can be decomposed as

The first-order and total Sobol’ indices introduced by Sobol’ [8] and Homma and Saltelli [20] for the subset of input variables $\bm{X}_{{\bm{\mathsf{u}}}}$ are defined by

Higher-order Sobol’ indices can be defined with the help of $S_{{\bm{\mathsf{u}}}}$. For example, the second-order or two-factor interaction Sobol’ index of $X_{1}$ and $X_{2}$ is given by

where we denote $S_{\left\{i\right\}}$ by $S_{i}$ for the sake of simplicity.

In the context of stochastic simulators defined in Equation 1, the input variables alone do not determine the value of the output. Iooss and Ribatet [9] extend $\bm{X}$ by adding the internal source of randomness represented by latent variables $\bm{Z}$, which turns the stochastic simulator into a deterministic one. In this case, all the input variables are gathered in $\left(\bm{X}_{{\bm{\mathsf{u}}}},\bm{X}_{\sim{\bm{\mathsf{u}}}},\bm{Z}\right)$, and thus the Sobol’ indices in Equation 3 can be naturally extended to

Note that $S_{{\bm{\mathsf{u}}}}$ has the same expression as in the case of deterministic simulators, but $S_{T_{\bm{\mathsf{u}}}}$ contains the additional variables $\bm{Z}$ [14]. Similarly, higher-order Sobol’ indices corresponding to interactions among components of $\bm{X}$ are defined in the same way as deterministic simulators, whereas interactions between components of $\bm{X}$ and $\bm{Z}$ involve $\bm{Z}$ in their definition. Since this is a direct extension, the Sobol’ indices defined in Equation 5 are referred to as classical Sobol’ indices in the sequel.

Another way to extend Sobol’ indices to stochastic simulators is to first eliminate the internal randomness by representing the response random variable $Y(\bm{x})$ by some summarizing statistical quantity, called here quantity of interest (QoI) denoted by ${\rm QoI}(\bm{x})$, such as the mean value $m(\bm{x})$, variance $v(\bm{x})$ [9], $\alpha$-quantile $q_{\alpha}(\bm{x})$ [21] and differential entropy $h(\bm{x})$ [12]. As a result, the stochastic simulator is reduced to a deterministic function ${\rm QoI}(\bm{x})$, and we can calculate the associated QoI-based Sobol’ indices as follows:

A third extension is defined by considering a stochastic simulator as a random field. For a fixed internal randomness $\bm{Z}(\omega)=\bm{z}$, the stochastic simulator is a deterministic function of the input variables, which corresponds to a trajectory. Hence, the associated Sobol’ indices are well-defined. Yet, they are random variables because of their dependence on $\bm{Z}$ [10], which results in the trajectory-based Sobol’ indices:

where the indices of the variance operators correspond to those variables to which these operators apply.

The three types of Sobol’ indices introduced above have different nature and focus. The classical Sobol’ indices defined in Equation 5 treat the latent variables $\bm{Z}$ as a set of separate input variables. As a result, indices of this type treat $\bm{Z}$ in the same way as $\bm{X}$. The first-order index $S_{{\bm{\mathsf{u}}}}$ indicates how much the output variance can be reduced (in expectation) if we can fix the value of $\bm{X}_{{\bm{\mathsf{u}}}}$. Besides, the classical Sobol’ indices can also quantify the influence of the intrinsic randomness as well as its interactions with input variables.

The QoI-based Sobol’ indices defined in Equation 6 help study a specific statistical quantity of the model response, which is a deterministic function of the inputs. Using a summary quantity to represent the random output might lead to a loss of information [10], unless this quantity itself is of interest. For example, we may want to find the variable(s) that has the largest effect(s) on the 95% quantile (i.e., $q_{\alpha}(\bm{X})$ for $\alpha=0.95$) of the model response. However, the importance (ranking) of the inputs $X_{i}$’s can be quite different depending on the choice of the QoI.

Unlike the previous two types of indices, trajectory-based Sobol’ indices presented in Equation 7 are random variables. This is because the latent variables $\bm{Z}$ and the input parameters $\bm{X}$ are treated differently: conditioned on a given $\bm{Z}=\bm{z}_{0}$, the stochastic simulator reduces to a deterministic function of $\bm{X}$, and we can calculate the associated (classical) Sobol’ indices. To evaluate the probability distribution of these trajectory-based Sobol’ indices requires that the same random seeds can be explicitly fixed in the simulator when running it for different values of $\bm{x}$. In a sense, trajectory-based Sobol’ indices emphasize the variation of the trajectory of stochastic simulators. They typically show the importance of each input variable in terms of its contribution to the variability of trajectories.

For stochastic simulators, the question “which input variable has the strongest effect?” is rather vague and cannot be answered by a single type of index. The analyst should properly pose the problem and select the appropriate sensitivity measures. If one aims at reducing the variance of the model output, the classical Sobol’ indices are of interest. If one is interested in some summary QoIs (which is often the case for applications, e.g., quantiles in reliability analysis), the QoI-based indices are more appropriate. Finally, if one can control the internal randomness and is interested in the variability of the model output for fixed intrinsic stochasticity, trajectory-based indices should be selected.

To further illustrate this discussion, let’s consider the stochastic simulator

where $X_{1}$, $X_{2}$ and $Z$ are independent random variables following standard normal distributions. The classical first-order Sobol’ indices are $S_{1}=0.5$ and $S_{2}=0$. Therefore, if we want to primarily reduce the variance of $Y$, $X_{1}$ should be investigated. For the response mean function $m(\bm{X})=X_{1}$, we have $S^{m}_{1}=1$ and $S^{m}_{2}=0$, which indicates that $X_{1}$ contributes fully to the variation of the mean function. In contrast, for the variance function $v(\bm{X})=X^{2}_{2}$, $S^{v}_{1}=0$ and $S^{v}_{2}=1$ reveals a different order. Regarding the trajectory-based indices, we have $S^{{\rm traj}}_{1}(Z)=\frac{1}{1+Z^{2}}$ and $S^{{\rm traj}}_{2}(Z)=\frac{Z^{2}}{1+Z^{2}}$ which are two random variables (due to the randomness in $Z$). The probability distribution of the two variables characterize how the randomness in $\bm{X}$ affects the function ${\mathcal{M}}_{s}(\bm{X},z)$ with $z$ being fixed as a single realization of $Z$. In summary, even for this simple toy example, the various indices provide different conclusions.

It is worth remarking that the classical Sobol’ indices $S_{{\bm{\mathsf{u}}}}$ in Equation 5 share some common properties with the mean-based Sobol’ indices in Eq. 6, when we consider the mean function ${\rm QoI}(\bm{x})\stackrel{{\scriptstyle\text{def}}}{{=}}m(\bm{x})={\mathbb{E}}\left[Y\mid\bm{X}=\bm{x}\right]$:

According to the law of total expectation, we have

As a result, $S_{{\bm{\mathsf{u}}}}$ can be rewritten as

This implies that both classical Sobol’ indices $S_{{\bm{\mathsf{u}}}}$ and mean-based Sobol’ indices $S^{m}_{{\bm{\mathsf{u}}}}$ provide the same ranking, as the numerators of Equations 8 and 10 are identical. However, it is worth emphasizing that $S_{{\bm{\mathsf{u}}}}$ is not equal to $S^{m}_{{\bm{\mathsf{u}}}}$ and they are measuring different quantities, since the denominator of Equation 10 is ${\rm Var}\left[Y\right]$ but that of Equation 8 is ${\rm Var}\left[m(\bm{X})\right]$.

As a summary for all three extensions, a stochastic simulator is essentially transformed into a reduced deterministic model at a certain stage. The classical Sobol’ indices include the latent variables as a part of the inputs. The QoI-based Sobol’ indices rely on a deterministic representation. The trajectory-based indices are random variables whose statistical properties can be studied at the cost of repeating a standard Sobol’ analysis for different realizations of the latent variables separately. As a result, the three types of sensitivity indices can be estimated by modifying only slightly the standard methods based on Monte Carlo simulation developed for deterministic simulators [4].

The classical Sobol’ indices involving $\bm{Z}$ (e.g., $S_{\bm{Z}}$, $S_{T_{{\bm{\mathsf{u}}}}}$ in Eq. 5) and the trajectory-based Sobol’ indices require controlling the latent variables $\bm{Z}$. In practical computations, this is achieved by fixing the random seed in the computational model [14, 10]. However, for certain types of stochastic simulators, or when the data are generated by physical experiments, it may be difficult to control or even identify $\bm{Z}$. For the sake of general applicability, we focus in this paper only on Sobol’ indices that can be estimated by manipulating $\bm{X}$, that is, the QoI-based Sobol’ indices and, to some extent, the classical Sobol’ indices.

Using Monte Carlo simulations to estimate these indices requires evaluating the simulator for various realizations of the input vector. In addition, it is generally necessary to evaluate the function ${\rm QoI}(\bm{x})$ for calculating the associated QoI-based Sobol’ indices. However, this function is not directly accessible due to the intrinsic randomness. Therefore this function is usually estimated by using replications; i.e., for each realization $\bm{x}$, the simulator is run many times, and ${\rm QoI}(\bm{x})$ is estimated from the output samples. Both factors call for a large number of model runs, which becomes impracticable for costly models. Therefore, the use of surrogate models is unavoidable.

In the sequel, we present the generalized lambda model as a stochastic surrogate. Such a model emulates the response distribution conditioned on $\bm{X}=\bm{x}$, which fully characterizes the statistical dependence between the inputs and output. Therefore, it can be used to estimate the considered Sobol’ indices.

The generalized lambda distribution is a flexible distribution family, which is designed to approximate many common distributions [18], e.g., normal, lognormal, Weibull and generalized extreme value distributions. A GLD is defined by its quantile function $Q(u)$ with $u\in[0,1]$, that is, the inverse of the cumulative distribution function $Q(u)=F^{-1}(u)$. In this paper, we consider the GLD of the Freimer-Kollia-Mudholkar-Lin (FKML) family [22] with four parameters, whose quantile function is defined by

where $\lambda_{1}$ is the location parameter, $\lambda_{2}$ is the scaling parameter, and $\lambda_{3}$ and $\lambda_{4}$ are the shape parameters. $\lambda_{2}$ is required to be positive to produce valid quantile functions (i.e., $Q$ being non-decreasing on $[0,1]$). Based on the quantile function defined in Equation 11, the probability density function (PDF) of a random variable $Y$ following a GLD is given by

where $\text{1}_{[0,1]}$ is the indicator function. A closed-form expression of $Q^{-1}$ is not available for arbitrary values of $\lambda_{3}$ and $\lambda_{4}$. Therefore, evaluating the PDF for a given $y$ usually requires solving the nonlinear equation (12) numerically.

GLDs cover a wide range of shapes determined by $\lambda_{3}$ and $\lambda_{4}$ [16]. For instance, $\lambda_{3}=\lambda_{4}$ produces symmetric PDFs, and $\lambda_{3}<\lambda_{4}$ ($\lambda_{3}>\lambda_{4}$) results in left-skewed (resp. right-skewed) distributions. Moreover, $\lambda_{3}$ and $\lambda_{4}$ are closely linked to the support and the tail behaviors of the corresponding PDF. More precisely, $\lambda_{3}$ and $\lambda_{4}$ control the left and right tail, respectively. Whereas $\lambda_{3}>0$ implies that the PDF support is left-bounded, $\lambda_{3}\leq 0$ implies that the distribution has a lower infinite support. Similarly, $\lambda_{4}>0$ implies that the PDF support is right-bounded, whereas it is $+\infty$ for $\lambda_{4}\leq 0$. In addition, for $\lambda_{3}<0$ ($\lambda_{4}<0$), the left (resp. right) tail decays asymptotically as a power law. Hence, GLDs can also provide fat-tailed distributions. The reader is referred to [16, 22] for a longer presentation of GLDs.

Consider a deterministic model ${\mathcal{M}}_{d}$ that maps a set of input parameters $\bm{x}=\left(x_{1},x_{2},\ldots,x_{M}\right)\in{\mathcal{D}}_{\bm{X}}\subset{\mathbb{R}}^{M}$ to the output $y\in{\mathbb{R}}$. Under the assumption that $Y={\mathcal{M}}_{d}(\bm{X})$ has finite variance, ${\mathcal{M}}_{d}$ belongs to the Hilbert space ${\mathcal{H}}$ of square-integrable functions with respect to the inner product $\langle u,v\rangle_{{\mathcal{H}}}={\mathbb{E}}\left[u(\bm{X})v(\bm{X})\right]=\int_{{\mathcal{D}}_{\bm{X}}}u(\bm{x})v(\bm{x})f_{\bm{X}}(\bm{x})\mathrm{d}\bm{x}$. If the joint distribution $f_{\bm{X}}$ satisfies certain conditions [23], the simulator ${\mathcal{M}}_{d}$ admits a spectral representation in terms of orthogonal polynomials:

where $\psi_{\bm{\alpha}}$ is a multivariate polynomial basis function indexed by $\bm{\alpha}\in{\mathbb{N}}^{M}$, and $c_{\bm{\alpha}}$ denotes the associated coefficient. The orthogonal basis can be obtained by using tensor products of univariate polynomials, each of which is orthogonal with respect to the probability measure $f_{X_{i}}(x_{i})\mathrm{d}x_{i}$ of the $i$-th variable $X_{i}$:

Details about the construction of this generalized polynomial chaos expansion can be found in [24, 25].

The PCE defined in Equation 13 contains an infinite sum of terms. However, in practice, it is only feasible to use a finite series as an approximation. To this end, truncation schemes are adopted to select a set of basis functions defined by a finite subset ${\mathcal{A}}\subset{\mathbb{N}}^{M}$ of multi-indices. A typical scheme is the hyperbolic (a.k.a. $q$-norm) truncation scheme [26] given by

where $p$ is the maximum degree of polynomials, and $q\leq 1$ defines the quasi-norm $\|\cdot\|_{q}$. Note that with $q=1$, we obtain the full basis of total degree less than $p$, which corresponds to the standard truncation scheme.

Because of the flexibility of GLDs, we assume that the response random variable $Y(\bm{x})$ of a stochastic simulator for a given input vector $\bm{x}$ can be well approximated by a GLD. Hence, the associated distribution parameters $\bm{\lambda}$ are functions of the input variables:

Under appropriate conditions discussed in Section 3.2, each component of $\bm{\lambda}(\bm{x})$ can be represented by a series of orthogonal polynomials. Because $\lambda_{2}(\bm{x})$ is required to be positive (see Section 3.1), the associated polynomial chaos representation is built on the natural logarithmic transform $\log\left(\lambda_{2}(\bm{x})\right)$. This results in the following approximations:

where ${\mathcal{A}}_{l}$ ($l=1,2,3,4$) is a finite set of selected basis functions for $\lambda_{l}$, and $c_{l,\bm{\alpha}}$’s are the coefficients. For the purpose of clarity, we explicitly express $\bm{c}$ in the PC approximations $\lambda^{\operatorname*{PC}}_{l}\left(\bm{x};\bm{c}\right)$ to emphasize that $\bm{c}$ are the model parameters yet to be estimated from data.

We assume that our costly stochastic simulator is evaluated once for each point $\bm{x}^{(i)}$ of the experimental design ${\mathcal{X}}$, and the associated model response $y^{(i)}$ is collected in ${\mathcal{Y}}$:

As already mentioned (and as emphasized by the notation $\omega^{(i)}$), no replications are required, and we do not control the random seed. To construct a GLaM from the available data $({\mathcal{X}},{\mathcal{Y}})$, both the truncated sets $\bm{{\mathcal{A}}}$ of basis functions and the coefficients $\bm{c}$ shall be determined. In this section, we summarize the method proposed in [17], which is designed to achieve both purposes without the need for replications.

Sometimes prior knowledge is available to set the basis functions. For example, when working with a standard linear regression problem, the data is supposed to be generated by

where $\epsilon$ has mean zero and is independent of $\bm{X}$. This case can be treated within the GLaM framework as follows: ${\mathcal{A}}_{1}$ contains the constant and linear term, and ${\mathcal{A}}_{2}$, ${\mathcal{A}}_{3}$, ${\mathcal{A}}_{4}$ contain only a constant term. Note that such a GLaM allows estimating the distribution of $\epsilon$, which is not required to be normal, whereas the usual linear regression framework assumes normally distributed $\epsilon$.

However, in general there is no prior knowledge that would help select $\bm{{\mathcal{A}}}$. Thus, we make the following assumptions to find appropriate hyperbolic truncation schemes defined in Equation 15 for each $\lambda_{i},i=1,\ldots,4$:

1. (A1)

   The response distribution of $Y(\bm{x})$ can be well-approximated by a generalized lambda distribution;

2. (A2)

   The shape of this distribution smoothly varies as a function of $\bm{x}$, so that the parameters $\lambda_{i}(\bm{x})$ can be well approximated by a low-order PCE

Because the shape of a GLD is controlled by $\lambda_{3}$ and $\lambda_{4}$, the associated hyperbolic truncation schemes ${\mathcal{A}}^{p,q,M}$ can be set with a small value of $p$, say $p=1$.

Moreover, the parameters $\lambda_{1}(\bm{x})$ and $\lambda_{2}(\bm{x})$ mainly affect the variation of the mean $m(\bm{x})$ and of the variance $v(\bm{x})$ as a function of the input $\bm{x}$, respectively. As a result, they may require possibly larger degree $p$. To this end, we modify the feasible generalized least-squares (FGLS) [27] to find suitable truncation schemes for the mean and variance function modeled as

Basically, FGLS iterates between a weighted least-square problem (WLS) to fit the mean function, and an ordinary least-square (OLS) analysis to estimate the variance function.

The details of the modified FGLS are presented in Algorithm 1. In this algorithm, the inputs $\bm{p}_{1}$ and $\bm{q}_{1}$ stand for the set of candidate degrees and $q$-norms that are tested to expand $\lambda_{1}(\bm{x})$, respectively. The same notation apply to $\bm{p}_{2}$ and $\bm{q}_{2}$ for $\lambda_{2}(\bm{x})$. Indeed, because of the low cost of least-square analysis, various combinations of $p$ and $q$ are tested for both $\lambda_{1}(\bm{x})$ and $\lambda_{2}(\bm{x})$.

More precisely, AOLS denotes adaptive ordinary least-squares with degree and $q$-norm adaptivity [28, 29]. This algorithm first builds a series of PCEs, each of which is obtained by applying ordinary least-squares with a truncation scheme ${\mathcal{A}}^{p,q,M}$ defined by a combination of $p\in\bm{p}$ and $q\in\bm{q}$. Then, it selects the PCE, therefore the associated truncation scheme, with the smallest leave-one-out errors (see [29] for details). $\operatorname{WLS}$ denotes the use of weighted least-squares, which takes the estimated variance $\hat{\bm{v}}$ as weight to re-estimate $\bm{c}_{m}$. In this procedure, the truncation set ${\mathcal{A}}_{m}$ for $m(\bm{x})$ is selected only once (before the loop), whereas a set of truncation schemes $\left\{{\mathcal{A}}^{i}_{v}:i=1,\ldots,N_{\rm FGLS}\right\}$ is obtained.

We finally select the one with the smallest leave-one-out error as the final truncated set ${\mathcal{A}}_{v}$ for $v(\bm{x})$. The number of iterations $N_{\rm FGLS}$ is defined by the user, typically $N_{\rm FGLS}=$ 5–10. After applying Algorithm 1, we set ${\mathcal{A}}_{1}={\mathcal{A}}_{m}$ and ${\mathcal{A}}_{2}={\mathcal{A}}_{v}$.

Once the basis functions are selected, we use the maximum (conditional) likelihood estimator to estimate $\bm{c}$:

where $\mathsf{L}\left(\bm{c}\right)$ is the conditional negative log-likelihood

with $f^{{\rm GLD}}$ being the probability density function of the GLD defined in Equation 12.

The advantages of the proposed estimator are twofold. On the one hand, the simulator is required to be evaluated only once (but not limited to one) on each point of the experimental design. Thereby, replications are not necessary (yet possible), and the method is versatile in this respect. On the other hand, if the underlying computational model can be exactly represented by a GLaM for a specific choice of $\bm{c}$, the maximum likelihood estimator is consistent (see proof in [17]).

In practice, the evaluation of $\mathsf{L}(\bm{c})$ is not straightforward because the PDF of generalized lambda distributions does not have an explicit form: it is necessary to solve nonlinear equations as shown in Eq. 12. Nevertheless, the nonlinear function $Q(u;\bm{\lambda})$ is monotonic and defined on $[0,1]$. Therefore, we proposed using the bisection method [30] to efficiently solve the nonlinear equations.

The first example is defined as follows:

where $X_{1},X_{2}\sim{\mathcal{U}}(0,2\pi)$, $X_{3}\sim{\mathcal{U}}(0.25,0.75)$ are independent input variables, and $Z\sim{\mathcal{N}}(0,1)$ denotes the latent variable that introduces the intrinsic randomness. The response distribution is a shifted lognormal distribution: the shift is equal to $\sin(x_{1})+7\sin^{2}\left(x_{2}\right)$ and the lognormal distribution is parameterized by ${\mathcal{L}}{\mathcal{N}}\left(\frac{x_{1}}{\pi},x_{3}\right)$. As a result, this stochastic simulator has a nonlinear location function and a strong heteroskedastic effect: the variance varies between 0.069 and 72.35. Besides, this example has a mild signal-to-noise ratio ${\rm SNR}=1.4$. This implies that the input variables can explain around 58% of the total variance of $Y$ (i.e., $S_{\left\{1,2,3\right\}}=0.58$).

Figure 1(b) compares the PDFs predicted by a GLaM built on an experimental design of $N=1{,}000$ with the reference response PDFs of the simulator. The results show that the developed algorithm correctly identifies the shape of the underlying shifted lognormal distribution. Moreover, the PDF supports and tails are also accurately approximated.

We consider the differential entropy $h(\bm{x})$ [12] as the QoI in this example. Because the analytical response distribution and entropy are known, we investigate the convergence of GLaM in terms of the conditional quantile function estimation Equation 29 and the entropy estimation Equation 28. The size of experimental design varies among $N\in\left\{250;500;1{,}000;2{,}000;4{,}000\right\}$. Note that the entropy of a GLD does not have a closed form. Therefore, we use $10^{4}$ Monte Carlo samples to estimate this quantity of a GLaM for each $\bm{x}$ in the test set.

In addition, we consider another model where we approximate the response distribution with a normal distribution. The mean and variance (as functions of $\bm{x}$) for such an approximation are chosen as the true mean and variance of the original. In other words, this model represents the “oracle” of Gaussian-type mean-variance models.

The results are summarized in Figure 2. The proposed method exhibits a clear convergence with respect to $N$ for both $Q(u;\bm{x})$ and $h(\bm{x})$ estimations. We observe in Figure 2(a) that the decay of $\varepsilon_{Q}$ has two regimes separated by $N=1{,}000$. For a small $N$, the error coming from the use of finite samples dominates the estimate accuracy. When we consider a large data set, the error mainly comes from the model mispecification, because the stochastic simulator cannot be exactly represented by a GLaM. This phenomenon is not significant for the entropy estimation, which demonstrates a relatively consistent decay.

The accuracy of the oracle normal approximation is reported with red dash lines in Figure 2. The error shown is only due to model misspecifications (because the true response distribution is not Gaussian) since we use the underlying true mean and variance. For both measures, the medians of the errors of GLaMs built on $N=250$ model runs are smaller than those of the normal approximation. For $N\geq 1{,}000$, the GLaM clearly outperforms the oracle of Gaussian-type mean-variance models. This example illustrates the limits of such Gaussian-type models in practice.

Finally, the errors of GLaMs are below $0.05$ for $N\geq 1{,}000$ indicating that the surrogate is able to explain over 95% of the variance of the target functions.

For sensitivity analysis, we focus on the classical first-order and the entropy-based total Sobol’ indices. Figure 3 and Figure 4 show the convergence of GLaMs for estimating these quantities of each input variable. The reference values are derived from Equation 30. As shown by the two figures, this toy example is designed to have $X_{2}$ as the most important variable according to the classical first-order Sobol’ indices, which also indicates that $X_{2}$ contributes the most to the variance of the mean function $m(\bm{X})$. In contrast, it has zero effect to the entropy. In comparison, $X_{1}$ is the dominant variable for the variation of the entropy $h(\bm{X})$. Because $X_{3}$ mainly controls the shape of the response distribution (especially the right tail), it has a minor first-order effect to the mean function, which leads to a very small value of $S_{3}$. In contrast, the entropy depends on the distribution shape, and thus $S^{h}_{T_{3}}$ is not negligible. The results reveal that GLaMs capture this characteristic and yield accurate estimates for both classical Sobol’ indices and entropy-based Sobol’ indices.

Similar to Figure 2, we also reported the sensitivity indices calculated by using Gaussian approximations with the true mean and variance. Because the classical first-order indices depend only on the mean and variance functions, the oracle Gaussian model gives the exact values. Therefore, we showed only the results for the entropy-based Sobol’ indices in Figure 4. It is clear that the Gaussian approximation with the true mean and variance demonstrates a significant bias. In contrast, GLaMs show nearly no bias and approximate much more accurately the reference values.

In this example, we perform the global sensitivity analysis for a Heston model used in mathematical finance [37]. The Heston model describes the evolution of a stock price $Y_{t}$. It is an extension of the geometric Brownian motion by modeling the volatility as a stochastic process $v_{t}$, instead of considering it as constant. Hence, the Heston model is a stochastic volatility model and consists of two coupled stochastic differential equations:

with

where $W^{1}_{t}$ and $W^{2}_{t}$ are two Wiener processes with correlation coefficient $\rho$, which introduce the intrinsic randomness of the stochastic model. The model parameters $\bm{x}=(\mu,\kappa,\theta,\sigma,\rho,v_{0})$ are summarized in Table 1. The range of the last five input variables are selected based on the parameters calibrated from real data (S&P 500 and Eurostoxx 50) [38]. The range of the first variable $\mu$ is set to $[0,0.1]$ to take the uncertainty of the expected return rate into account. Without loss of generality, we set $Y_{0}=1$.

In this example, we are interested in the stock price after one year, i.e., $Y_{t}(\bm{x})$ with $t=1$. A closed form solution to Equation 31 is generally not available. To get samples of $Y_{1}(\bm{x})$, we simulate the entire time evolution of $Y_{t}$ and $v_{t}$ for a given $\bm{x}$ using Euler integration scheme with $\Delta t=0.001$ over the time interval $[0,1]$. Note that when simulating the bivariate process $(Y_{t},v_{t})$, a problem may happen: since $v_{t}$ follows a Cox–Ingersoll–Ross process [38], the simulation scheme can generate negative values for $v_{t}$. To overcome the problem, we apply the full truncation scheme, which replaces the update of $v_{t}$ by $\max\left(v_{t},0\right)$ [38].

Figure 5 shows two response PDFs predicted by a surrogate built upon $N=2{,}000$ model runs. The reference histograms are obtained from $10^{4}$ repeated model runs with the same input parameters. We observe that the variance of the response distribution is not constant, e.g., 0.065 and 0.027 for the two illustrated PDFs. Moreover, the PDF shape varies: it changes from symmetric to left-skewed distributions depending on the model parameters. This would be difficult to approximate with a simple distribution family such as normal or lognormal. In contrast, GLaMs are able to accurately capture this shape variation, because of the flexibility of generalized lambda distributions.

Even though a closed form distribution of $Y_{1}(\bm{x})$ does not exist, the mean function $m(\bm{x})={\mathbb{E}}\left[Y_{1}(\bm{x})\right]$ can be derived analytically:

As a result, we use $\varepsilon_{m}$ defined in Equation 29 with ${\rm QoI}(\bm{x})=m(\bm{x})$ to assess the convergence of the surrogate. In addition, we also consider the expected payoff of an European call option. The payoff $C(\bm{x})$ and the expected payoff $m_{C}(\bm{x})$ of an European call option are defined by

where $K$ is the strike price and set to $1$ in the following analysis. In finance, $m_{C}$ not only is important for making investment decisions but also has a very similar form to the option price [39]. For the Heston model, numerical methods based on the Fourier transform have been developed to calculate the expected payoff without the need for Monte Carlo simulations [37]. For the GLaM surrogate, this quantity can also be calculated numerically (see Section 6.1). As a second performance index, we compute the associated error denoted by $\varepsilon_{C}$ (Equation 29) for the convergence study.

Figure 6 shows box plots of the errors $\varepsilon_{m}$ and $\varepsilon_{c}$ for $N\in\{500;1{,}000;\allowbreak 2{,}000;\allowbreak 4{,}000;\allowbreak 8{,}000;16{,}000\}$. Both $\varepsilon_{m}$ and $\varepsilon_{c}$ are relatively large for $N\leq 2{,}000$. This is mainly due to the fact that the variability of the model response is dominated by the intrinsic randomness: the model parameters $\bm{X}$ altogether are only able to explain about 2% of the variance of the output (i.e., $S_{\left\{1,\ldots,6\right\}}=0.02$). In other words, the stochastic simulator has a very small signal-to-noise ratio ${\rm SNR}=0.02/(1-0.02)\approx 0.02$. Since GLDs are flexible, a few data scattered in a moderately high dimensional space may not provide enough information of the response distribution variation. We observe that for $N\leq 1{,}000$, the selection procedure proposed in Algorithm 1 can choose $\lambda^{\operatorname*{PC}}_{1}(\bm{x})$ and $\lambda^{\operatorname*{PC}}_{2}(\bm{x})$ being only constant. Such a model is too simple and thus fails to capture the variations of the scalar quantities. Consequently, it is necessary to have enough data to achieve an accurate estimate: when increasing the size of $N$ of the LHS design, we observe a clear decay of the errors.

We now study the convergence for the Sobol’ indices estimations. According to Equation 33, the mean function depends only on the first input variable $X_{1}$, which contributes little (2%) to the total variance of $Y_{1}(\bm{X})$. This implies that the classical Sobol’ indices are not informative (they are either 0 or very close to 0). However, we cannot ignore the variability of the input variables because the response distribution demonstrates a clear dependence on the input parameters, as shown in Figure 5. Therefore, we focus on the accuracy of the expected-payoff-based total Sobol’ indices, denoted by $S^{C}_{T_{{\bm{\mathsf{u}}}}}$.

As a second quantity of interest, we also calculate the total Sobol’ indices associated to the $95\%$-superquantile, referred to as $S^{\rm sq}_{T_{{\bm{\mathsf{u}}}}}$. Superquantiles are known as the conditional value-at-risk, which is an important risk measure in finance [40]. The $\alpha$-superquantile of a random variable $Y$ is defined by

where $q_{\alpha}$ is the $\alpha$-quantile of $Y$. This quantity corresponds to the conditional expectation of $Y$ being larger than its $\alpha$-quantile. For the Heston model, this quantity does not have an analytical closed form, whereas ${\rm sq}_{\alpha}$ of a GLD can be derived analytically (see Section 6.1).

We use $10^{5}$ Monte Carlo samples to evaluate (numerically) the function $m_{C}(\bm{x})$ to obtain a reference value for each $S^{C}_{T_{{\bm{\mathsf{u}}}}}$. To calculate the Sobol’ indices associated to the 95%-superquantile ${\rm sq}_{95}(\bm{X})$, it is necessary to evaluate the function ${\rm sq}_{95}(\bm{x})$. Because it cannot be analytically derived for the Heston model, we use $10^{4}$ replications to calculate the sample 95%-superquantile $\hat{{\rm sq}}_{95}(\bm{x})$ as an estimate for ${\rm sq}_{95}(\bm{x})$. Then, we treat it as a deterministic function and use $10^{4}$ samples to estimate each Sobol’ index. This indicates that a total number of $7\times 10^{8}$ model runs are performed to obtain the six reference 95%-superquantile-based total Sobol’ indices. Because only $10^{4}$ samples are used to estimate each $S^{\rm sq}_{T_{{\bm{\mathsf{u}}}}}$, we use bootstraps [41] to calculate the 95% confidence interval to account for the uncertainty of the Monte Carlo simulation.