Density modification based reliability sensitivity analysis
Paul Lemaître, Ekatarina Sergienko, Aurélie Arnaud, Nicolas Bousquet, Fabrice Gamboa, Bertrand Iooss
aa r X i v : . [ m a t h . S T ] M a r Density modification based reliability sensitivity analysis
P. Lemaˆıtre a b ∗∗ E. Sergienko c d
A. Arnaud a N. Bousquet a F. Gamboa d B. Iooss a d a EDF R&D, 6 Quai Watier - 78401 Chatou ; b INRIA Sud-Ouest, 351 cours de la lib´eration - 33405 Talence ; c IFP EN, 1 avenue de Bois-Pr´eau - 92852 Rueil-Malmaison ; d Institut de Math´ematiques de Toulouse, 118 route de Narbonne - 31062 Toulouse
August 14, 2018
Abstract
Sensitivity analysis of a numerical model, for instance simulating physical phenomena, is useful toquantify the influence of the inputs on the model responses. This paper proposes a new sensitivity index,based upon the modification of the probability density function (pdf) of the random inputs, when thequantity of interest is a failure probability (probability that a model output exceeds a given threshold).An input is considered influential if the input pdf modification leads to a broad change in the failureprobability. These sensitivity indices can be computed using the sole set of simulations that has alreadybeen used to estimate the failure probability, thus limiting the number of calls to the numerical model.In the case of a Monte Carlo sample, asymptotical properties of the indices are derived. Based onKullback-Leibler divergence, several types of input perturbations are introduced. The relevance of thisnew sensitivity analysis method is analysed through three case studies.
In the context of structural reliability, computer models are used in order to assess the safety of industrialsystems relying on complex physical phenomena. For instance, an electric operator would like to predictthe level of a potential river flood in order to determine the height of a dyke preventing any disaster. Inthis example, the computer model (simulating the hydraulic model) has some uncertain input variables (flowrate, river length, water height, etc.), that are modelled by random variables. In this paper, the computercode is a ”black-box” deterministic numerical model and the study focuses on one of its output. Due to therandomness of the model inputs, this output is a random variable more or less sensitive to the uncertaintyof the input variables.Sensitivity analysis (SA) is a tool used to explore, understand and (partially) validate numerical models.It aims at explaining the outputs regarding the input uncertainties ([14]). We use the “global SA” definitiongiven by Saltelli et al. [16] wherein the whole variation range of the inputs is considered. The application ofsuch an approach can be model simplification (by removing irrelevant modelling elements), input variablesranking or research prioritization. There is a wide range of SA techniques, regarding what type of problemthe experimenter faces with ([10]). For instance, screening methods are to be applied when there is alarge number of inputs, and few models assumptions. From a quantitative point of view, the most populartechniques are variance-based methods and the so-called Sobol’ indices ([16, 17]). These are based uponHoeffding decomposition of L function and functional variance decomposition [1].It should be noticed that most SA methods focus on real-valued continuous numerical output variables.When the output is a binary value (e.g. when the numerical model returns “faulty system” or “safe system”), ∗∗ Corresponding author. Email: [email protected]
1A techniques are underdeveloped.Some basic techniques can be quoted, such as Monte-Carlo filtering ([16])which consists in measuring differences between a “safe” sample and a “faulty” sample via standard statisticaltests. In structural reliability analysis, some sensitivity factors resulting from the First or Second OrderReliability Methods (FORM/SORM, [11]) can also be used to classify the impact of the inputs on the failureprobability. More recent works give methods combining the two objectives: estimating a failure probabilityand assessing the influence of the input uncertainty on this probability ([12, 13]).In this paper, a real-valued numerical model denoted by G : R d → R is considered. This model mayfurther be called the “failure function”. In practice, each run of G can be CPU time consuming. We areinterested in the (rare) event G ( X ) < G ( X ) ≥ X = ( X , ..., X d ) T is a d -dimensional continuous random variable whose joint probability densityfunction (pdf) is denoted f . For i = 1 , · · · , d , let f i denotes the distribution of X i (the marginal pdf). Wemake the assumption that all components of X are independent. The quantity of interest is the system failureprobability: P = Z { G ( x ) < } f ( x ) d x . The aim of this work is the quantification of the influence of each variable X i on this probability.Let us ask the question: what are the engineers motivations when he perform a SA on his/her black-boxmodel that produces a binary response? We provided an overview of the general objectives of SA: variableranking, model simplification, model understanding. But from our discussions with practitioners, we haveidentified three ”engineer motivations”: • the practitioner wants to determine which are the inputs that impact the most the failure event theinputs distributions being set and supposed to be perfectly known. This amounts to an absolute rankingobjective. • P will be impacted by the choice of the input distributions; the engineer wants to assess the influenceof this choice on the FP. Therefore the objective here is to quantify the sensitivity of the model outputto the family or shape of the inputs. • In practice, input distributions are estimated from data, thus leading to uncertainty on the values of thedistribution parameters. The practitioner wants to assess the influence of the distribution parameterson P . Therefore the objective here is to test the sensitivity of the model to the parameters of theinputs.In most studies, sensitivity indices for failure probabilities are defined in strong correspondence with agiven method of estimation (e.g. [11, 13]). Their interpretation is consequently limited. We propose in thisarticle to define new generic sensitivity indices. Our sensitivity index is based upon density modification, andis adapted to failure probabilities. A methodology to estimate such indices is derived. For simplicity reasons,a classical Monte Carlo framework is considered in the following. Additionally, the sensitivity index can becomputed using the sole set of simulations that has already been used to estimate the failure probability P ,thus limiting the number of calls to the numerical model.The outline of the article is the following: first we define a generic strategy of input perturbation inSection 2, based upon maximum entropy rules. We then present our index and its theoretical properties inSection 3, altogether with the estimation methodology. The behaviour of the indices is examined in Sec-tion 4 through numerical simulations in various complexity settings, involving toy examples and a realisticcase-study. Comparisons with two reference sensitivity analysis methods (FORM indices and Sobol’ indices)highlight the relevance of the new indices in most situations. The main advantages and remaining issues arefinally discussed in the last section of the article, that introduces avenues for future research. Our sensitivity analysis method requires to define a perturbation for each input. In general, and especiallyin preliminary reliability studies, there is no prior rule allowing to elicit a specialized perturbation for eachinput variable. We thus would like to propose a simple perturbation methodology, allowing the practitioner2o answer the questions itemized in the Introduction. Furthermore, we make the implicit hypothesis that theextreme values of the inputs lead to the (rare) failure event.Given a unidimensional input variable X i with pdf f i , let us call X iδ ∼ f iδ the corresponding perturbedrandom input. This perturbed input takes the place of the real random input X i , in a sense of modellingerror : what if the correct input were X iδ instead of X i ?More precisely, we suggest to define a perturbed input density f iδ as the closest distribution to the original f i in the entropic sense and under some constraints of perturbation. Information-theoretical arguments ([6])led us to choose the Kullback-Leibler (KL) divergence between f iδ and f i as a measure of the discrepancy tominimize under those constraints. Given the hypotheses and the needs, we focus on linear constraints thatcan be interpreted in term of moments perturbation. This will lead to a quantification of the impact on P ofeach variable. This will also provide results on the sensitivity of P to the choice of input distributions. Wewill later present the perturbations corresponding to a mean shift and a variance shift.Recall that between two pdf p and q we have KL ( p, q ) = Z + ∞−∞ p ( y ) log p ( y ) q ( y ) dy if log p ( y ) q ( y ) ∈ L ( p ( y ) dy ) . (1)Let i = 1 , · · · , d , the constraints are expressed as follows in function of the modified density f mod : Z g k ( x i ) f mod ( x i ) dx i = δ k,i ( k = 1 · · · K ) . (2)Here, for k = 1 , · · · , K , g k are given functions and δ k,i are given real. These quantities will lead to aperturbation of the original density. The modified density f iδ considered in our work is: f iδ = argmin f mod | ( ) holds KL ( f mod , f i ) (3)and the result takes an explicit form ([7]) given in the following proposition. Proposition 2.1
Let us define, for λ = ( λ , · · · , λ K ) T ∈ R K , ψ i ( λ ) = log Z f i ( x ) exp " K X k =1 λ k g k ( x ) dx , (4) where the last integral can be finite or infinite (in this last case ψ i ( λ ) = + ∞ ). Further, set Dom ψ i = { λ ∈ R K | ψ i ( λ ) < + ∞} . Assume that there exists at least one pdf f m satisfying (2) and that Dom ψ i is an openset. Then, there exists a unique λ ∗ such that the solution of the minimisation problem (3) is f iδ ( x i ) = f i ( x i ) exp " K X k =1 λ ∗ k g k ( x i ) − ψ i ( λ ∗ ) . (5)The theoretical technique to compute λ is provided in appendix A. Here are presented two kinds ofperturbations used further on. Mean shifting
The first moment is often used to parameterize a distribution. Thus the first perturbationpresented here is a mean shift, that is expressed with a single constraint: Z x i f mod ( x i ) dx i = δ i . (6)In term of SA, this perturbation should be used when the user wants to understand the sensitivity of theinputs to a mean shift - that is to say “what if the mean of input X i were δ i instead of E [ X i ]?”. Proposition 2.2
Considering the constraint (6), under the assumptions of Proposition 2.1 the expressionof the optimal perturbed density is f iδ i ( x i ) = exp( λ ∗ x i − ψ i ( λ ∗ )) f i ( x i ) (7) where λ ∗ is such that equation (6) holds. ψ i ( λ ) = log Z f i ( x i ) exp( λx i ) dx i = log ( M X i ( λ )) (8)where M X i ( u ) is the moment generating function (mgf) of the i − th input. With this notation, λ ∗ is suchthat Z x i exp ( λ ∗ x i − log ( M X i ( λ ∗ ))) f i ( x i ) dx i = δ i , which leads to Z x i exp ( λ ∗ x i ) f i ( x i ) dx = δ i M X i ( λ ∗ ) . This can be simplified to: M ′ X i ( λ ∗ ) M X i ( λ ∗ ) = δ i . (9)This equation may be easy to solve when the expression of the mgf of the input X i and of its derivative isknown. Variance shifting
In some cases, the expectation of an input may not be the main source of uncertainty.One might be interested in perturbing its second moment. This case may be treated considering a couple ofconstraints. The perturbation presented is a variance shift, therefore the set of constraints is: (R x i f mod ( x i ) dx i = E [ X i ] , R x i f mod ( x i ) dx i = V per ,i + E [ X i ] . (10)The perturbed distribution has the same expectation E [ X i ] as the original one and a perturbed variance V per ,i = Var [ X i ] ± δ i . Proposition 2.3
Under the assumptions of Proposition 2.1, for the constraint (10), the expression of theoptimal perturbed density is: f iδ i ( x i ) = exp( λ ∗ x + λ ∗ x − ψ i ( λ ∗ )) f i ( x i ) where λ ∗ and λ ∗ are so that equation (10) holds. As an example, the two kind of perturbations previously presented are provided for two families of inputs(Gaussian and Uniform) in figure 1. The perturbations are respectively a mean and variance increasing. Itis noticeable (and will be proved further on) that the shape is conserved for the Gaussian distribution whenshifting the mean or the variance. On the other hand, when increasing its mean, the Uniform distribution ispacked down on the right-hand boundary of its support. When increasing its variance, the density is packeddown on both boundaries of its support.
Perturbation of Natural Exponential Family
In general, when perturbating the input densities, theshape is not conserved. However in the specific case of Natural Exponential Family (NEF), the followingproposition can be derived.
Proposition 2.4
Assume that the original random variable X i belongs to the NEF, i.e. its pdf can be writtenas: f i,θ ( x i ) = b ( x i ) exp [ x i θ − η ( θ )] where θ is a parameter from a parametric space Θ , b ( . ) is a function that depends only of x i and η ( θ ) = log Z b ( x ) exp [ x i θ ] dx i is the cumulant distribution function. Considering the assumptions of Proposition 2.1, the optimal pdfsproposed respectively in Proposition 2.2 and Proposition 2.3 are also distributed according to a NEF. The proof comes from theorem 3.1 in [7]. The details of computation are given for a mean shift and avariance shift in Appendix B. 4igure 1: Mean shifting (left) and variance shifting (right) for Gaussian (upper) and Uniform (lower) distri-butions. The original distribution is plotted in solid line, the perturbed one is plotted in dashed line.5
Definition, estimation and properties of a sensitivity index
Given a unidimensional input variable X i with pdf f i and the corresponding perturbed random input X iδ ∼ f iδ . The perturbed failure probability becomes: P iδ = Z { G ( x ) < } f iδ ( x i ) f i ( x i ) f ( x ) d x (11)where x i is the i th component of the vector x . Independently of the mechanism chosen for the perturbation(see previous section for proposals), a good sensitivity index S iδ should have intuitive features that makeit appealing to reliability engineers and decision-makers. We believe that the following definition can fulfilthese requirements. Definition 3.1
Define the Density Modification Based Reliability Sensitivity Indices (DMBRSI) the quantity S iδ : S iδ = (cid:20) P iδ P − (cid:21) { P iδ ≥ P } + (cid:20) − PP iδ (cid:21) { P iδ
P i.e. if the remaining ( epistemic ) uncertainty on the modelling X i ∼ f i can increase the failure risk. In this case, the uncertainty on the concerned variable should be moreaccurately analysed. Conversely, if if P iδ < P , P can be interpreted as a conservative assessment of thefailure probability, with respect to variations of X i . In such a case, deeper modelling studies on X i appearless essential.Thirdly, given its sign, the absolute value of S iδ has simple interpretation and provides a level of theconservatism or non-conservatism induced by the perturbation. A value of α > P iδ = (1 + α ) P . If S iδ = − α < P iδ = (1 / (1 + | α | )) P .The postulated ability of S iδ to enlighten the sensitivity of P to input perturbations must be testedin concrete cases, when an estimator ˆ P N of P can be computed using an already available design of N numerical experiments. In this paper, N is assumed to be large enough such that statistical estimationstands within the framework of asymptotic theory. Besides, we assume for simplicity a standard Monte Carlodesign of experiments, according to which ˆ P N = P Nn =1 { G ( x n ) < } /N where the x , · · · , x N are independentrealisations of X . The strong Law of Large Numbers (LLN) and the Central Limit Theorem (CLT) ensurethat for almost all realisations ˆ P N −−−−→ N →∞ P and p N/ [ P (1 − P )]( ˆ P N − P ) L −−−−→ N →∞ N (0 , . The Monte Carlo framework allows P iδ to be consistently estimated without new calls to G , through a“reverse” importance sampling mechanism:ˆ P iδN = 1 N N X n =1 { G ( x n ) < } f iδ ( x ni ) f i ( x ni ) . This property holds in the more general case when P is originally estimated by importance sampling ratherthan simple Monte Carlo, which is more appealing when G is time-consuming [3, 8]. This generalization isdiscussed further in the text (Section 5). The following lemma ensures the asymtotic behaviour of such anestimator. Lemma 3.1
Assume the usual conditions(i) Supp ( f iδ ) ⊆ Supp ( f i ) , ii) Z Supp ( f i ) f iδ ( x ) f i ( x ) dx < ∞ ,then ˆ P iδN −−−−→ N →∞ P iδ and √ Nσ − iδN (cid:16) ˆ P iδN − P iδ (cid:17) L −−−−→ N →∞ N (0 , . The exact expression of σ − iδN is given inAppendix C, equation (15). It can be consistently estimated by ˆ σ iδN = 1 N N X n =1 { G ( x n ) < } (cid:18) f iδ ( x ni ) f i ( x ni ) (cid:19) − ˆ P iδN . The proof of this Lemma is given in Appendix C.1.
The asymptotic properties of any estimator of S iδ will depend on the correlation between ˆ P N and ˆ P iδN . Thenext proposition summarizes the features of the joint asymptotic distribution of both estimators. Proposition 3.1
Under assumptions (i) and (ii) of Lemma 3.1, √ N (cid:20)(cid:18) ˆ P N ˆ P iδN (cid:19) − (cid:18) PP iδ (cid:19)(cid:21) L −−−−→ N →∞ N (0 , Σ iδ ) where Σ iδ is given in Appendix C, Equation (16) and can be consistently estimated by ˆΣ iδ = (cid:18) ˆ P N (1 − ˆ P N ) ˆ P iδN (1 − ˆ P N )ˆ P iδN (1 − ˆ P N ) ˆ σ iδN (cid:19) . The proof of this Proposition is given in Appendix C.2.
Given ( ˆ P N , ˆ P iδN ), the plugging estimator for S iδ isˆ S i δN = " ˆ P iδN ˆ P N − { ˆ P iδN ≥ ˆ P N } + " − ˆ P N ˆ P iδN { ˆ P iδN < ˆ P N } . (12)In corollary of Proposition 3.1, applying the continuous-mapping theorem to the function s ( x, y ) = (cid:2) yx − (cid:3) y ≥ x + h − xy i y Assume that assumptions (i) and (ii) of Lemma 3.1 hold and further that P = P iδ , wehave √ N h ˆ S i δN − S i δ i L −−−−→ N →∞ N (cid:0) , d T Σ d (cid:1) (13) with d = (cid:18) ∂s∂x ( P, P iδ ) , ∂s∂y ( P, P iδ ) (cid:19) T for x = y , and ∂s∂x ( x, y ) = − y { y ≥ x } /x − y { y In this section, the methodology presented throughout the article is tested on two academic cases and a morerealistic industrial model. The new indices are compared to the results provided by two reference methods,FORM indices (or importance factors) and Sobol’ indices. First, we briefly present these two methods. The First Order Reliability Method (FORM)[11] is an estimation technique for a failure probability basedupon approximating the failure domain and solving an optimization problem. In practice, it is considered astandard solution in structural reliability, mainly because of its low cost. This method proceeds in four steps: • transformation of the input variables space into the standard space (a space where all the variables arestandard independent Gaussian); • search of the closest failure point to the origin of the standard space (also referred as the design point)via an optimisation algorithm; • approximation of the failure surface by an hyperplane; • failure probability estimation based on the geometry of the failure domain.Importance factors are byproducts of this method. These sensitivity measures aims at quantifying theimportance of a variable on the failure probability. From the design point u ∗ one writes: u ∗ = β HL α ∗ where β HL is the distance between the origin of the standard space and u ∗ and α ∗ is the normalised vectorof direction. Then for each variable U i , one can obtain the importance factor α ∗ i . In the SA framework, let us have X = ( X , ..., X d ), a random vector where the variables are mutuallyindependent and Y = G ( X ), the output of a deterministic code G (). Thus a functional decomposition ofthe variance is feasible, often referred as functional ANOVA: Var[ Y ] = P di =1 V i ( Y ) + P di 1) for i = 1 , .., 4. The failure functionis defined as: G ( X ) = k − X i =1 a i X i where k and a = ( a , a , a , a ) are the parameters of the model. For this numerical example, parametersare set with values k = 16 and a = (1 , − , , P can be given: P = φ − k/ vuut X i =1 a i = φ (cid:18) − √ (cid:19) ≃ . φ ( . ) is the standard normal cumulative distribution function.It is expected that the influence of X i on P uniquely depends on | a i | . The larger the absolute value of thecoefficient, the larger the expected influence. The aim for choosing one non-influential variable X (because a = 0) is to assess if the SA methods can identify this variable as non-influential on the failure probability. In this ideal hyperplane failure surface case, FORM provides an approximated value ˆ P F ORM = 0 . Variable X X X X Importance factor 0 . 018 0 . 679 0 . 302 0 Table 1: Importance factors for hyperplane function The first-order and total indices are displayed in Table 4.2.2. The interpretation of the results is that X and X concentrate most of the variance of the indicator function. At first order, 25% of its variance isexplained by X without any interaction. It should be noted that the total index for X is null, assessingthat this variable does not impact the failure probability. The C.o.v. of the total indices estimators are small,meaning that this method is reproducible and that 6 × points are enough to estimate in an efficient waythe indices S T i . On the other hand, some C.o.v. values for small first order indices are quite high. This effectis well-known in SI estimation problems and can be corrected by improved formulas [15]. Sobol’ Index S S S S S T S T S T S T Mean 0 . . . . . − . . . . . . . . 012 0 . . 013 0 Table 2: Sobol’ indices for hyperplane function The method presented throughout this article is applied on the hyperplane function. As explained in section2, several ways to perturb the input distributions exist. For this case, a mean shifting is first applied, then avariance shifting with fixed mean. A simple calculus gives that the perturbed pdf are Gaussian, respectivelywith the constraint mean and variance 1 for the mean shifting perturbation (see Table 7), and with mean 0and the constraint variance for the variance shifting perturbation.9he MC estimation gives ˆ P = 0 . function calls. For the mean shifting (see (6)), thedomain variation for δ ranges from − δ = 0 cannot be considered as aperturbation. For the variance shifting (see (10)), the variation domain for V per ranges from 1 / 20 to 3 with28 points, where V per = 1 is not a perturbation. The estimated indices are plotted respectively in Figure 2for mean shifting and in Figure 3 for variance shifting. The 95% confidence intervals are plotted around theindices, using the presented asymptotic formulas in Section 3. Mean perturbation indices The indices c S iδ behave in a monotonic way given the importance of theperturbation. The slope at the origin is directly related to the value of a i . For influential variables ( X and X ), the increasing or the decreasing is faster than linear, whereas the curve seems linear for the slightlyinfluential variable ( X ). Modifying the mean with a positive amplitude slightly rises the failure probabilityfor X , highly decreases it for X and increases it for X (Figure 2). The effects are reversed with similaramplitude for negative δ . It can be seen that X has no impact on the failure probability for any perturbation.Those results are consistent with the expression of the failure function. One can see that the confidenceintervals (CI) associated to X and X are fairly well separated, except for the small absolute value of δ .On the other hand, the CI associated to X and X are not separated until absolute value of δ higher than0 . 6. It can be concluded from the observation of the CI is that the impact of variable X and X cannot bedifferentiated, unless a broad change of the mean occurs. Variance perturbation indices Increasing the variance of input X and X increases the failure prob-ability, whereas it decreases when decreasing the variance (Figure 3). Modifying the variance of X and X have no effect on the failure probability. The increasing of the indices is linear for X and X , and thedecreasing of the indices is faster than linear, especially for X . Considering the CI, one can see that theyare well separated for variable X and X , assessing the relative importance of these variables. On the otherhand, the CI associated to X and X are not separated and contain 0. The DMBRSI has brought the following conclusions: when shifting the mean (that is to say the centraltendency in this case), the most influential variable is X , followed by X . X is slightly influential while X is not influential at all. When shifting the variance, variable X is more influential than variable X .Variables X and X have no impact when shifting the variance. We argue that all these information aremuch richer than the ones provided by importance factors and by Sobol’ indices. Indeed, the information areprovided about regions of the input space leading to failure event. This is, in our opinion, more of interestto the practitioner than a ”simple” variable ranking. The Ishigami function is a common test case in SA since it presents high degree of non-linearity with inter-actions between the variables. A modified version of the Ishigami function will be considered in this paper.One has G ( X ) = sin ( X ) + 7 sin ( X ) + 0 . X sin ( X ) + 7where X is a 3 − dimensional vector of independent marginals uniformly distributed on [ − π, π ] . In Figure 4,the failure points (where G ( x ) < 0) are plotted in a 3-d scatterplot.There are 614 failure points on a MC sample of 10 points. Therefore the failure probability here isroughly ˆ P = 6 . × − . The complex repartition of the failure points can be noticed. Those points lay ina zone defined by the negative values of X , the extremal and mean values of X (around − π , 0 and π ), andthe extremal values of X (around − π and π ). The algorithm FORM converges to an incoherent design point (6 . , . , 0) in 50 function calls, giving anapproximate probability of ˆ P F ORM = 0 . 54. The importance factors are displayed in Table 3. The badperformance of FORM is expected given that the failure domain consists in six separate domains and that10igure 2: Estimated indices c S iδ for hyperplane function with a mean shiftingFigure 3: Estimated indices [ S i,V f for hyperplane function with a variance shifting11igure 4: Ishigami failure points from a MC sample12he function is highly oscillant, leading to optimization difficulties. The design point is aberrant, thereforethe FORM results of SA are incorrect. Variable X X X Importance factor 1 e − Table 3: Importance factors for Ishigami function The first-order and total indices, computed with 5 × function calls are displayed in Table 6. The smallvalues of first order indices show that no variable has impact on the variance of the indicator of failure on itsown. The three total indices have relatively high and similar values. This states that all the variables highlyinteract with each other to cause system failure. The SI method is thus non-discriminant in this case. Thelow C.o.v. show that the method is reproducible. Sobol’ Index S S S S T S T S T Mean 0 . . . . . . . . . . . . Table 4: Sobol’ indices for Ishigami function The method presented throughout this article is applied on the thresholded Ishigami function. As for thehyperplane test case, a mean shifting and a variance shifting are applied. The modified distribution whena mean shift is applied on a uniform distribution is given in Table 7. The modified pdf when shifting thevariance and keeping the same expectation is proportional to a truncated Gaussian when decreasing thevariance. When increasing the variance, the perturbed distribution is a symmetrical distribution with 2modes close to the endpoints of the support (see figure 1). As previously, the same MC sample of size 10 (also used to produce Figure 4) is used to estimate the indices with both perturbations. For the mean shifting(see (6)), the variation domain for δ ranges from − V f rangesfrom 1 to 5 with 40 points. Let us recall that the original variance is Var[ X i ] = π / ≃ . . The estimatedindices are plotted respectively in Figure 5 for mean shifting and in Figure 6 for variance shifting. Mean perturbation indices A perturbation of the mean for X and X will increase the failure proba-bility, though the impact for the same mean perturbation is stronger for X ( [ S , − and d S , approximatelyequal respectively 9 . X show that a mean shiftbetween − − d S , approximatively equals − . ). Therefore, Figure 5 leadsto two conclusions. First, the failure probability can be strongly reduced when increasing the mean of thefirst variable X (this is also provided by Figure 4 wherein all failure points have a negative value of X ).Second, any change in the mean for X or X will lead to an increase of the failure probability. The CI arewell separated, except in the − X contains 0 betweenvalues of δ from − . . 5, thus the associated indices might be null in these case. This has to be takeninto account when assessing the relative importance of X . Variance perturbation indices Figure 6 (upper part) shows that a change in the variance has littleeffect on X and X , though the change is of opposite effect on the failure probability. However, consideringthat the indices \ S ,V per ,i and \ S ,V per ,i lie between − . . 4, one can conclude that the variance of thesesvariables are not of great influence on the failure probability. On the other hand, Figure 6 (lower part) showsthat any reduction of Var [ X ] strongly decreases the failure probability, and that an increase of the variance13igure 5: Estimated indices c S iδ for the thresholded Ishigami function with a mean shiftingslightly increases the failure probability. This is relevant with the expression of the failure surface, as X isfourth powered and multiplied by the sinus of X . A variance decreasing as formulated gives a distributionconcentrated around 0. Decreasing Var [ X ] shrinks the concerned term in G ( X ). Therefore it reduces thefailure probability. The CI associated to X are broadly separated from the others. The goal of this test case is to assess the risk of a flood over a dyke for the safety of industrial installations [5].This comes down to model the level of a flood. As a function of hydraulical parameters, many of them beingrandomized to account for uncertainty. From a simplification of the Saint-Venant equation, a flood risk modelis obtained. The quantity of interest is the difference between the level of the dyke and the height of water.If this quantity is negative, the installation is flooded. Hydraulical parameters are the following: Q the flowrate, L the watercourse section length studied, B the watercourse width, K s the watercourse bed frictioncoefficient (also called Strickler coefficient), Z m and Z v respectively the upstream and downstream bottomwatercourse level above sea level and H d the dyke height measured from the bottom of the watercourse bed.The water level model is expressed as: H = QK s B q Z m − Z v L . Therefore the following quantity is considered: G = H d − ( Z v + H ) . Among the model inputs, the choice is made that the following variables are known precisely: L = 5000(m), B = 300 (m), H d = 58 (m), and the following are considered to be random. Q (m . s − ) follows apositively truncated Gumbel distribution of parameters a = 1013 and b = 558 with a minimum value of 0. K s (m / s − ) follows a truncated Gaussian distribution of parameters µ = 30 and σ = 7 . 5, with a minimumvalue of 1. Z v (m) follows a triangular distribution with minimum 49, mode 50 and maximum 51. Z m (m)follows a triangular distribution with minimum 54, mode 55 and maximum 56. A 10 MC sample gives anestimation of the failure probability ˆ P = 8 . × − . 14igure 6: Estimated indices \ S i,V per for the thresholded Ishigami function with a variance shifting. The upperfigure is a zoom where the \ S i,V per axis lies into [ − . , . \ S i,V per . The curves cross for the value of V per that corresponds to the original variance, namely π / .4.1 FORM The algorithm FORM converges to a design point (1 . , − . , . , − . 18) in 52 function calls, giving anapproximate probability of ˆ P F ORM = 5 . × − . The importance factors are displayed in Table 5. Variable Q K s Z v Z m Importance factor 0 . 246 0 . 725 0 . 026 0 . Table 5: Importance factors for flood case FORM assesses that K s is of extremely high influence, followed by Q that is of medium influence. Z v hasa very weak influence and Z m is negligible. It can be noticed that the estimated failure probability is twiceas small as the one estimated with crude MC, but remains in the same order of magnitude. The first-order and total indices are displayed in Table ?? . It can be seen that the estimates of some indicesare negative despite the fact that Sobol’s indices are theoretically positive. The estimation can indeed producenegative results for values close to 0. Sobol’ Index S S S S T S T S T Mean 0 . . . . . . . . . . . . Table 6: Sobol’ indices for Ishigami function Considering the first order indices, Z v and Z m are of null influence on their own. Q is considered tohave a minimal influence (1% of the variance of the indicator function) by itself, and K s explains 24% of thevariance on its own. When considering the total indices, it can be noticed that both Z v and Z m have a weakimpact on the failure probability. On the other hand, Q has a major influence on the failure probability. K s total index is close to one, therefore K s explains (with or without any interaction with other variables)almost all the variance of the failure function.Let us compare the informations provided by the Sobol’ indices with the information provided by theimportance factors. One cannot conclude from the total Sobol’ indices that Z m is not influent whereas theimportance factors assess that this variable is of negligible influence. Additionally, the total Sobol’ indexassociated to K s and Q state that both these variables are of high influence whereas the importance factorsstate that K s is of high influence and Q is of medium influence. The method presented throughout this article is applied on the flood case. Only the mean shifting will beapplied here. The modified pdf are given in Table 7 (Appendix D)and a numerical trick is used to deal withtruncated distributions, as stressed in Appendix E. One can notice that the different inputs follow variousdistributions (unlike the other examples), thus the question of ”equivalent” perturbation arises. It will bediscussed further in Section 5. Here the choice has been made to shift the mean relatively to the standarddeviation, hence including the spread of the various inputs in their respective perturbation. So for any input,the original distribution is perturbed so that its mean is the original’s one plus δ times its standard deviation, δ ranging from − Z v , strongly for Q , and diminishes it slightly for Z m and strongly for K s . This goes the opposite way whendecreasing the mean. In terms of absolute modification, K s and Q are of same magnitude, even if K s has aslightly stronger impact. On the other hand, the effects of mean perturbation on Z m and Z v are negligible.The CI associated to Q and K s are well separated from the others, except in a δ = − . . Z v and Z m overlap. Thus even though the indices seem to have different value, it is notpossible to conclude with certainty about the influence of those variables.16igure 7: Estimated indices c S iδ for the flood case with a mean perturbation The method presented in this paper gives relevant complementary information in addition of traditional SAmethods applied to a reliability problem. Traditional SA methods provide variable ranking, whereas theproposed method provides an indication on the variation in the probability of failure given the variation ofparameter δ . This is useful when the practitioner is interested on which configurations of the problem lead toan increase of the failure probability. This might also be used to assess the conservatism of a problem, if everyvariations of the input lead to decrease in the probability of failure. Additionally, it has two advantages: • The ability for the user to set the most adapted constraints considering his/her problem/objective, • The MC framework allowing to use previously done function calls, thus limiting the CPU cost of theSA, and allowing the user to test several perturbations.We argue that with an adapted perturbation, this method can fulfil the three presented engineers objective. The question of ”equivalent” perturbation arises from cases where all inputs are not identically distributed.Indeed, problems may emerge when some inputs are defined on infinite intervals and when other inputsare defined on finite intervals (such as uniform distributions). For instance, consider a two-dimensionalmodel with one Gaussian distribution and one uniform distribution as inputs. Thus, a mean shift will bea translation for the first input, whereas it will lead to a Dirac distribution in one endpoint for the otherinput. Hence, a mean shift cannot be considered as an ”equivalent” perturbation. A ”relative mean shift”seems promising idea. But if we consider a model with two Gaussian inputs of equal variance 1 and of meanrespectively 1 and 10000. Then, a relative mean shift of 10% will result in Gaussian distributions with mean17espectively 1.1 and 11000, and still variance 1. This counter-example shows that relative mean shift mightnot be an adequate perturbation in terms of ”equivalence”. In the examples given throughout this paper, the perturbations of the inputs left the support of thosevariables unperturbed. However, the practitioner might be interested in the sensitivity on the boundaries ofthe support. The proposed method will be applied with support perturbations in further tests. However, westress that given the estimation method (reverse importance sampling), it is mandataory that the support ofthe perturbed density is included in the support of the original density. Thus one cannot perturb the inputsso that the perturbed support is wider than the original one. Finally, these first results provide some avenues for future research: • Adapting the estimator of the indices S iδ in term of variance reduction and of number of function calls.Further work will be made with importance sampling methods, and possibly subset methods. The useof sequential methods [4] may also be tested. • Second, there is a need to find a way to perturb ”equivalently” several distributions of different natures.To answer this problem, perturbations that are not based upon moment constraints can be proposed.For instance, quantile constraints might be considered. As far as we noticed, in most cases the values ofthe input leading to the failure event comes from the tails of the input distributions. What if these tailswere badly modeled ? Therefore a perturbation based on the quantiles is proposed. More precisely,suppose that the weigth of the left tail is increased. That is to say that the value q becomes for themodified density, for instance the δ quantile. This writes: Z ] −∞ ; q ] ( x ) f mod ( x ) dx = δ. We believe that this kind of perturbation is equivalent with respect to inputs of different natures.A perturbation based on an entropy constraint might also be proposed. The differential entropy ofa distribution can be seen as a quantification of uncertainty [2]. Thus an example of (non-linear)constraint on the entropy can be: − Z f X iδ ( x ) log f X iδ ( x ) dx = − δ Z f X i ( x ) log f X i ( x ) dx. Yet further computations have to be made to obtain a tractable solution of the KL minimizationproblem under the above constraint. Acknoweldgements Part of this work has been backed by French National Research Agency (ANR) through COSINUS program(project COSTA BRAVA noANR-09-COSI-015). We thank Dr. Daniel Busby (IFP EN) for several discus-sions and Emmanuel Remy (EDF R&D) for proofreading. We also thank an anonymous reviewer and theassociated editor for their helpful comments. All the statistical parts of this work have been performed withinthe R environment, including the mistral and sensitivity packages. References [1] A. Antoniadis. Analysis of variance on function spaces. Math. Operationsforsch. Statist. Ser. Statist. ,15(1):59–71, 1984. 182] B. Auder and B. Iooss. Global sensitivity analysis based on entropy. In Safety, Reliability and RiskAnalysis - Proceedings of the ESREL 2008 Conference , pages 2107–2115. CRC Press, september 2008.[3] R.J. Beckman and M.D McKay. Monte-Carlo estimation under different distributions using the samesimulation. Technometrics , 29(2):153–160, 1987.[4] J. Bect, D. Ginsbourger, L. Li, V. Picheny, and E. Vazquez. Sequential design of computer experimentsfor the estimation of a probability of failure. Statistics and Computing , 22(3):1–21, 2011.[5] P. Bernardara, E. de Rocquigny, N. Goutal, A. Arnaud, and G. Passoni. Uncertainty analysis in floodhazard assessment: hydrological and hydraulic calibration. Canadian Journal of Civil Engineering ,37(7):968–979, 2010.[6] T.M. Cover and J.A. Thomas. Elements of information theory 2nd edition (Wiley series in telecommu-nications and signal processing). 2006.[7] I. Csisz´ar. I-divergence geometry of probability distributions and minimization problems. The Annalsof Probability , 3(1):146–158, 1975.[8] T.C. Hesterberg. Estimates and confidence intervals for importance sampling sensitivity analysis. Math-ematical and Computer Modelling , 23(8):79–85, 1996.[9] T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. Reliability Engineering & System Safety , 52(1):1–17, 1996.[10] B. Iooss. Revue sur l’analyse de sensibilit´e globale de mod`eles num´eriques. Journal de la Soci´et´eFran¸caise de Statistique , 152(1):1–23, 2011.[11] M. Lemaire, A. Chateauneuf, and J.C. Mitteau. Structural reliability . Wiley Online Library, 2009.[12] J. Morio. Influence of input PDF parameters of a model on a failure probability estimation. SimulationModelling Practice and Theory , 19(10):2244–2255, 2011.[13] M. Munoz Zuniga, J. Garnier, E. Remy, and E. de Rocquigny. Adaptive directional stratificationfor controlled estimation of the probability of a rare event. Reliability Engineering & System Safety ,92(12):1691–1712, 2011.[14] A. Saltelli. Sensitivity analysis for importance assessment. Risk Analysis , 22(3):579–590, 2002.[15] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, and S. Tarantola. Variance based sensi-tivity analysis of model output. design and estimator for the total sensitivity index. Computer PhysicsCommunications , 181(2):259–270, 2010.[16] A. Saltelli, S. Tarantola, F. Campolongo, and M. Ratto. Sensitivity analysis in practice: A guide toassessing scientific models. 2004. Chichester, England: John Wiley & Sons , 46556:48090–9055.[17] I.M. Sobol. Sensitivity analysis for non-linear mathematical models. Mathematical Modelling and Com-putational Experiment , 1(4):407–414, 1993.[18] A.W. Van der Vaart. Asymptotic statistics . Number 3. Cambridge Univ Pr, 2000. AppendicesA Computation of Lagrange multipliers Let H be the Lagrange function: H ( λ ) = ψ i ( λ ) − K X k =1 λ k δ k . λ ∗ = arg min H ( λ ) . The expression of the gradient of H with respect to the j th variable is ∇ j H ( λ ) = R g j ( x ) f i ( x ) exp( P Kk =1 λ k g k ( x )) dx exp ψ i ( λ ) − δ j . Similarly, the expression of the second derivative of H with respect to the h th and the j th variables is D hj H ( λ ) = R g h ( x ) g j ( x ) f i ( x ) exp( P Kk =1 λ k g k ( x )) dx exp ψ i ( λ ) − R g j ( x ) f i ( x ) exp( P Kk =1 λ k g k ( x )) dx exp ψ i ( λ ) R g h ( x ) f i ( x ) exp( P Kk =1 λ k g k ( x )) dx exp ψ i ( λ ) . This method has been used in this paper for computing the optimal vector λ ∗ when a variance shifting wasapplied. The integrals were evaluated with Simpson’s rule. B Proofs of the NEF properties In this Appendix, the details of the calculus for the Proposition 2.4 are provided. NEF specificities : If the original density f i ( x ) is a NEF, then under a set of K linear constraints on f ( x ), one has : f ( x ) = b ( x ) exp [ xθ − η ( θ )] , thus : f δ ( x ) = f ( x ) exp " K X k =1 λ k g k ( x ) − ψ ( λ ) The regularization constant from (4) can be written as: ψ ( λ ) = log Z b ( x ) exp " xθ + K X k =1 λ k g k ( x ) − η ( θ ) dx (14)If the integral on (14) is finite, f δ exists and is a density. Mean shifting With a single constraint formulated as in (6), (14) becames : ψ ( λ ) = log Z b ( x ) exp [ xθ + λx − η ( θ )] dx = log Z b ( x ) exp [ x ( θ + λ ) − η ( θ ) + η ( θ + λ ) − η ( θ + λ )] dx if η ( θ + λ ) is well defined. ψ ( λ ) = ( η ( θ + λ ) − η ( θ )) + log (cid:20)Z b ( x ) exp [ x ( θ + λ ) − η ( θ + λ )] (cid:21) dx = η ( θ + λ ) − φ ( θ )since b ( x ) exp [ x ( θ + λ ) − η ( θ + λ )] = f θ + λ ( x )with notation from (2.4), is a density of integral 1. Thus f δ ( x ) = b ( x ) exp [ xθ − φ ( θ )] exp [ λx − η ( θ + λ ) + η ( θ )]= b ( x ) exp [ x [ θ + λ ] − η ( θ + λ )] = f θ + λ ( x )Thus the mean shifting of a NEF of CDF η ( . ) results in another NEF with mean η ′ ( θ + λ ) = δ (constraint)and variance η ′′ ( θ + λ ). 20 ariance shifting With a single constraint formulated as in (10), using (14), the new distribution has fordensity: f δ ( x ) = b ( x ) exp (cid:2) xθ + xλ + x λ − ψ ( λ ) − η ( θ ) (cid:3) Since λ is known or computed, and θ is also known, consider the variable change z = √ λ x assuming λ isstrictly positive (the variable change is z = √− λ x if λ is strictly negative). Thus, f δ ( x ) = b ( z √ λ ) exp (cid:2) z (cid:3) exp (cid:20) z √ λ ( θ + λ ) − ψ ( λ ) − η ( θ ) (cid:21) = exp (cid:20) η (cid:18) ( θ + λ ) √ λ (cid:19) − η ( θ ) − ψ ( λ ) (cid:21) c ( z ) exp (cid:20) z ( θ + λ ) √ λ − η (cid:18) ( θ + λ ) √ λ (cid:19)(cid:21) with c ( z ) = b ( z √ λ ) exp (cid:2) z (cid:3) . By (4), ψ ( λ ) = log Z b ( x ) exp (cid:2) xθ + xλ + x λ − η ( θ ) (cid:3) dx = log Z b ( z √ λ ) exp (cid:2) z (cid:3) exp (cid:20) ( θ + λ ) √ λ z − η ( θ ) + η (cid:18) ( θ + λ ) √ λ (cid:19) − η (cid:18) ( θ + λ ) √ λ (cid:19)(cid:21) dx = (cid:18) η (cid:18) ( θ + λ ) √ λ (cid:19) − η ( θ ) (cid:19) + log Z c ( z ) exp (cid:20) ( θ + λ ) √ λ z − η (cid:18) ( θ + λ ) √ λ (cid:19)(cid:21) dx = η (cid:18) ( θ + λ ) √ λ (cid:19) − η ( θ )Thus one has : f δ ( x ) = c ( z ) exp (cid:20) z ( θ + λ ) √ λ − η (cid:18) ( θ + λ ) √ λ (cid:19)(cid:21) thus the variance shifting of a NEF results in another NEF parameterized by ( θ + λ ) √ λ . C Proofs of asymptotic properties C.1 Proof of Lemma 3.1 Under assumption (i) , we have Z Supp( f iδ ) { G ( x ) < } f iδ ( x i ) f i ( x i ) f ( x ) d x ≤ Z Supp( f iδ ) f iδ ( x i ) dx i = 1 . So that, the strong LLN may be applied to ˆ P iδN . Defining σ iδ = Var (cid:20) { G ( X ) < } f iδ ( X i ) f i ( X i ) (cid:21) , (15)one has σ iδ = Z Supp( f i ) { G ( x ) < } f iδ ( x i ) f i ( x i ) d Y j = i f j ( x j ) d x − P iδ < ∞ under Condition (ii) .Therefore the CLT applies: √ N σ − iδ (cid:16) ˆ P iδN − P iδ (cid:17) L −→ N (0 , . Under assumption (ii) , the strong LLN applies to ˆ σ iδN . So that, the final result is straightforward usingSlutsky’s lemma. 21 .2 Proof of Proposition 3.1 First, note that E h b P c P iδ i − P P iδ = E " N N X n =1 { G ( x n ) < } ! N X n =1 { G ( x n ) < } f iδ ( x ni ) f i ( x ni ) ! − P P iδ = 1 N E N X n =1 (cid:2) { G ( x n ) < } (cid:3) f iδ ( x ni ) f i ( x ni ) + N X n =1 N X j = i { G ( x n ) < } { G ( x j ) < } f iδ ( x ji ) f i ( x ji ) − P P iδ = 1 N [ N P iδ + N ( N − P P iδ ] − P P iδ = 1 N ( P iδ − P P iδ ) . Assuming the conditions under which Lemma 1 is true, the bivariate CLT follows withΣ iδ = (cid:18) P (1 − P ) P iδ (1 − P ) P iδ (1 − P ) σ iδ (cid:19) . Each term of this matrix can be consistently estimated, using the results in Lemma 1 and Slutsky’s lemma. D Summary Table with modified distributions for mean shiftE Numerical trick to work with truncated distribution In the case where a mean shifting is considered on a left truncated distribution. We present a tip that canhelp to compute λ ∗ .The studied trucated variable Y T has distribution f Y T . Let us denote Y ∼ f Y the corresponding non-truncated distribution. The truncation occurs for some real value a . This truncation may happen for somephysical modelling reason. One has: f Y T ( y ) = 11 − F ( a ) 1 [ a, + ∞ [ ( y ) f Y ( y ) . The formal definition of M Y T ( λ ) the mgf of Y T for some λ is: M Y T ( λ ) = 11 − F Y ( a ) Z + ∞ a f Y ( y ) exp [ λ y ] dy. Let us recall that we are looking for λ ∗ such as: δ = M ′ Y T ( λ ∗ ) M Y T ( λ ∗ ) = R + ∞ a yf Y ( y ) exp [ λ y ] dy R + ∞ a f Y ( y ) exp [ λ y ] dy . (16)When the expression does not take a practical form, one can use numerical integration to estimate theintegral terms. Unfortunately, for some heavy tailed distribution (for instance Gumbel distribution), thisnumerical integration might be complex or not possible. This is due to the multiplication by an exponentialof y . The following tip helps to avoid such problems. Denoting M Y ( λ ) the mgf of the non-truncateddistribution, one can remark that: M Y ( λ ) = Z + ∞−∞ f Y ( y ) exp [ λ y ] dy = Z a −∞ f Y ( y ) exp [ λ y ] dy + Z + ∞ a f Y ( y ) exp [ λ y ] dy Thus another expression for M Y T ( λ ) is: M Y T ( λ ) = 11 − F Y ( a ) (cid:20) M Y ( λ ) − Z a −∞ f Y ( y ) exp [ λ y ] dy (cid:21) . . )is the cdf of the standard normal distribution, and φ ( . ) is its pdf. Original distribution Modified distribution Modified pdf f iδ Link between λ ∗ and δ NEF( θ ) NEF ( θ + λ ∗ ) f iδ ( x i ) = b ( x i ) exp [ x i [ θ + λ ∗ ] − η ( θ + λ ∗ )] η ′ ( θ + λ ∗ ) = δ Special case of NEF: N ( µ, σ ) N ( δ, σ ) f iδ ( x i ) = σ √ π exp (cid:20) − (cid:16) x i − δσ (cid:17) (cid:21) λ ∗ = δ − µσ Uniform distribution: U [ a,b ] ∝ truncated exponential f iδ ( x i ) = λ ∗ e λ ∗ b − e λ ∗ a [ a,b ] ( x i ) e λx i δ = b − a ) e λ ∗ b ( λ ∗ b − ) + e λ ∗ a ( − λ ∗ a ) λ ∗ ( e λ ∗ b − e λ ∗ a )Left Tr Gaussian N T ( µ, σ, a ) N T ( µ + σ λ ∗ , σ, a ) f iδ ( x i ) = [ a, + ∞ [ ( x i )1 − F ( a ) 1 σ √ π exp (cid:20) − (cid:16) x i − µ − σ λ ∗ σ (cid:17) (cid:21) δ = µ + σ λ ∗ − σ φ (cid:18) a − ( µ + σ λ ∗ ) σ (cid:19) − Φ (cid:18) a − ( µ + σ λ ∗ ) σ (cid:19) Triangle T ( a, b, c ) – f iδ ( x i ) = exp( x i λ ∗ − ψ ( λ ∗ )) f ( x i ) δ = ( a − λ ∗ ) e λ ∗ a ( b − c )+( b − λ ∗ ) e λ ∗ b ( c − a )+( c − λ ∗ ) e λ ∗ c ( a − b ) e λ ∗ a ( b − c )+ e λ ∗ b ( c − a )+ e λ ∗ c ( a − b ) Left Tr Gumbel G T ( µ, β, a ) – f iδ ( x i ) = exp( x i λ ∗ − ψ ( λ ∗ )) f ( x i ) δ = M ′ Y ( λ ∗ ) − R a −∞ yf Y ( y ) exp [ λ ∗ y ] dyM Y ( λ ∗ ) − R a −∞ f Y ( y ) exp[ λ ∗ y ] dy with : M Y ( λ ∗ ) = Γ (1 − β ) exp [ λ ∗ µ ] M ′ Y ( λ ∗ ) = Γ (1 − β ) exp [ λ ∗ µ ] (cid:2) µ − β ̥ (0) (1 − λ ∗ ) (cid:3) he integral term is much smaller in the left heavy tailed distribution case. Therefore the numerical integra-tion (for instance using Simpson’s method) is much more precise or became possible.The same goes for M ′ Y T ( λ ) which has alternative expression: M ′ Y T ( λ ) = 11 − F Y ( a ) (cid:20) M ′ Y ( λ ) − Z a −∞ yf Y ( y ) exp [ λ y ] dy (cid:21) . Finally, another form of 16 is: δ = M ′ Y ( λ ) − R a −∞ yf Y ( y ) exp [ λ y ] dyM Y ( λ ) − R a −∞ f Y ( y ) exp [ λ y ] dy . (17)This alternative expression may lead to more precise estimations of λ ∗ when M Y ( λ ) and M ′ Y ( λλ