Model’s output variance can increase when input variance decreases : a sensitivity analysis paradox?

Pascal Pernot1,4,a, Michèle Désenfant2,4, et François Hennebelle3,4 1Laboratoire de Chimie Physique, UMR8000 CNRS/Univ. Paris-Sud, Orsay, France 2Laboratoire National de Métrologie et d’Essais (LNE), 1 rue Gaston Boissier, 75724 Paris Cedex 15, France 3Laboratoire Electronique, Informatique et Image, UMR6306 CNRS/Université de Bourgogne, Avenue des Plaines de l’Yonne, 89000 Auxerre, France 4Réseau "Mesures, Modèles et Incertitudes", Mission pour l’Intersciplinarité CNRS


Introduction
The current paradigm in uncertainty propagation would be "the smaller the uncertainty on the inputs of a model, the smaller the uncertainty on this model's output(s)".This is certainly true in the context of the Law of Propagation of Uncertainty (LPU) [1], valid for models approximately linear.
For nonlinear models, uncertainty estimation by propagation of distributions [2] can produce less intuitive results.In fact, for some models, the inverse behavior is observed : reducing the uncertainty u X of an input parameter X can increase the uncertainty u Y on the model's output Y = f (X). 1  This was observed for the Ishigami function -a classical benchmark model for sensitivity analysis methods-by the authors of the variance gradients (VG) sensitivity analysis method [3].The ability of the VG method to detect such cases is remarkable, because the most popular sensitivity analysis methods, based on variance decomposition, are blind to the phenomenon.
The Ishigami function being rather involved, we found interesting to search elementary model functions for which negative VG indices can be obtained and appreciate if these models can be of interest in measurement uncertainty a. e-mail: pascal.pernot@u-psud.fr1.We guarantee the reader the immediate attention of incredulous coworkers when (s)he states this proposition, especially from those trained to LPU.
studies.We also wanted to assess the role of the inputs probability density function (PDF).
Faced to this counter-intuitive phenomenon, one can also question the pertinence to use the variance as a descriptor of uncertainty, and we studied the behavior of a more robust statistics, the (symmetric) 95 percent interquantile range IQ 95 .
In the next section we briefly review the VG concept and present the computational details.In Section 3, we provide examples of elementary models presenting negative VG indices, or not, and study them by propagation of distributions [2].These case studies nourish the general discussion in Section 4.

Variance Gradients
Variance gradients have been introduced [3] as sensitivity measures for nonlinear functions of n independent random variables Y = f (X 1 , . . ., X n ).VG indices G X i are defined as They are dimensionless and vary in ] − ∞, +∞[.For non-linear models, the VG indices do not sum to one, i G X i 1.Besides the possible negativity of VG indices, this is a major difference with sensitivity analysis measures based on variance decomposition.However, for linear models, Eq. 1 recovers the standard decomposition of variance, i.e.
is the portion of the output variance due to X i [1].In this case, G X i indices vary in [0, 1] and sum to 1.For all model types, G X i quantifies the effect of a change in the variance of X i on u 2 Y : i.e., a change of 100 × p percent of u 2 X i would cause a change of 100 × p × G X i percent of u 2 Y .For non-linear models, this relation is expected to hold for small variations (|p| 1), whereas for linear models, it is valid for any p ≥ −1.This points to a convenient and little known interpretation of the usual analysis of variance coefficients for linear models.
For non-linear models, the VG indices provide information that cannot be obtained by sensitivity measures based on variance decomposition.For instance, as reported by Campanelli et al. [3], and verified below (Sect.3.4), VG indices provide a correct ranking when considering small inputs variance reduction.By contrast, global sensitivity analysis tools such as the Sobol method [4], inform us on the output variance reduction to be expected from a full cancellation of the variance of each input.

Numerical applications
The results presented next were obtained by Monte Carlo uncertainty propagation [2].The VG indices were computed using Monte Carlo samples generated for the propagation of distributions.All calculations have been made in R [5], with functions gumS1 and vgSA from our teaching package rgumlib (https ://github.com/ppernot/rgumlib).Sobol indices have been calculated with the soboljansen function of package sensitivity [6].A R code example for the VG analysis of onedimensional models is provided in the Appendix.

Numerical experiments
Several models have been selected to display the main features linked to negative VG indices.They are mostly one-dimensional elementary functions (symmetric-exponential, hyperbolic tangent and sine).We also included the Ishigami function as an example of a multi-dimensional model.

The Symmetric-exponential model
Let us first consider the following model where X is uniformly distributed, with half-width a = √ 3u X , and centered in 0 (Fig. 1).
For this model, we derived analytical results for the quantities of interest, in order to validate the Monte Carlo  estimations :

Numerical results
Uncertainty propagation has been performed for a series of u X values by Monte Carlo, with 10 5 runs per point.The variations of G X and u Y with u X are presented in Fig. 2-(top).They are in perfect agreement with the analytical results.The symmetric 95 percent coverage interval IQ 95 estimated from the Monte Carlo samples has also been reported in Fig. 2-(top).
The u Y curve presents two regimes : for small values of u X , u Y increases with u X (which would be the "intuitive" LUP-compliant behavior), then it passes by a maximum and decreases.Simultaneously, the VG index G X decreases continuously from 1 to negative values.The maximum of u Y corresponds to the value of u X where G X changes its sign (u X 1.9).The IQ 95 statistics follows a similar behavior as u Y .As the PDF of Y is non-symmetric, we also evaluated the behavior of Shortest Probability Intervals [7], with the same conclusion (not shown).
Similar results have been obtained for the Gaussian model Y = exp(−X 2 ) (not shown), and can be extended to any single peak/dip shape function with a flat baseline.Y by 1 percent, confirming the interpretation of VG indices (Eq.3).

Output Probability Density Function
We built histograms of Y samples to approximate PDF(Y) and collected them to design a 2D density plot, as a function of u X and Y (Fig. 3).The PDF of Y is unimodal for all values of u X .As u X increases, the maximum probability density moves from the high values of Y to the lower ones.For large values of u X , there is an accumulation of the propability density at the lower values of Y which induces a decrease of variance.We will consider the value Y = 0 as an accumulation point for PDF(Y).

Effect of X distribution shape
To evaluate the impact of the input's PDF, Monte Carlo estimates of G X , u Y and IQ 95 have been generated with a normal law of mean 0 and standard deviation u X , The results are reported in Fig. 2-(bottom).All quantities present similar behavior as observed for the uniform distribution.The maximum of u Y is reached for a value of u X very close to the one observed for the uniform distribution.For the present model, the shape of PDF(X) has a minor influence on the results of VG analysis.

The Hyperbolic tangent function
The previous example hints at the role of accumulation points in variance reduction.We consider now a model with two distinct asymptotic accumulation points in Y space, the hyperbolic tangent function (Fig. 4) Monte Carlo estimation of the quantities of interest are reported in Fig. 5 : u Y increases monotonically and reaches an asymptote, while G X decreases monotonically to zero.We plotted also a line depicting the variance increase in the LPU/GUM approximation, verifying that the linear approximation is valid for small values of u X .
As u X increases, the points of the output sample symmetrically accumulate towards the limits of Y.The asymptotic shape with u X of the PDF of Y should then be well represented by a Beta distribution lim Beta(0.5(Y − 1); α, β).
Knowing that the asymptotic variance for the standard Beta distribution when α and β vanish is 0.25, one can estimate the asymptotic value for u 2 Y to be close to 1, in agreement with the Monte Carlo simulations.IQ 95 reaches more rapidly its asymptotic value of 2, also in agreement with the properties of the Beta distribution.

The Sine model
Let us now consider the Sine model (Fig. 6) This model differs from the previous ones in the sense that there is no asymptotic accumulation point(s).Analytical expressions of the mean, variance and VG index are Monte Carlo estimation of the quantities of interest are reported in Figs.7-(top) and 8.The Monte Carlo values match perfectly the analytical results (Eqns.16 and 17).
The value of u Y increases with u X until a = √ 3u X becomes slightly larger than 3π/4 and then displays damped oscillations with period a = π around its asymptotic value of 1/ √ 2. G X oscillates with the same period and a phase shift of −π/4.For large values of u X , G X is dominated by the cos 2a term, without damping of its oscillations.IQ 95 presents oscillations of smaller amplitude.
Looking at the 2D density plot of PDF(Y) (Fig. 8), the first maximum of u Y corresponds to an accumulation of the probability density around the extremal values of Y, such as for the tanh model.However, for larger values of u X , more points are sampled from the central values of Y, leading to a decrease of u Y .The oscillatory pattern is damped, as each new period adds a decreasing relative contribution to PDF(Y).

Effect of X distribution shape
Monte Carlo estimates of G X , u Y and IQ 95 have been generated with a normal distribution for X (Eq.9).The results are reported in Fig. 7-(bottom).The oscillatory pattern is almost fully damped.G X does not become significantly negative (i.e.beyond the sampling uncertainty).The spread of the normal distribution overlaps simultaneously regions of decrease and increase observed for the uniform distribution, with a practically null budget.
This model shows that local accumulation points can be a source of negative VG indices, but the latter are more sensitive to the shape of the distribution of X than for asymptotic accumulation points.

The Ishigami function
In their seminal paper, Campanelli et al. [3] applied the VG method to the Ishigami function and obtained a negative index.We reproduce here this experiment and analyze the results in the light of the preceding models.
The VG analysis reveals a negative coefficient for X 1 (Table 1).Its occurrence can be directly interpreted by the presence of local accumulation points, as seen in the Sine model.Considering the dependence of Y on X 2 , negative VG indices are also likely to occur periodically for this variable, for instance at The ranking of the variables obtained by the first order Sobol indices and by the VG indices are remarkably different (Table 1).Sobol indices suggest that reducing the variance of X 1 or X 2 should have the largest impact, whereas VG indices suggest that there is more to be expected from reducing the variance of X 3 .As explained above, this illustrates the different purposes of both indices : Sobol indices refer to a full cancellation of each input's variance, whereas VG indices refer to a small reduction of their variance.
To assess the validity of the VG ranking, we performed a 10 percent variance reduction of all variables individually 2 .The results are reported in Table 1.In each case, the relative variation of variance Δu 2 Y /u 2 Y is computed (1) from the Monte Carlo samples ; and (2) by the VG estimator, i.e.Δu 2 Y /u 2 Y 0.1×G X i .The agreement is very good in all cases, and the most efficient reduction (18 %) is indeed observed for X 3 , which has the smallest first order Sobol index.

Discussion
Examples above show that negative VG indices are not a property of exotic models, and that the prediction of variance variation derived from the VG indices are confirmed by simulations.
We discuss next two points of interest : (1) the link of negative VG indices with the presence of accumulation points in PDF(Y) ; and (2) the need to complement the GUM-Supp1 with basics of sensitivity analysis.

Importance of accumulation points
We have shown the role of accumulation points in the PDF(Y) for the occurrence of negative VG indices, accumulation points linked to asymptotic or local null derivatives of the model.We have not performed an exhaustive exploration of models, but salient features of the few models studied above enable to characterize nonlinear models suitable to display negative VG indices.
The main point is that the model should present null partial derivatives in the variation range of X.These null derivatives correspond to values of Y for which peaks (accumulation points) might appear in the PDF of Y.We have seen that a single asymptotic accumulation point will lead to negative VG indices.By contrast, if the distribution of Y is drawn towards two asymptotic accumulation points, u Y is expected to increase monotonously with u X .The presence of local accumulation points might also induce negative VG indices.Therefore, the sampling of any nonmonotonous model is liable to produce negative VG indices.
Note that these features will be most prominent for input variables with sharply bounded distributions, e.g., uniform PDFs.In particular, the presence of local accumulation points might be totally damped by a normal input PDF.

Scope
It is to be expected that users studying non-linear models in the presence of large inputs uncertainty by Monte Carlo uncertainty analysis have fair chances to encounter them.The combination of non-linear models and large input uncertainties has been reported for instance in the chemical modeling of planetary atmospheres [8], or in astrochemistry [9].
In metrology, where uncertainties are more controlled, negative VG indices should be infrequent.

Interest of VG for the GUM Supp1
We have shown above that the decrease of output variance with the increase of input variance is not an artefact observed for exotic models.We might still question if this is an artifact of the variance itself as an estimator of uncertainty, and, perhaps, revise the definition of uncertainty.We have observed on all the example models that more robust statistics, such as IQ 95 display the same behavior as the variance.In fact, we have seen that the shape of the PDF(Y) itself is affected by the presence of null derivatives in the model.We have therefore to consider this phenomenon as normal, even if reconciliation with the LPU paradigm is impossible.This can be considered as an added value of the propagation of distributions.
It seems important to us that users of the GUM-Supplement 1 [2] be warned of this possibility.Variance Gradients should be proposed as a standard tool of sensitivity analysis (an elementary implementation in R is provided in the Appendix).Beyond the merit of revealing negative VG indices, it provides in all cases a predictive assessment of the importance of variables in the uncertainty budget.The impact of small reductions of variance of the uncertain parameters might be of higher interest to the measurement community than the rankings provided by other sensitivity analysis methods.The top sub-table provides the mean value, standard deviation (SD), variance gradients (VG) and Sobol first order (S) and total order (T) indices for the original model.The three following sub-tables report results for a reduction of 10 percent on the variance of each parameter X 1 , X 2 and X 3 .The relative variation of the Monte Carlo output variance, Δu 2 Y /u 2 Y , is reported for each case and compared to the prediction of the VG analysis.

Mean
Modified values of interest are in bold characters.

Figure 2 .
Figure 2. Symmetric-exponential model : output uncertainty u Y , symmetric 95 percent coverage interval IQ 95 and VG index G X as a function of the input distribution standard uncertainty u X , where X is uniformly distributed (top) or normally distributed (bottom) and centered in 0. Lines : analytical ; bullets : Monte Carlo.

Figure 3 .
Figure 3. Symmetric-exponential model : 2D density plot of the probability density function of Y as a function of u X .Cuts at selected u X values are presented in the right panel.The corresponding values of the variance u Y and VG index G X are reported above the cuts.

Figure 4 .
Figure 4. Curve for the Tanh model and probability density function PDF(X) for u X = 1/ √ 3.

Figure 5 .
Figure 5. Tanh model : u Y , IQ 95 and G X as a function of the input uncertainty u X , where X is uniformly distributed and centered in X = 0.0.

Figure 6 .
Figure 6.Model curve Y = sin(X) and probability density function g(X) for a = 2.

Figure 7 .
Figure 7. Sine model : u Y , IQ 95 and G X as a function of u X ; (top) X uniformly distributed and centered in 0 ; (bottom) X normally distributed and centered in 0. The full line represents the value of u Y estimated by combination of variances (GUM).The vertical dashed lines represent values of a = √ 3u X multiples of π/4.

Figure 8 .
Figure 8. Sine model : 2D density plot of the probability density function of Y as a function of u X .Cuts at selected u X values are presented in the right panel.The corresponding values of the variance u Y and VG index G X are reported above the cuts.
Monte Carlo distribution propagation, Sensitivity Analysis and Variance Gradients tests for the Ishigami function.