Uncertainty and sensitivity analysis of functional risk
curves based on Gaussian processes
Bertrand Iooss and Loı̈c Le Gratiet

arXiv:1704.00624v2 [math.ST] 25 Jul 2017

EDF Lab Chatou, 6 Quai Watier, 78401 Chatou, France

Abstract
A functional risk curve gives the probability of an undesirable event as a function of the value of a critical parameter of a considered physical system. In
several applicative situations, this curve is built using phenomenological numerical models which simulate complex physical phenomena. To avoid cpu-time
expensive numerical models, we propose to use Gaussian process regression to
build functional risk curves. An algorithm is given to provide confidence bounds
due to this approximation. Two methods of global sensitivity analysis of the
models’ random input parameters on the functional risk curve are also studied. In particular, the PLI sensitivity indices allow to understand the effect of
misjudgment on the input parameters’ probability density functions.
Keywords: Computer experiments, Metamodel, Gaussian process, Sobol’
indices, Structural Reliability, Non Destructive Testing, Probability of
Detection

1. Introduction
In industrial practice, the estimation of a functional risk curve (FRC) is
often required as a quantitative measure of the safety of a system [1]. A FRC
(also called “functional risk criterion” or “fragility curve” in some application
domains [2]) gives the probability of an undesirable event as a function of the
value of a critical parameter of a considered physical system. The estimation of
this curve sometimes relies on deterministic phenomenological computer models
which simulate complex physical phenomena. Uncertain input parameters of
Preprint submitted to Reliability Engineering & System Safety

July 26, 2017

this computer code are modeled as random variables while one of its scalar
output variables becomes the studied random variable of interest.
As it is based on a probabilistic modeling of uncertain physical variables
and their propagation through a numerical model, this problem can be directly
related to the uncertainty management methodology in numerical simulation
[1, 3]. This methodology proposes a generic framework of modeling, calibrating,
propagating and prioritizing uncertainty sources through a numerical model (or
computer code). Indeed, investigation of complex computer code experiments
has remained an important challenge in all domains of science and technology,
in order to carry on simulations as well as predictions, uncertainty analysis
or sensitivity studies [4, 5]. In this paper, the physical numerical model G is
expressed as
Y = G(a, X) = G(a, X1 , . . . , Xd ) ,

(1)

with a ∈ R the critical input parameter of the FRC and X ∈ Rd the random
input vector of dimension d and Y ∈ R a scalar model output. In this paper,
parameters a and X are considered as independent. However the presented
results hold when a and X are dependent and the distribution of X conditionally
to a is known.
However, standard uncertainty treatment techniques require many model
evaluations and a major algorithmic difficulty arises when the computer code
under study is too time expensive to be directly used. For cpu-time expensive
models, one solution consists in replacing the numerical model by a mathematical approximation, called a response surface or a metamodel. Several statistical
tools based on numerical design of experiments, efficient algorithms for uncertainty propagation and metamodeling concepts will then be useful [4]. In this
paper, Gaussian process regression [6] is used as a metamodeling technique (surrogate model of the computer code) and applied in the particular context of a
FRC as a quantity of interest. Numerous studies have shown that this model
provides a powerful statistical framework to compute an efficient predictor of a
deterministic computer code response [7, 8].

2

Associated to the estimation of the quantity of interest after the uncertainty
propagation step, the sensitivity analysis step is performed to determine those
parameters that mostly influence the model response [9]. In particular, global
sensitivity analysis methods take into account the overall uncertainty ranges of
the input parameter (see [10] for a recent review). Several works have focused on
estimating global sensitivity indices (especially the variance-based ones, known
as the Sobol’ indices) using Gaussian process models [11, 12, 8]. In this paper, we
present global sensitivity indices attached to the whole FRC (firstly introduced
in [13] in a particular context), while showing how to develop them with a
Gaussian process model. We also focus on a recently developed method, the
Perturbed-Law based sensitivity Indices (PLI) of [14], which seems promising
in terms of computational efficiency and interpretation.
Two examples of FRC largely used in some industrial safety practices are
given in the next section. The third section describes the Gaussian process way
to model and estimate a FRC. The fourth section develops the sensitivity indices
adapted to FRC.

2. Motivating examples and classical methods
2.1. Probability of detection curves
In several industries, the Probability Of Detection (POD) curve is a standard
tool to evaluate the performance of Non Destructive Testing (NDT) procedures.
The goal is to assess the inspection capability for the detection of flaws in the
inspected structure. For instance, new aeronautic regulations require to perform
appropriate damage tolerance assessments to address the potential for failure
due to critical parts of material and manufacturing variability. This imposes
enhanced requirement on POD sizing. Another example is the eddy current
non destructive examination process which is used to ensure integrity of steam
generator tubes in nuclear power plants [15, 16, 13]. Experimental campaigns
and simulated experiments are used in order to provide POD curves, which

3

can be used to demonstrate the performance of inspection process to regulatory
authorities.
In practice, the high costs of the implementation of experimental POD campaigns combined with the continuous increase in the complexity of configuration
make them sometimes unaffordable. To overcome this problem, it is possible to
resort to numerical simulation of NDT processes. This approach is commonly
called MAPOD for “Model Assisted Probability of Detection” [17, 18, 13].
More precisely, in the POD context, the problem is formulated as follows.
Given a threshold s > 0, a flaw is considered to be detected when Y > s with Y
the signal amplitude and s a detection threshold. Therefore the one dimensional
POD curve is denoted by:
∀a > 0

POD(a) = P (Y > s | a) ,

(2)

where a is the POD parameter of interest (for example the size flaw) and X are
the other input parameters. When simulated experiments are used to build the
POD (see for example [13]), Y is a scalar output of a numerical model (see Eq.
1)), where a is determined by its bounds and X is a random vector defined by
its joint probability density function (pdf).
In Sections 3 and 4 of this paper, we will illustrate the mathematical methods
for FRC building and analysis on a case study on the eddy current examination
of the wear induced on steam generator tubes by anti-vibration bars [15]. For
details, the reader is referred to [16, 13]. The study involves the finite element
computer code C3D [19]. The NDT process under study is the examination
of anti-vibration bars (AVB) subject to wear of steam generator tubes inside a
nuclear power plant [15]. The computer code used simulates by finite elements
the eddy-current propagation phenomena. The critical input parameter a is the
maximum of the two flaw sizes while five other input parameters (vector X)
of the computer code are also considered (with their associated pdf): E (pipe
thickness), ebav1 and ebav2 (lengths of the gap between the AVB and the first
and second flaws), h11 and h12 (first and second flaw heights). The output Y
of the computer code is the signal amplitude, and the threshold of detection s
4

is fixed at a specific value.
2.2. Seismic fragility curves
The “fragility curve” is a popular functional risk criterion, commonly used in
many engineering fields [2]. It describes the probability that the actual damage
to a structure exceeds a damage threshold, when the structure is assigned to a
specified load intensity. For instance, the seismic fragility curve, which concerns
systems subject to earthquakes, is of particular interest in nuclear safety studies
(see [20]).
In the case of seismic risk assessment the load is usually expressed as a scalar
characteristic of a seismic signal, typically the horizontal Peak Ground Acceleration (PGA), which is a common choice in civil engineering. Several parameters
and phenomena, distinct from the PGA, also influence the load: for each PGA
value, the occurrence of the damage event is random. The fragility curve may
then be interpreted as the cumulative distribution function of the “structural
capability”, i.e. the maximum load the structure under investigation can bear
without damage [21, 22]. Fragility curves are useful tools in structural analysis
as they provide more complete information than the usual failure probability
(established for a reference value of the load only).
In standard practice, the assessment is made either following an approach
entirely based on available expertise or by the statistical analysis of actually
observed or simulated data. As actual damage data may be scarce due to the
rarity of severe earthquakes liable to generate damages on highly safe structures,
observations are generated by mock-up or (most often) numerical experiments.
In this context, a fragility function is expressed as
∀a > 0

Frag(a) = P (Y > s | a) ,

(3)

where Y is the variable which characterizes the structural damage, s = Ymax is
the maximal capacity of the structure and a is an intensity measure characterizing the ground motion severity (e.g. peak ground acceleration). When the
fragility curve is estimated by numerical simulation (see for example [20]), Y is
5

the numerical model output (see Eq. (1)), a is determined by its bounds and
the random vector X is defined by its joint pdf.
2.3. Classical methods
The standard practice for establishing a fragility curve consists in assuming
a lognormal distribution for Y [23]. The classical approach for determining a
POD curve is called the Berens model and relies on assuming a simple linear
relation (with Gaussian residuals) between log(Y ) and a [24]. In fact, both cases
are equivalent and lead to the following form for the FRC:


a−α
P (Y > s|a) = Φ
,
β

(4)

where Φ is the standard Gaussian distribution function, and α and β are two
constants to be estimated (by the maximum likelihood estimation for example)
as functions of the Y -data sample and s. The parameter a can also be expressed
in terms of a logarithm in this expression. We note that if the Berens method
is applied in a FRC framework (e.g. with Y the maximal displacement of a
building and a the pick ground acceleration of a seismic signal) the resulting
PoD curve is the classical FRC.
As shown in [13], this model can easily be improved by considering a BoxCox transformation [25] instead of a simple logarithmic one on Y : Ỹ =

Y λ −1
λ .

λ is determined by maximum likelihood as the real number that offers the best
linear regression of Ỹ regarding the parameter a. The model is now based on
Ỹ and is defined as
Ỹ (a) = β0 + β1 a + ,

(5)


with  the model error such as  ∼ N 0, σ2 . The maximum likelihood method
provides the estimators βˆ0 , βˆ1 and σˆ . Hence the model implies the following


result: ∀a > 0, Ỹ (a) ∼ N βˆ0 + βˆ1 a, σˆ 2 . Then the value of the FRC can be
estimated from βˆ0 , βˆ1 and σˆ , as its confidence interval from the property of the
maximum likelihood estimators. Note that it is possible to relax the Gaussian
hypothesis on  by using non-parametric approaches [26, 13]

6

3. Gaussian-process based functional risk curve estimation
Classical methods described in Section 2.3 use a simple model between Y
and a without modeling the functional dependence of Y on the other uncertain
variables X. As in our case data yN (sample of Y of size N ) are obtained from
numerical models, the input variables X are controlled and can be introduced in
the FRC determination process. However, with cpu-time expensive numerical
models, model evaluations can be somewhat limited and only small samples of
Y can be obtained. As shown in [27, 21, 13], the use of a metamodel, which is a
mathematical approximation of the computer code [4], is useful. We highlight
that it is important in a safety framework to evaluate accurately the uncertainty
on the surrogate model. We propose in this section a general formulation for
FRC that uses Gaussian process regression as a metamodel [28, 6]. We propose
in this section a general formulation for FRC that uses the Gaussian process
regression metamodel [28, 6]. Through the Gaussian assumption, this metamodel provides a convenient framework to perform uncertainty quantification
on the quantity of interest estimate. This formulation has been introduced in
[13] without the full mathematical algorithm which is developed in the present
paper.
Let X be the random vector of influential and uncertain parameters of the
computer model G(·) and a be the parameter of interest (the abscissa of the
FRC). The prior knowledge on G(a, X) is modeled by Y (a, X) and defined as
follows .
Y (a, X) = β0 + β1 a + Zσ2 ,θ (a, X),

(6)

where Zσ2 ,θ is a centered Gaussian process. We make the assumption that
Zσ2 ,θ is second order stationary with variance σ 2 and a parametric covariance
function (which corresponds to the kernel in other statistical learning methods).
As covariance choice, one can cite the popular Matérn 5/2 parametrized by its
lengthscale θ). Thanks to the maximum likelihood method, we can estimate the
values of the so far-unknown parameters: β0 , β1 , σ 2 and θ (see for instance [7]
for more details).
7

Gaussian process regression (also known as the kriging process) provides an
estimator of G(a, x) which is called the kriging predictor and written Yb (a, x).
In addition to the kriging predictor, the kriging variance σY2 (a, x) quantifies
the uncertainty induced by estimating Y (a, x) with Yb (a, x). The predictive
distribution is given by Y (a, x) conditioned by yN :


∀x (Y (a, x) | yN ) = N Yb (a, x), σY2 (a, x)

(7)

where Yb (a, x) (the kriging mean) and σY2 (a, x) (the kriging variance) can both
be explicitly estimated.
Obtaining the FRC Ψ(a) consists in replacing Y = G(a, X) by its Gaussian
process metamodel (7) in (2). Hence we can estimate the value of Ψ(a), for
a > 0 from:
Ψ(a) = P (YN (a, X) > s|a) ,

(8)

where YN (a, X) is a Gaussian process which follows the distribution (Y (a, X) |
yN ) and s is a fixed threshold value. Two sources of uncertainty have to be
taken into account in (8): the first coming from the parameter X and the
second coming from the Gaussian distribution in (7).
From (8), the following estimate for Ψ(a) can be deduced:
"
!#
b (a, X)
s
−
Y
b
Ψ(a)
= EX 1 − Φ
.
σY (a, X)

(9)

This equation corresponds to the mean FRC with respect to the Gaussian process metamodel. The expectation in (9) is estimated using a classical Monte
Carlo integration procedure.
In order to estimate the uncertainty on the FRC estimation, we start from
its integral expression:
Z
Ψ(a) = PX (G(a, X) > s|a) =

1G(a,x)>s f (x)dx,

(10)

where f (x) is the joint pdf of X (independent on a). The first uncertainty
source on Ψ(a) comes from the numerical evaluation of the integral inside the
Ψ(a) definition. This evaluation is done by the Monte Carlo method:
n

Ψ(a) ≈ ΨM C (a) =

8

1X
1
,
(i)
n i=1 G(a,x )>s

(11)

where (x(i) )i=1,...,n is a sample of the random variable X. The required size of
this sample can be too large to use this formula in practice (case of a costly
evaluation of the code). We then replace the code by its Gaussian process
approximation YN (a, X):
n

ΨM C (a) ≈ ΨM C,GP (a) =

1X
1Y (a,X)>s .
n i=1 N

(12)

Let us recall that YN (a, X) is a Gaussian process with known mean and variance.
To compute the integration error, the central limit theorem is used and gives
for n → ∞:
√

n (ΨM C,GP (a) − ΨGP (a)) −→ N (0, ΨGP (a)(1 − ΨGP (a))) ,

(13)

where
Z
ΨGP (a) =

1YN (a,X)>s f (x)dx.

(14)

ΨGP (a) corresponds to the FRC for the process YN (a, X). The algorithm uses
a significantly large value of n in order to deal with a valid Gaussian approximation.
The second uncertainty source on Ψ(a) estimation comes from Gaussian
process approximation. We simulate m realizations (y (j) (a, x))j=1,...,m from
(Y (a, x)|yN ) in order to evaluate the variability of ΨM C,GP (a) which comes
from Gaussian process approximation. Thus, we compute:
n

(j)
ψM C,GP (a)

1X
=
1 (j) (i) .
n i=1 y (a,x )>s

(15)

(j)

Then, for each ψM C,GP (a), j = 1, . . . , m, we compute the Monte Carlo error
using the central limit theorem cited before in (13) where the variance of the
(j)

(j)

Gaussian distribution in (13) is estimated with ψM C,GP (a)(1 − ψM C,GP (a)),
j = 1, . . . , m.
From (13), we compute for each j = 1, . . . , m a sample of size nCLT . From
h
i
(j,k)
this double Monte Carlo method, we obtain a sample ψM C,GP (a) j=1,...,m
k=1,...,nCLT

of size m × nCLT . Then the FRC estimate is given by :
X
1
(j,k)
ψM C,GP (a),
m × nCLT j=1,...,m
k=1,...,nCLT

9

h
i
(j,k)
and the uncertainty is evaluate from ψM C,GP (a)

j=1,...,m

by using classical

k=1,...,nCLT

quantile estimators.
From this double Monte Carlo method, we are able to compute the FRC
uncertainty due to the Gaussian process metamodel and due to the numerical
integration process. Figure 1 illustrates the algorithm on a POD curve estimation in an NDT application (see section 2.1). A Gaussian process metamodel has
been built on N = 100 numerical simulations of eddy current examinations of
steam generator tubes [13]. We visualize the confidence interval induced by the
Monte Carlo (MC) estimation (n = 10, 000, nCLT = 100, 000), the one induced
by the Gaussian process approximation (m = 3, 000) and the total confidence
interval (including both approximations: GP+MC). It could be noticed that
h
i
(j,k)
from the double Monte-Carlo sample ψM C,GP (a) j=1,...,m , one can evaluk=1,...,nCLT

ate with an ANOVA decomposition the contribution of the Gaussian process
and the parameters X on the FRC estimate mean squared error (see [12]).

Figure 1: Example of FRC estimated with a Gaussian process model (from [13]). The present
curve is a POD curve. The lines represent unilateral confidence intervals (LB means “lower
bound”). The shaded areas represent the minimal and the maximal observations of POD(a)
when only the uncertainty on the Gaussian process is considered (light gray) and when both
the Gaussian process and the Monte-Carlo errors are considered (dark gray).

10

Choosing m, the number of conditional Gaussian process realizations to
simulate, is a difficult issue of this process. As explained in [12] and other
papers, there are some computational constraints (in memory and cpu time)
when simulating conditional Gaussian processes on large samples. The choice
of m = 3000 has been chosen as a sufficient value (that we avoid to be too large)
in order to estimate a 95% lower bound (due to the Gaussian process error).
Then n has been tuned in order to have a smaller error induced by the Monte
Carlo sampling process than the one induces by the Gaussian process.

4. Global sensitivity analysis of functional risk curves
The objective of sensitivity analysis is to determine those input variables
that mostly influence the model response [9, 29, 30]. Global sensitivity analysis methods take into account the overall uncertainty of the input parameters.
Previous works on the POD [31] and seismic fragility context [32] have only considered sensitivity analysis on the variance and distribution of the model output
variable Y . In this section, we propose two global sensitivity indices attached
to the whole FRC: one already introduced in [13] (aggregated Sobol’ indices)
and a new one based on a completely different idea (perturbed-law based indices). Therefore, the influence of the parameter of interest a, as it is the FRC
abscissa, will not be considered as an input variable whose sensitivity index has
to be computed.
4.1. Aggregated Sobol’ indices
We first focus on the variance-based sensitivity formulation [33], which is one
the most popular tools that has been proven to be robust and easily interpretable
[10]. The associated sensitivity indices are often called Sobol’ indices.
In the case of independent inputs, using the Sobol-Hoeffding decomposition
[34], the variance of the numerical model Y = G(X1 , . . . , Xd ) is decomposed
into the following sum:
Var(Y ) = V =

d
X
i=1

X

Vi +

1≤i<j≤d

11

Vij + . . . + V1...d

(16)

with Vi = Var[E(Y |Xi )], Vij = Var[E(Y |Xi , Xj )] − Vi − Vj , etc. Then, ∀i, j =
1 . . . d, i < j, the Sobol’ indices are given by
Si =

X
Vi
Vij
, Sij =
, . . . , and Ti = Si +
Sij +
V
V
j6=i

X

Sijk + . . . . (17)

j<k,j6=i,k6=i

The first-order Sobol’ index Si measures the individual effect of the input Xi
on the variance of the output Y , while the total Sobol’ index Ti measures the
Xi effect and all the interaction effects between Xi and the other inputs (as
V−i
with
the second-order effect Sij ). Ti can be reinterpreted as Ti = 1 −
V
V−i = Var[E(Y |X−i )] and X−i is the vector of all inputs except Xi . Each Sobol’
index is then interpreted in terms of percentage of total variance explanation.
In order to define similar sensitivity indices for the whole FRC (Eq. (9)),
[13] defines the following quantities:


b (a,X)
s−Y
σY (a,X)

ΨX (a)

=

PX (Y > s | a, X) = 1 − Φ

ΨXi (a)

=

PX (Y > s | a, Xi ) ,

ΨX−i (a)

=

PX (Y > s | a, X−i ) ,

D

=

EX kΨ(a) − ΨX (a)k2
hR
i
2
EX
(Ψ(a) − ΨX (a)) da ,

=



,

(18)

with k.k the euclidean norm. As in [35] (see also [36]) which deals with functional
model outputs, our objective is to obtain a single sensitivity index for each input
by averaging the variability of the function, here the FRC depending on a. The
specificity here is that we deal with a FRC function. Therefore, [13] defines the
FRC aggregated Sobol’ indices by:
SiFRC

=

FRC

=

Ti

EX kΨ(a) − ΨXi (a)k2
,
D
EX kΨX (a) − ΨX−i (a)k2
.
D

(19)

SiFRC (resp. TiFRC ) gives the first-order (resp. total) effect of Xi on the mean
FRC.
These FRC aggregated Sobol’ indices are easily computed by the metamodels. In particular, the
allows to replace PX (Y > s | a) by
" kriging metamodel!#
s − Yb (a, X)
the expectation EX 1 − Φ
in the FRC expressions of (18). In
σY (a, X)
12

order to compute the Sobol’ indices, we use the classical pick-freeze formulas
(see for example [34]). The confidence intervals due to the finite Monte Carlo
sample are then obtained by bootstrap. Let us noticed that it is also possible
to estimate the error on the aggregated Sobol’ indices due to the kriging error
as in [12] for the classical Sobol’ indices. This computationally heavy process
has not been applied here.
The sensitivity analysis results on our use case are given in Figure 2 with
the FRC aggregated Sobol’ indices of the five physical input parameters (called
E, ebav1 , ebav2 , h11 and h12 ). ebav1 is the main influent parameter on the FRC
(the POD curve of Figure 1), followed by h12 and ebav2 . E and h11 have no
influence.

0.8
0.6
0.4
0.2
0.0

0.0

0.2

0.4

T

S

0.6

0.8

1.0

Total effect sensitivity indices

1.0

Main effect sensitivity indices

E

ebav1

ebav2

h11

h12

E

ebav1

ebav2

h11

h12

Figure 2: First order (left) and total (right) FRC aggregated Sobol’ indices (from [13]).

The FRC aggregated Sobol’ indices quantify the sensitivity of each input on
the overall curve. However, we could be interested in the sensitivities on the
FRC at a specific value of the parameter of interest a. This can be directly
done by replacing Y by ΨX (a) in the equations (16) and (17). If we are now
interested by the sensitivities at a specific probability value of the FRC, we have
to study the inverse function of the FRC: Ψ−1
X (p) with p a given probability.
Similarly to the previous case, the Sobol’ indices can be obtained by replacing
13

Y by Ψ−1
X (p) in the equations (16) and (17). Figure 3 gives these sensitivity
indices on our data for p = 0.90 which corresponds to the quantile of a at order
90% (noted a90 ). We see that the influences of the inputs on a90 are very close
to the ones on the overall POD curve. It could be explained by the fact that
the different inputs have rather linear effects on the output of the model (the
amplitude).

0.8
0.6
S
0.4
0.2
0.0

0.0

0.2

0.4

S

0.6

0.8

1.0

a90 − Total effect sensitivity indices

1.0

a90 − Main effect sensitivity indices

E

ebav1

ebav2

h11

h12

E

ebav1

ebav2

h11

h12

Figure 3: First order (left) and total (right) Sobol’ indices on a90 (from [13]).

4.2. Perturbed-law based sensitivity indices
We propose now to quantify the impact on the FRC of a perturbation of the
input parameters pdfs by answering to the following question: “What would be
the FRC if the pdf of the ith input Xi had been modified?”. In this approach, all
the input parameters are modeled by random variables, and their input probability densities are supposed to be unknown. Let us remark that a negligible
sensitivity index of an input will not allow to fix this input, but just to say that
its pdf has no influence on the FRC.
We use the so-called Perturbed-Law based sensitivity Indices (PLI) measures
recently introduced in [14] (see also [37]). We start from the integral-form of

14

the FRC:
Z
PX (G(a, x) > s) =

1G(a,x)>s f (x) dx,

where f (x) is the joint pdf of X. Modifying the pdf fi (xi ) of Xi gives us fi,δ (xi ),
the perturbed pdf of Xi . After this perturbation, the FRC, denoted PXi ,δ (·)
instead of PX (·), can be written as:
Z
PXi ,δ (G(a, x) > s) =

1G(a,x)>s

fi,δ (xi )
f (x) dx,
fi (xi )

(20)

The PLI measures only consist of the comparison of the FRCs before and
after the perturbation, and are defined by:

PXi ,δ (G(a, X) > s) − PX (G(a, X) > s)





PXi ,δ (G(a, X) > s)





if
P
(G(a,
X) > s) ≥ PX (G(a, X) > s) ,
Xi ,δ


Si,δ =


 PXi ,δ (G(a, X) > s) − PX (G(a, X) > s)





PX (G(a, X) > s)



 if PX ,δ (G(a, X) > s) < PX (G(a, X) > s) .

(21)

i

A negative Si,δ means that the FRC is smaller after the perturbation, while a
positive Si,δ means that the FRC has increased. The estimation of PXi ,δ (G(a, X) > s)
is based on reverse importance sampling [38]. Asymptotical properties of the
estimators give also some confidence intervals on the PLI measures.
In [14], the numerical model is directly used to estimate the PLI measures
by large Monte Carlo samples. We propose here to estimate the PLI measures
using the Gaussian process metamodel and integrating its error in the estimates.
The mean FRC which we consider is given by (9):
"
!# Z
!
s − Yb (a, X)
s − Yb (a, x)
Ψ(a) = EX 1 − Φ
= 1−Φ
f (x) dx. (22)
σY (a, X)
σY (a, x)
After the perturbation of the pdf of the ith input, the mean FRC is given by:
!
Z
s − Yb (a, x) fi,δ (xi )
Ψi,δ (a) = 1 − Φ
f (x) dx.
(23)
σY (a, x)
fi (xi )

15

PLI measures are then given by:

Ψi,δ (a) − Ψ(a)





Ψi,δ (a)

Si,δ (a) =



Ψi,δ (a) − Ψ(a)



Ψ(a)

if

Ψi,δ (a) ≥ Ψ(a),
(24)

if

Ψi,δ (a) < Ψ(a).

The last element of the PLI method is the definition of the perturbations
fi,δ (xi ). [14] choose to perturb a statistical characteristic (for example the mean,
or the variance, or a quantile, . . . ) of Xi in order to be “as close as possible”
to the initial pdf fi (xi ). The dissimilarity measure between fi,δ (xi ) and fi (xi ),
which contains the required properties, is the Kullback-Leibler divergence:
Z +∞
fi,δ (xi )
KL(fi,δ (xi ), fi (xi )) =
fi,δ (xi ) log
dxi .
(25)
fi (xi )
−∞
It implies that fi,δ (Xi ) and fi (X) have the same definition domain (as a consequence, the input domain bounds cannot be changed). The Kullback-Leibler
divergence has the advantage to be easily minimized [39], in order to obtain
explicit solutions for fi,δ (Xi ) for classical pdf (e.g. Gaussian).
Figure 4 gives some examples of perturbations for the uniform pdf on [0, 1]
(mean is 0.5 and variance is 1/12). Figure 5 gives some examples of perturbations for the standard Gaussian pdf (mean is 0 and variance is 1).
Figure 6 gives an example of a PLI-based sensitivity analysis on our NDT
test case (see section 2.1) which aims to estimate POD curves and the associated
sensitivity indices to its five physical input parameters (called E, ebav1 , ebav2 ,
h11 and h12 ), the defect size a being fixed here at a given value. The input
parameters pdfs are all uniform on [0, 1]. Their non-perturbed mean is then
1/2. Then, the mean of each input is modified in the range of δ ∈ [0.1, 0.9]. For
each δ value and each influential parameter, Figure 6 shows the PLI estimates
Ŝδ . Large absolute values of Ŝδ imply a large impact of the perturbation on the
FRC. Moreover, the sign of Ŝδ indicates if the new probability has decreased
(negative case) or increased (positive case). Confidence intervals shown in Figure
6 are obtained using a Gaussian asymptotic property of the estimates Ŝi,δ [14].
In Figure 6, one can see for instance that increasing the h12 input mean largely
16

Variance shifting for the Uniform distribution
3.0

10

Mean shifting for the Uniform distribution

2.0

2.5

δ = 0.5 12
δ = 0.75 12
δ = 1.25 12
δ = 1.5 12
δ = 1 12

1.5

fδ(x)
0

0.0

0.5

2

1.0

4

fδ(x)

6

8

δ = 0.1
δ = 0.2
δ = 0.3
δ = 0.4
δ = 0.5

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

x

0.4

0.6

0.8

1.0

x

Figure 4: Examples of perturbations for the uniform pdf on [0, 1]. Left: the mean of the pdf
is perturbed (δ is the new mean value). Right: the variance of the pdf is perturbed (δ is the
new variance value).

0.8
0.6

δ = 0.5
δ = 0.75
δ = 1.25
δ = 1.5
δ=1

0.2

0.2

fδ(x)

0.3

δ = − 0.5
δ = − 0.25
δ = 0.25
δ = 0.5
δ=0

0.0

0.0

0.1

fδ(x)

Variance shifting for the Gaussian distribution

0.4

0.4

Mean shifting for the Gaussian distribution

−4

−2

0

2

4

−4

x

−2

0

2

4

x

Figure 5: Examples of perturbations for the standard Gaussian pdf. Left: the mean of the pdf
is perturbed (δ is the new mean value). Right: the standard-deviation of the pdf is perturbed
(δ is the new standard-deviation value).

increases the FRC, while decreasing the ebav1 input mean largely decreases the
FRC. In contrary, perturbations of the E mean have no effect on the FRC.
In order to have a global view of the FRC sensitivity for different defect sizes
a, Figure 7 gives the PLI graphs for each input, with variations of a and the
δ perturbation. Clearly, the inputs ebav1 and h12 have strong impacts on the

17

Figure 6: Graphical example for the PLI measures on the POD curves. Each color corresponds
to one input parameter. The dots show the PLI values in function of the perturbation δ on
the input pdf mean, while the dotted curves correspond to the 95% confidence intervals on
the PLI estimates.

FRC, but especially in the a < 0.3 range. In contrary, the inputs E and h11 have
no influence. These elements are key points for the engineers to understand the
effects of the physical parameters on the POD curve.

5. Conclusion
In this paper, several methods of uncertainty and sensitivity analysis of
model outputs have been developed on the particular quantity of interest of the
functional risk curve (FRC). A FRC provides the probability of an undesirable
event as a function of the value of a critical input parameter of a considered
physical system. Focus has been on the use of the Gaussian process metamodels
in order to build FRCs from numerical experiments. This approach is useful
when the computer model is expensive to evaluate in such a way that only a
small sample of the model output can be obtained. In addition to the mean

18

^
Sδ − parameter : ebav1
0.3

0.02

0.2

0.3
0.2

8

0.8

0

2

0.1

2
0.1

0.8

−0.02

0.1

^
Sδ − parameter : E

0.1
6

0.2

0.2

0.14
8
0.0
6

0.6

0.1
0

0.0

0.4

2

0.1

0.0

0.04
0

δ

δ

0.6

0.0

−0.0

0.0

2

−0.04

0.4

−0.08

0.2

4

−0.1

−0.12

0.0

6

−0.14
0.0

8

0.2

0

2

0.1
2

0.0

0

0.1

.06

−0

4

0.0

0.06

−0

2

.1

0.0

−0.1

0.0

0.2

0

−0.2
0.10

0.15

0.20

0.25

0.30

0.35

−0.2

0.40

0.10

0.15

0.20

0.25

a

0.30

0.35

0.40

a

^
Sδ − parameter : ebav2

^
Sδ − parameter : h11
0.3

−0.04

0.3

0.8

0.04

0.8
0.2

0.0

0.2

2

−0.02

0.6

0.1
0

0

δ

δ

0.6

0.0

0.4

0.1
0

0.0

0.4

0.02

−0.02
0.04

−0.1

−0.02

−0.1

0.06

0.2

0.2
−0.04

−0.04

−0.06

0.10

0.15

0.20

0.25

−0.2

0.30

0.35

−0.2

0.40

0.10

0.15

0.20

0.25

a

0.30

0.35

0.40

a

^
Sδ − parameter : h12
0.3

0.1

4

0.2

0.22

8

2

0.1

0.2

0.8

0.0

0.16

0.14

0.2

0.12
6
0.0

0.08

δ

0.6

0.1

0.04

0

−0.02

0.0

−0.04

0.4

−0.06
−0.08

−0

.1

−0.1
−0.12

0.2

−0.1

4

−0.2

0

0.10

0.15

0.20

0.25

0.30

0.35

0.40

a

Figure 7: PLI measures on the POD curves for different perturbations on the pdf means
and different values of a. δ = 0.5 corresponds to the non perturbed case. The grey color
corresponds to no change in the FRC, while the red color (resp. blue) corresponds to an
increase (resp. decrease) of the FRC.

risk curve, the metamodel allows to obtain the confidence bands via conditional
Gaussian process simulations.

19

One important advantage of using a metamodel is that sensitivity analysis is
facilitated. We have defined two kinds of sensitivity indices related to the FRC
as the quantity of interest. First, we have formulated the FRC aggregated Sobol’
indices, which are variance-based measures. Second, based on perturbation of
the pdf of each model input variable, sensitivity indices (called PLI) of the model
inputs on the FRC are also proposed. These allow to understand the effect of
misjudgment on the pdf of each input parameters.
An example taken from simulated NDE inspections highlights the added
value of the FRC, which in this context correspond to POD curves. The two
proposed sensitivity analysis methods provide strongly complementary information to experts interpreting the NDE results. However, FRCs are used in many
other engineering contexts, e.g. in the seismic fragility assessment (as shown
in section 2.2) and in the evaluation of hydraulic works reliability subject to
extremely high water levels [2]. These tools would also be useful in separating
the effects of stochastic input variables and epistemic input parameters. For
instance, this problem is highlighted in [40] and [41].
Finally, two main mathematical perspectives are identified from this work.
First, adaptive designs can be developed: in the spirit of [42], it would consist
of using the Gaussian process model in order to define some SUR (“Sequential
Uncertainty Reduction”) criteria on the FRC as the quantity of interest. Optimizing such a criterion would provide a new set of input parameter values, which
would be run with the computer code in order to decrease the POD confidence
interval. Second, [43] have started a preliminary study that considers FRC as
a random distribution function. This framework allows to deal with stochastic
computer codes instead of the deterministic ones of this paper (see for example
Kleijnen [44]).

Acknowledgments
Part of this work has been backed by French National Research Agency
(ANR) through project ByPASS ANR-13-MONU-0011. We acknowledge three

20

anonymous reviewers for their deep analysis of the paper and the useful comments they provide. We are also grateful to Alberto Pasanisi for giving helpful
discussions and Andreas Schumm for his help with the English language.

References
References
[1] E. de Rocquigny, Modelling Under Risk and Uncertainty: An Introduction
to Statistical, Phenomenological and Computational Methods, Wiley, 2012.
[2] A. Pasanisi, Uncertainty analysis and decision-aid: methodological, technical and managerial contributions to engineering and R&D studies, Mémoire
d’Habilitation à Diriger des Recherches de l’Université de Technologie de
Compiègne, 2014.
[3] M. Baudin, A. Dutfoy, B. Iooss, A. Popelin, Open TURNS: An industrial software for uncertainty quantification in simulation, in: R. Ghanem,
D. Higdon, H. Owhadi (Eds.), Springer Handbook on Uncertainty Quantification, Springer, 2017, pp. 2001–2038.
[4] K.-T. Fang, R. Li, A. Sudjianto, Design and modeling for computer experiments, Chapman & Hall/CRC, 2006.
[5] R. Smith, Uncertainty quantification, SIAM, 2014.
[6] C. Rasmussen, C. Williams, Gaussian processes for machine learning, MIT
Press, 2006.
[7] A. Marrel, B. Iooss, F. Van Dorpe, E. Volkova, An efficient methodology for
modeling complex computer codes with Gaussian processes, Computational
Statistics and Data Analysis 52 (2008) 4731–4744.
[8] L. Le Gratiet, S. Marelli, B. Sudret, Metamodel-based sensitivity analysis: Polynomial chaos expansions and Gaussian processes, in: R. Ghanem,

21

D. Higdon, H. Owhadi (Eds.), Springer Handbook on Uncertainty Quantification, Springer, 2017, pp. 1289–1325.
[9] A. Saltelli, K. Chan, E. Scott (Eds.), Sensitivity analysis, Wiley Series in
Probability and Statistics, Wiley, 2000.
[10] B. Iooss, P. Lemaı̂tre, A review on global sensitivity analysis methods,
in: C. Meloni, G. Dellino (Eds.), Uncertainty management in SimulationOptimization of Complex Systems: Algorithms and Applications, Springer,
2015, pp. 101–122.
[11] J. Oakley, A. O’Hagan, Probabilistic sensitivity analysis of complex models:
A Bayesian approach, Journal of the Royal Statistical Society, Series B 66
(2004) 751–769.
[12] L. Le Gratiet, C. Cannamela, B. Iooss, A Bayesian approach for global
sensitivity analysis of (multifidelity) computer codes, SIAM/ASA Journal
of Uncertainty Quantification 2 (2014) 336–363.
[13] L. Le Gratiet, B. Iooss, T. Browne, G. Blatman, S. Cordeiro, B. Goursaud,
Model assisted probability of detection curves: New statistical tools and
progressive methodology, Journal of Nondestructive Evaluation 36 (2017)
8.
[14] P. Lemaı̂tre, E. Sergienko, A. Arnaud, N. Bousquet, F. Gamboa, B. Iooss,
Density modification based reliability sensitivity analysis, Journal of Statistical Computation and Simulation 85 (2015) 1200–1223.
[15] L. Maurice, V. Costan, E. Guillot, P. Thomas, Eddy current NDE performance demonstrations using simulation tools, AIP Conference Proceedings
1511 (2013) 464–471.
[16] T. Browne, L. Le Gratiet, G. Blatman, S. Cordeiro, B. Goursaud, B. Iooss,
L. Maurice, Building Probability of Detection curves via metamodels, in:

22

12th International Conference on Applications of Statistics and Probability in Civil Engineering, ICASP12, Vancouver, Canada, The University of
British Columbia, 2015, p. 573.
[17] R. Thompson, A unified approach to the Model-Assisted determination of
Probability of Detection, in: 34th Annual Review of Progress in Quantitative Nondestructive evaluation, July 2007, Golden, Colorado, American
Institute of Physics, Melville, New-York, 2008, pp. 1685–1692.
[18] P. Calmon, Trends and stakes of NDT simulation, Journal of Non Destructive Evaluation 31 (2012) 339–341.
[19] L. Maurice, V. Costan, P. Thomas, Axial probe eddy current inspection
of steam generator tubes near anti-vibration bars: performance evaluation
using finite element modeling, Proceedings of JRC-NDE, Cannes, France
(2013) 638–644.
[20] I. Zentner, Numerical computation of fragility curves for NPP equipment,
Nuclear Engineering and Design 240 (2010) 1614–1621.
[21] I. Zentner, E. Borgonovo, Construction of variance-based metamodels for
probabilistic seismic analysis and fragility assessment, Georisk 8 (2014)
202–216.
[22] G. Damblin, M. Keller, A. Pasanisi, P. Barbillon, E. Parent, Approche
décisionnelle bayésienne pour estimer une courbe de fragilité, Journal de la
Société Française de Statistique 155 (2014) 78–103.
[23] M. Shinozuka, M. Feng, J. Lee, T. Naganuma, Statistical analysis of
fragility curves, Journal of Engineering Mechanics 126 (2000) 1224–1231.
[24] A. Berens, NDE reliability data analysis, in: Metals Handbook, 9th edition,
Vol. 17, 1988, pp. 689–701.
[25] G. Box, D. Cox, An analysis of transformations, Journal of the Royal Statistical Society 26 (1964) 211–252.
23

[26] C. Mai, K. Konakli, B. Sudret, Seismic fragility curves for structures using
non-parametric representations, Frontiers of Structural and Civil Engineering 11 (2017) 169–186.
[27] S. Demeyer, F. Jenson, N. Dominguez, Modélisation d’un code numérique
par un processus gaussien - Application au calcul d’une courbe de probabilité de dépasser un seuil, in: Proceedings of 44èmes Journées de Statistique, Bruxelles, Delgique, 2012.
[28] J. Sacks, W. Welch, T. Mitchell, H. Wynn, Design and analysis of computer
experiments, Statistical Science 4 (1989) 409–435.
[29] E. Borgonovo, E. Plischke, Sensitivity analysis: A review of recent advances, European Journal of Operational Research 248 (2016) 869–887.
[30] B. Iooss, A. Saltelli, Introduction: Sensitivity analysis, in: R. Ghanem,
D. Higdon, H. Owhadi (Eds.), Springer Handbook on Uncertainty Quantification, Springer, 2017, pp. 1103–1122.
[31] F. Rupin, G. Blatman, S. Lacaze, T. Fouquet, B. Chassignole, Probabilistic
approaches to compute uncertainty intervals and sensitivity factors of ultrasonic simulations of a weld inspection, Ultrasonics 54 (2014) 1037–1046.
[32] E. Borgonovo, I. Zentner, A. Pellegri, S. Tarantola, E. de Rocquigny, On the
importance of uncertain factors in seismic fragility assessment, Reliability
Engineering and System Safety 109 (2013) 66–76.
[33] I. Sobol, Sensitivity estimates for non linear mathematical models, Mathematical Modelling and Computational Experiments 1 (1993) 407–414.
[34] C. Prieur, S. Tarantola, Variance-based sensitivity analysis: Theory and
estimation algorithms, in: R. Ghanem, D. Higdon, H. Owhadi (Eds.),
Springer Handbook on Uncertainty Quantification, Springer, 2017, pp.
1217–1239.

24

[35] F. Gamboa, A. Janon, T. Klein, A. Lagnoux, Sensitivity analysis for multidimensional and functional outputs, Electronic Journal of Statistics 8 (1)
(2014) 575–603.
[36] A. Marrel, N. Saint-Geours, Sensitivity analysis of spatial and/or temporal phenomena, in: R. Ghanem, D. Higdon, H. Owhadi (Eds.), Springer
Handbook on Uncertainty Quantification, Springer, 2017, pp. 1327–1357.
[37] E. Borgonovo, B. Iooss, Moment-independent and reliability-based importance measures, in: R. Ghanem, D. Higdon, H. Owhadi (Eds.), Springer
Handbook on Uncertainty Quantification, Springer, 2017, pp. 1265–1287.
[38] T. Hesterberg, Estimates and confidence intervals for importance sampling
sensitivity analysis, Mathl. Comput. Modelling 23 (1996) 79–85.
[39] P. Lemaı̂tre, Analyse de sensibilité en fiabilité des structures - Sensitivity
analysis in structural reliability, Thèse de l’Université Bordeaux I, 2014.
[40] J. Helton, J. Johnson, W. Obekampf, C. Salaberry, Representation of analysis results involving aleatory and epistemic uncertainty, International Journal of General Systems 39 (2010) 605–646.
[41] A. Pasanisi, M. Keller, E. Parent, Estimation of a quantity of interest in
uncertainty analysis: Some help from Bayesian decision theory, Reliability
Engineering and System Safety 100 (2012) 93–101.
[42] J. Bect, D. Ginsbourger, L. Li, V. Picheny, E. Vazquez, Sequential design of computer experiments for the estimation of a probability of failure,
Statistics and Computing 22 (2012) 773–793.
[43] T. Browne, J.-C. Fort, Redéfinition de la POD comme fonction de
répartition aléatoire, in: Proceedings of 47èmes Journées de Statistique,
Lille, France, 2015.
[44] J. Kleijnen, Design and analysis of simulation experiments, Second Edition,
Springer, 2015.
25

