{'fraction_non_alphanumeric': 0.03693983334426875, 'fraction_numerical': 0.012039892395512106, 'mean_word_length': 4.449231319270647, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Global sensitivity analysis of rare event probabilities

By their very nature, rare event probabilities are expensive to compute; they are also delicate to estimate as their value strongly depends on distributional assumptions on the model parameters. Hence, understanding the sensitivity of the computed rare event probabilities to the hyper-parameters that define the distribution law of the model parameters is crucial. We show that by (i) accelerating the calculation of rare event probabilities through subset simulation and (ii) approximating the resulting probabilities through a polynomial chaos expansion, the global sensitivity of such problems can be analyzed through a double-loop sampling approach. The resulting method is conceptually simple and computationally efficient; its performance is illustrated on a subsurface flow application and on an analytical example.


Introduction
Quantifying rare event probabilities is often needed when modeling under uncertainty [5,21,27,28]. Rare events are commonly associated with system failures or anomalies which pose a risk; it is thus imperative that rare event probabilities be computed reliably. For the sake of concreteness, we consider q to be a scalarvalued quantity of interest (QoI) whose inputs are drawn from the sample space Θ ⊆ R d with associated sigma algebra F and probability measure P. For a given thresholdτ , the corresponding rare event probability is defined as Pτ = P(q(θ) >τ ), where θ ∈ Θ is a random vector whose entries represent uncertain model parameters. Rare event probabilities are notoriously challenging to compute; indeed, basic Monte-Carlo simulations of (1) are inefficient in this context for the simple reason that few samples actually hit the rare event domain. Several methods have been proposed to compute Pτ more efficiently, ranging from importance sampling and Taylor series approximations to subset simulation, the latter of which we use in this article, see for instance [5] and Section 3. The evaluation of the rare event probability (1) requires the distribution law governing the model parameters θ. In practice, such a law is typically assumed. Clearly, Pτ depends on these assumptions; should they be misguided, the resulting rare event probability is likely be misleading. We let ξ denote a set of hyper-parameters charecterizing the distribution law of θ. It is crucial to understand the sensitivity of Pτ to ξ. In this article, we develop an efficient method to quantify, through global sensitivity analysis (GSA), the robustness of Pτ to the choice of hyper-parameters characterizing the distribution law of the model parameters.
To account for the uncertainty in ξ, we model the corresponding hyper-parameters as random variables. The rare event probability takes the form Pτ (ξ) = P({q(θ) >τ } | ξ). (2) A number of recent studies have considered how to assess the sensitivity of rare event estimation procedures to uncertain inputs and/or to the distributions of these inputs. There is a general consensus that the naive double-loop approach -whereby for each realization of ξ multiple samples of θ are used to estimate Pτ -is Draft infeasible but for the simplest of problems. An early work [21] combines rare event estimation techniques with the traditional Monte Carlo approach for GSA of the hyper-parameters. Several studies introduce new sensitivity measures [9,12,13,17] which are tailored to make the rare event SA process more tractable.
Others perform sensitivity analysis in the joint space of both input parameters and hyper-parameters [9,13,30,31]. These methods increase computational efficiency through use of local SA methods [9], surrogate models [13], kernel density estimates [31], and Kriging [30]. A thorough overview of current methods at the intersection of SA and rare event simulation can be found in [8] .
Our main contribution is to show that a double-loop approach can in fact be not only feasible but computationally expedient in order to perform GSA of Pτ (ξ) with respect to ξ. This may seem counterintuitive since, while informative, this type of second level sensitivity analysis is expensive. Our approach is however structurally simpler than most of the previously cited work and achieves computational efficiency through a combination of fast methods for rare event simulations together with the use of surrogate models. Specifically, we rely on subset simulation [5] to estimate rare event probabilities and approximate Pτ (ξ) using a polynomial chaos expansion (PCE), see respectively in Section 3 and Section 4. GSA is performed through a variance-based approach: crucially, the Sobol' indices for appropriate approximations to Pτ (ξ) can then be obtained "for free" through analytical formulae. To demonstrate the efficiency gains of the proposed method, we present an illustrative example in Section 2 and deploy our approach on it in Section 5.1. In Section 5.2, we apply the method to a Darcy flow problem requiring multiple estimates of the rare event probability to show feasibility in a more computationally demanding framework. We discuss additional challenges, perspectives and future work in Section 6.

A motivating example
We consider the following illustrative example [5,22,26] throughout the article where q is the QoI in (1) and θ = θ 1 · · · θ d with independent normally distributed entries θ i ∼ N (µ i , σ 2 i ), i = 1, . . . , d. It is elementary to check that, for any values of the hyper-parameters ξ = µ 1 . . . µ d σ 2 For a given ξ, the rare event probability is simply We model the uncertainty in the hyper-parameters by considering them as independent uniformly distributed random variables with a 10 percent perturbation around their respective nominal values. Figure 1 illustrates the case d = 5 with ξ nom = 1 2 3 4 5 10 8 6 4 2 as the nominal value for ξ. In particular, Figure 1, right, shows how the uncertainty in Pτ changes asτ varies. Asτ increases, i.e., as the event becomes rarer, the uncertainty in Pτ -measured through its coefficient of variation -increases. We contend that this latter behavior is generic for rare event simulations, establishing the need for methods allowing the quantification of the effects of hyper-parameter choices on the uncertainty in Pτ .
Note that as the event becomes less rare, µ will grow resulting in the diminishing of the bound on the coefficient of variation. We point out that the inequality (6) can be obtained directly from the more general Bhatia-Davis inequality [6].

Rare event simulation
Monte Carlo simulation is a standard way of approximating the rare event probability Pτ defined in (1).
In the case of rare events, i.e., of small probabilities Pτ , the basic MC estimator (8) becomes computationally inefficient. Indeed, consider the coefficient of variation δ(P M C τ ) of the above estimator and observe In other words, ensuring a given accuracy requires N ≈ 1 Pτ δ 2 . For increasingly rare events, i.e. small Pτ , the error in (8) will increase accordingly. Standard MC methods are thus poor candidates for rare event estimation.

Draft
The subset simulation method We rely on the subset simulation (SS) method [4,25] to accelerate rare event computation. This approach decomposes the rare event estimation problem into a series of "frequent event" estimation problems that are more tractable; it has been observed that this may reduce the coefficient of variation by more than an order of magnitude over standard MC [4,5,25]. This corresponds to a substantially lower computational burden for estimating rare event probabilities.
Consider the rare event domain F = {θ ∈ Θ | q(θ) >τ } and a sequence of nested subsets of F . . , L with τ 1 < τ 2 < · · · < τ L =τ . The rare event probability Pτ can thus be decomposed into a product of conditional probabilities with, by convention, F 0 = Θ. Computing Pτ according to (10) requires an efficient and accurate method for estimating the L conditional probabilities. We use a modification of the Metropolis-Hastings algorithm to accomplish this [4]. This modified Metropolis algorithm (MMA), which belongs to the family of Markov Chain Monte Carlo (MCMC) methods, draws samples from a conditional distribution and either accepts or rejects the samples based on a chosen acceptance parameter. One uses samples that belong to F i−1 as seeds for estimating the conditional probability P(F i |F i−1 ). We refer the interested reader to [5], which provides a high level discussion of MMA as well as other variants of SS; a more thorough analysis of MMA and MCMC algorithms can be found in [22]. Choosing a proper sequence of thresholds {τ i } L i=1 is a major challenge of the SS method. Since one has little prior knowledge of the PDF of q(θ), it is often not feasible to prescribe the sequence of thresholds a priori. Instead, one may require that P(F i | F i−1 ) = p 0 , i = 1, . . . , L − 1, for some chosen quantile probability p 0 [4]. We can then iteratively estimate the proper threshold at each "level" of the algorithm; the SS estimator of (10) takes the form where the final conditional probability P(F L | F L−1 ) is estimated via the MMA procedure mentioned earlier.
Although p 0 = 0.1 is a standard choice in engineering applications, there has been significant work done to determine optimal values for p 0 ; this, in general, depends on the QoI under consideration. It has been shown that, for practical purposes, the optimal p 0 lies in the interval [0.1, 0.3] and that, within this interval, the efficiency of SS is insensitive to the particular choice of p 0 [32]. With the approach for computing the sequence of thresholds in (10), each τ i is a random variable, estimated via a finite number of conditional samples. Consequently the number of levels or iterations necessary to terminate SS is also random. For a sufficiently large number of samples, the number of levels is given in [26] as Implementation of subset simulation For completeness, we provide an algorithm outline for the SS method in Algorithm 1. We assume Gaussian inputs in the examples considered in this article; the SS method can however be applied to non-Gaussian input distributions, see Appendix B of [20] for details. Additional information on the implementation of the SS algorithm, including the MMA implementation, is for instance available in [4,26,28]. As this MCMC implementation reuses the input parameters from each previous level to estimate the threshold for the next level, this method does not require any burn-in samples to draw from the conditional distribution; it begins by sampling from the previous rare event domain. On the theoretical side, the SS algorithm is asymptotically unbiased andP SS τ converges almost surely to the true rare event probability, Pτ . For a detailed convergence analysis of SS and derivation of its statistical properties, see [4].

Algorithm 1 Subset Simulation
Input: Rare event threshold,τ , samples per level, NSS, quantile probability, p0, routine that evaluates QoI, q(θ) Output: Estimate of rare event probability:P SS τ 1: Draw NSS samples of θ from appropriate distribution 2: Compute NSS samples of the QoI, compute τ1 as the p0 quantile 3: Save the NSS · p0 inputs such that q(θ) > τ1 as seeds for the next level 4: i ← 1 {i indicates the current level} 5: while τi <τ do 6: i ← i + 1 Using MCMC seeds, evaluate the QoI and compute τi as the p0 quantile 9: Save the NSS · p0 inputs such that q(θ) > τi as seeds for the next level 10: end while 11: Using seeds from FL−1, sample the QoI and estimate P(FL | FL−1) using MC 12: Computational cost We turn now to the computational cost of estimating Pτ using SS. The computational cost is measured in terms of the number of function evaluations required to run the algorithm. As the number of levels L is random, so is the computational cost associated with SS. For simplicity, we assume for our cost analysis that a sufficient number of samples has been used so that L does not vary. The total number of QoI evaluations required by SS is L · N SS , where N SS is a user-defined parameter that determines the number of samples per intermediate level of the iteration. Say, for example, the true rare event probability is 10 −6 and we wish to estimate Pτ with a coefficient of variation within δ = 0.1. For standard MC sampling, we would need N ≥ 1/(δ 2 · Pτ ) = 10 8 samples of the QoI. Take the SS method with a quantile probability of p 0 = 0.1. Then, according to (12), we would have L = 7, corresponding to 7 levels of conditional probabilities. The coefficient of variation for each of the conditional probabilities is more difficult to quantify, however, as in the case of the standard MC estimator, they are proportional to 1/p 0 ; see [4]. In this case, one would expect to see a significant reduction in the cost of estimating Pτ with SS.
We lastly emphasize the power of SS for estimating rare event probabilities in the context of QoIs with high-dimensional inputs. Not only does SS improve upon the slow convergence rates of standard MC by a wide margin, it also inherits the property of having a convergence rate independent of input dimension.

Surrogates for GSA of rare event probabilities
We seek to apply variance-based GSA to Pτ (ξ), defined in (2), with respect to components of ξ. To mitigate the computational expense of performing such analysis, we combine the SS algorithm and surrogate models, in the form of polynomial chaos expansions (PCEs). We assume ξ to be an M -dimensional vector with independent entries. The procedure, which amounts to double-loop sampling, is outlined below: • For each j ∈ {1, . . . , N samp }, estimate Pτ (ξ (j) ) using SS; denote these estimates by P • Compute the Sobol' indices ofPτ (ξ).
Instead of using SS for computing Pτ (ξ (j) ), one may be tempted to apply a surrogate further "upstream" by computing a surrogate modelq ξ (j) (θ) for q(θ) from samples {q(θ (k) )} n k=1 drawn from law of θ as determined by ξ (j) . This surrogate model of q can then be used for fast approximation of the rare event probability Pτ (ξ (j) ). This procedure, however, has two major pitfalls: (i) an expensive surrogate modeling procedure must be carried out for each j ∈ {1, . . . , N samp } and, more importantly, (ii) surrogate models are typically poorly suited to the task of rare event estimation. Indeed, surrogates typically fail to capture the tail behavior of the distribution of the QoI q, making them unsuitable for rare event simulations. This shortcoming Draft is well-documented in the uncertainty quantification literature [18,23] although efforts are being made to tailor the surrogate model construction process for the efficient estimation of rare event probabilities [18,19].
PCE surrogate for rare event probability Our approach leverages the properties of PCE surrogates for fast estimation of Sobol' indices [11,16]; it also takes advantage of the regularity of the mapping ξ → Pτ (ξ). Specifically, assuming the PDF of ξ satisfies certain (mild) differentiability and integrability conditions, one can show that Pτ (ξ) is a differentiable function of ξ; see [3,Proposition 3.5].
The PCE of Pτ (ξ) is defined asPτ where Ψ 0 , · · · , Ψ N PC belong to a family of orthogonal polynomials and β 0 , · · · , β N PC are the (scalar) PCE coefficients. The specific family of polynomials is chosen to guarantee orthogonality with respect to the PDF of ξ; see, e.g., [16]. We use a total order truncation scheme for the PCE: the multivariate polynomial basis contains all possible polynomial basis elements up to a total polynomial order r. In this case, N PC in (13) satisfies The coefficients β 0 , . . . , β N PC can be computed in a number of ways, including non-intrusive spectral projection or regression [7,11,16]. A regression based approach is preferred here because the evaluations of Pτ are noisy due to sampling errors incurred in the SS procedure. We estimate the vector β = [β 0 , β 1 , · · · , β N PC ] from function evaluations P (j) τ = SS(Pτ (ξ (j) )), j = 1, . . . , N samp , by solving the penalized least squares problem In (14), the penalty parameter λ acts as a sparsity control on the recovered PCE coefficients. We generate the realizations {ξ (j) } Nsamp j=1 of the hyper-parameter vector through Latin hypercube sampling; for further details on the implementation of sparse regression for PCE, see [14,15]. The numerical results in Section 5 are obtained using the SPGL1 solver [29].
GSA of Pτ using the PCE surrogate As is well-known, the Sobol' indices of a PCE surrogate can be computed analytically. For example, the first order Sobol' indices, S i (Pτ ), i = 1, . . . , M , of Pτ can be approximated as follows: where K i denotes the set of all PCE terms that depend only on ξ i . Sobol' indices for arbitrary subsets of variables, as well as total indices, can be computed in an analogous manner [2,16]. In practice, PCE surrogates with modest accuracy are often sufficient to obtain reliable estimates of Sobol' indices. While the above approach for GSA of Pτ does require repeated simulations of the QoI q during the calls to the SS algorithm, it still provides orders of magnitude speedup over the standard "pick and freeze" MC methods, also known as Saltelli sampling, for computing the Sobol' indices of Pτ [24]. Indeed, a fixed sample {ξ (j) } Nsamp j=1 with modest N samp is sufficient to compute the PCE surrogate from which the Sobol' indices can be computed at a negligible computational cost. Moreover, the sparse regression approach for estimating PCE coefficients is forgiving of noisy function evaluations. Therefore, large sample sizes are not needed in the calls to the SS algorithm. We demonstrate the merits of the proposed approach in our computational results presented in Section 5.

Numerical results
We summarize, in Section 5.1, the computational results for the motivating example from Section 2; a more challenging model problem involving flow through porous media is considered in Section 5.2.

Results for the analytic test problem
We consider the example from Section 5.1 and study Pτ withτ = 3. To establish a baseline for the values of the Sobol' indices of Pτ (ξ), we compute the total order Sobol' indices directly from (5) using Saltelli sampling. The reference Sobol' indices are computed with 10 6 samples for each of the conditional terms; convergence was numerically verified. We plot the reference total indices in Figure 2 for comparison. We now compare the reference indices with those obtained through the PCE surrogate when Pτ (ξ) is computed analytically using Equation (5). We allocate 10 3 samples of Pτ (ξ) each for the Saltelli sampling method and sparse regression PCE method. The Saltelli method requires N (d + 1) samples [24] and so we divide the budget of 10 3 samples equally among each conditional term. Each PCE coefficient can be estimated using the full set of 10 3 samples. For a fair comparison, we use Latin hypercube sampling for both the PCE and Saltelli method. We also use a total PCE order of 3 and the penalty parameter λ = 5 × 10 −2 . Given that the set of total indices is computed, in each method, using a finite number of samples, each index is a random variable with an associated distribution. We compare the standard deviation of each total index for the two GSA methods. In each case, we compute 10 3 realizations of the full set of total indices and compare their respective standard deviations in Figure 2.   Figure 2 illustrates the higher accuracy, or lower variance, of PCE with sparse regression over Saltelli sampling: the standard deviation of the largest Sobol' index is roughly 3 times smaller with sparse regression than it is with Saltelli sampling. This gap in accuracy appears to diminish for smaller indices, although the methods do not show comparable accuracy until the indices are below 0.1. As Pτ can be expressed analytically, there may be additional benefits of the sparse regression method to be seen when one considers performing GSA on a rare event probability with noise due to sampling. We note that the total order of the PCE basis and the penalty parameter λ, which are user-defined parameters, can be changed without the need for additional runs of SS. These parameters can be cross validated in a post-processing step after the rare event simulation step, providing flexibility in this approach without adding any significant computational burden.
When combining PCE-based GSA with SS for estimating Pτ (ξ), there is a tradeoff between the inner loop cost of estimating Pτ via SS and the outer loop of aggregating Pτ samples to build the PCE. In Figure  3 shows that for N SAMP = 10 3 , we only need a modest N SS to approximate the Sobol' indices. Indeed, for N SS = 500, we are able to resolve the total indices very well. We also examine the case of N SAMP = 10 4 in Figure 3 (bottom). Again, we are able to resolve the total indices well using only N SS = 500 and are able to achieve the correct ordering for as little as N SS = 100.
These results indicate that (i) a modest number of samples allocated to SS is enough to get a rough estimate of Pτ and (ii) a moderate number of realizations of Pτ (ξ) is then sufficient for accurate GSA. In other words, given rather poor estimations of Pτ , we are still able to extract accurate GSA results, due to the fact that the sparse regression technique is robust to noisy QoI evaluations.

Subsurface flow application
We consider the equations for single-phase, steady state flow in a square domain D = [0, 1] 2 : where κ is the permeability, µ is the viscosity, and p is the pressure. The boundaries Γ 1 , Γ 2 , and Γ 3 indicate the left boundary, the right boundary, and the top/bottom boundaries, respectively. The Darcy velocity is defined as v = − κ µ ∇p. In the present study, we let µ = 1. The source of uncertainty in this problem is in the permeability field, which we model as a random field. We consider the flow of particles through the medium and focus on determining the probability that said particles do not reach the outflow boundary in a given amount of time. This problem has been used previously as a test problem for rare event estimation in [28] as it pertains to the long-term reliability of nuclear waste repositories. Our goal is to perform GSA with respect to the hyper-parameters that define the distribution law of the permeability field.
The statistical model for the permeability field Following standard practice [10,28], we model the permeability field as a log-Gaussian random field: where x ∈ D and ω belongs to sample space that carries the random process. Here,ā is the mean of the random field, σ a is a scalar which controls the pointwise variance of the field, and z is a centered (zero-mean) random process. We let the covariance function of z be given by where x and y denote the correlation lengths in horizontal and vertical directions. The random field is represented via a truncated Karhunan-Lòeve expansion (KLE): In this representation, θ 1 , . . . , θ N KL are independent standard normal random variables and (λ i , e k ), k = 1, . . . , N KL , are the leading eigenpairs of the covariance operator of the process. Our setup for the uncertain log-permeability field follows the one in [10]: we use permeability data from the Society for Petroleum Engineers [1] to defineā. Once we truncate the KLE, the random vector θ = [θ 1 θ 2 · · · θ N KL ] fully describes the uncertainty in the log-permeability field.
To ensure that the KLE accurately models the variability of the infinite-dimensional field, we examine the eigenvalue decay of the covariance operator with the goal to truncate the KLE so that at least 90% of the average variance of the field is maintained. For x = y = 0.4, which are the smallest correlation lengths considered in the present study, we require at least N KL = 126. The number of retained KL modes then determines the dimensionality of the rare event estimation problem, and is henceforth fixed at 126. The dimension independent properties of SS is advantageous in this regime.
For illustration, we plot two realizations of the random field, with the corresponding pressure and velocity fields obtained by solving the governing PDE (16), in Figure 4. In our computations, we solve the PDE using piecewise linear finite elements in Matlab's finite element toolbox with 50 mesh points in each direction.
The QoI and rare events under study The position x of a particle moving with the flow through the medium is determined by the following ODE dx dt = v, where v is the Darcy velocity. We consider a single particle with initial position at x 0 = 0 0.5 . The solution x of (20) depends not only on time but also on θ due to dependence of v on θ, i.e., x = x(t, θ). We take the QoI q as the hitting time, i.e., the time it takes a particle to travel through the medium from left to right We aim to determine the rare event probability Pτ = P(q >τ ). The parameters x , y , and σ a parametrize the uncertainty in the permeability field; we consider them as hyper-parameters and set ξ = [ x y σ a ] . We set the nominal values of the hyper-parameters ξ nom = [0.4 0.4 0.8] . We simulate realizations of the permeability field at these nominal hyper-parameters and plot the distribution of q. Each of these realizations requires one PDE solve and one ODE solve. As illustrated in Figure 5, the distribution for q corresponds to a heavy-tailed distribution. We select as the thresholdτ = 4.5 and consider quantifying the sensitivity of Pτ (ξ) with respect to the hyper-parameters defining the KLE. Our hyper-parameters are drawn from a uniform distributed centered at ξ nom with a spread of plus or minus 10% of ξ nom . We use these SS estimations of Pτ (ξ) in order to build the corresponding PCE surrogate, where the polynomial basis is truncated at a total polynomial order of 5. Note the decision of where to truncate the PCE basis does not need to be made prior to estimating the set of Pτ (ξ) samples. The samples for the hyper-parameters are drawn using a Latin hypercube sampling scheme. We use 10 3 estimations of Pτ (ξ) to construct the PCE surrogate. Again, we use sparse regression to recover the PCE coefficients, while promoting sparsity in the set of PCE coefficients, and so mitigating the effects of noise induced by SS. In Figure 6, we use two different values of λ when promoting sparsity in order to illustrate the effect of λ on the results. Note that when λ is made smaller, the PCE coefficients decrease in magnitude, promoting a sparser PCE spectrum. In both cases, the ordering of the total Sobol' indices remains consistent, Draft  and thus, conclusions with respect to parameter sensitivity are unaffected. For this experiment, we therefore conclude that choosing λ by trial and error is sufficient. Should one encounter a scenario where the GSA results are more sensitive to λ, more systematic approaches are possible [7,15].
We lastly return to the key point made in Section 5.1, that the proposed method is capable of producing reliable GSA results, while using a modest number of inner and outer loop samples (N SS and N samp , respectively). In Figure 7, we report results corresponding to N SS = 500. In the left panel of the Figure we study the effect of N samp on the PDF of the PCE surrogate. In the right panel, we plot the Sobol' indices corresponding to each of the computed surrogates. The results in Figure 7 should also be compared with those in Figure 6, where larger values of N SS and N samp were used. This experiment indicates that Pτ and the Sobol' indices themselves can be well-approximated with a modest number of samples in both the inner and outer loops. In this case, using both N SS and N samp on the order of 10 2 is sufficient for obtaining accurate GSA results. The combined cost of this method is thus reduced by a significant margin compared with the similar results in Figure 6. The efficiency gains of this method indicate the potential for deployment on problems which would otherwise be intractable. Figure 7: Distributions of Pτ for N SS = 500 and varying values of N samp . For each N samp , we build the PCE surrogate and approximate its PDF with 10 5 samples. The set of Pτ samples used for differing N samp is nested within sets of larger samples. Corresponding total indices are included, computed directly from the PCE surrogates.

Conclusion and future work
We have shown that the feasibility of the standard double-loop approach for GSA of rare event probabilities can be significantly extended beyond simple applications. This requires appropriate acceleration methods; in our case, this is achieved through subset simulation and the choice of a surrogate model allowing for the analytical calculation of Sobol' indices. This approach is conceptually simple and does not require the development of new, ad hoc sensitivity concepts. While we have extended the range of applicability of the double-loop approach, we acknowledge that more research is needed to deal with computationally expensive, high-dimensional problems.
The efficiency of our method crucially depends on working with surrogate models for which sensitivity measures -here, Sobol' indices -can be computed cheaply or "for free"; this clearly and strongly limits the type of GSA which can be carried out by the approach. More generally, if q is the original QoI and ifq is the resulting QoI for a given surrogate model, more work is needed to understand the relationship between the approximation error q −q and the resulting GSA error S(q) − S(q) where S(·) is some sensitivity measure; more explicitly, there may be room for the development "cheap" surrogate models with moderate approximation errors and small GSA errors. Additionally, both our sensitivity analysis method as well as surrogate modeling approach rely on the assumption that the hyper-parameters are independent. In some cases one might be interested in GSA of rare event probabilities to both hyper-parameters and additional parameters in a model that might be uncertain and possibly correlated. Therefore, another interesting line of inquiry is to consider GSA of rare event probability with respect to correlated parameters. Further study may also include extensions of our approach to other moment-based QoIs (e.g. CDF approximation, skewness, kurtosis) and the use of perturbation-based methods for GSA [17] as opposed to considering a discrete set of hyper-parameters.