# Local Sensitivity of Failure Probability through Polynomial Regression and Importance Sampling

Marie Chiron, Jérôme Morio, Sylvain Dubreuil

► To cite this version:

Marie Chiron, Jérôme Morio, Sylvain Dubreuil. Local Sensitivity of Failure Probability through Polynomial Regression and Importance Sampling. Mathematics, 2023, 11 (20), pp.4357. 10.3390/math11204357. hal-04628868

HAL Id: hal-04628868

https://hal.science/hal-04628868v1

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Local Sensitivity of Failure Probability through Polynomial Regression and Importance Sampling

Marie Chiron, Jérôme Morio* and Sylvain Dubreuil

ONERA/DTIS, Université de Toulouse, F-31055 Toulouse, France; marie.chiron@onera.fr (M.C.); sylvain.dubreuil@onera.fr (S.D.)

* Correspondence: jerome.morio@onera.fr

**Abstract:** Evaluating the failure probability of a system is essential in order to assess its reliability. This probability may significantly depend on deterministic parameters such as distribution parameters or design parameters. The sensitivity of the failure probability with regard to these parameters is then critical for the reliability analysis of the system or in reliability-based design optimization. Here, we introduce a new approach to estimate the failure probability derivatives with respect to deterministic inputs, where the bias can be controlled and the simulation budget is kept low. The sensitivity estimate is obtained as a byproduct of a heteroscedastic polynomial regression with a database built with simulation methods. The polynomial comes from a Taylor series expansion of the approximated sensitivity domain integral obtained with the Weak approach. This new methodology is applied to two engineering use cases with the importance sampling strategy.

**Keywords:** failure probability; reliability-based sensitivity analysis; local sensitivity; heteroscedastic polynomial regression

**MSC:** 65C05; 26A24

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Chiron, M.; Morio, J.; Dubreuil, S. Local Sensitivity of Failure Probability through Polynomial Regression and Importance Sampling. Mathematics 2023, 11, 4357. https://doi.org/10.3390/math11204357

Academic Editors: Sinan Xiao and Pan Wang

Received: 11 September 2023

Revised: 13 October 2023

Accepted: 15 October 2023

Published: 20 October 2023

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

In many scientific fields, a complex system is often modeled with a function $g$ expected to simulate the behavior of the system. The output $Y = g(s, X)$ of this function is the observed response, where $s \in \mathbb{R}^p$ and $X \in \mathbb{R}^d$ denote, respectively, the deterministic inputs and the random inputs of $g$. The deterministic inputs $s$ can either be distribution parameters of the random variable $X$ or design parameters of $g$. Function $g$ is not analytically known but can be called for any input values $(s, X)$ where $g$ is defined. The failure of the system is then characterized by a negative value of the observed response, $Y$, and $g$ is referred to as the limit state function (lsf). The estimation of the failure probability $P_f = \mathbb{P}(Y < 0)$ provides crucial information about the reliability of the system but it greatly depends on the settings of the different inputs. The study of the influence of the inputs of $g$ on $P_f$ is the purpose of reliability-based sensitivity analysis (RSA) [1].

In this article, we uniquely focus on local RSA. The local method consists of computing the derivatives of $P_f(s)$ according to $s$ for a fixed value of $s$. The local sensitivity of $P_f$, with regard to the distribution parameters, provides valuable information about the influence of the probabilistic model selected for the random inputs on the system's failure occurrence [2]. Whereas the local sensitivity of $P_f$ with regard to the design parameters can be of great use for reliability-based design optimization (RBDO) [3].

The derivatives of $P_f$ with regard to the design parameters of $g$ necessarily lead to surface integrals [4]. However, the derivatives of $P_f$, with regard to distribution parameters of $X$, are either domain integrals in the original input space [5] or surface integrals if a standardized input space is considered. Domain integrals are much easier to handle than surface integrals; nevertheless, the literature on failure probability estimation in
standardized spaces is more luxuriant than in the original input space [6]. The main challenge of the estimation of the derivatives with respect to s, no matter its nature, is to increase as little as possible the simulation budget required for the estimation of the probability of failure $P_f(s)$. Consequently, it is assumed that the computation of this sensitivity reuses as much as possible the evaluations of the limit state function needed for the computation of $P_f$. In other words, it is uncommon to estimate $P_f$ and its derivatives in different input spaces, as the limit state functions differ. In this paper, a new approach to compute the failure probability sensitivity in the standard normal space is presented. This new approach relies on simulation methods and a heteroscedastic polynomial regression. In order to reuse the evaluations of the lsf g, in this article, we only focus on estimating the failure probability $P_f(s)$ using importance sampling [7] and the derivatives of $P_f(s)$, with respect to any deterministic inputs s.

The article is organized as follows. In Section 2, the local RSA background is detailed with an emphasis on the Weak approach [8], as the proposed approach is greatly inspired by it. Sections 3 and 4 describe the proposed algorithm. Section 3 focuses on the derivation of the new expression of the sensitivities, while Section 4 presents the heteroscedastic polynomial regression. Section 5 presents two numerical examples, with standard normal inputs. The results are summarized and conclusions are drawn in Section 6.

## 2. Local RSA with Respect to Deterministic Inputs

### 2.1. Integral Expression of the Local RSA

In the standard normal space, the dependence of $P_f$ on $s = [s_1, \dots, s_p]$ is necessarily contained in the limit state function g, while the standardized joint distribution $f_X$ is parameter-free. The failure probability is, thus, written as follows

$$P_f(s) = \int_{D_f(s)} f_X(\mathbf{x}) \, d\mathbf{x} = \int_{\mathbb{R}^d} \mathbb{I}_{D_f(s)}(\mathbf{x}) f_X(\mathbf{x}) \, d\mathbf{x}, \quad (1)$$

where $D_f(s) = \{\mathbf{x} \in \mathbb{R}^d | g(\mathbf{s}, \mathbf{x}) \le 0\}$ is the failure domain in the standard space and $\mathbb{I}$ is an indicator function. Introducing an auxiliary density, $h$, the failure probability $P_f(s)$ can then be estimated by the importance sampling (IS) method [7] leading to the following estimate

$$\widehat{P_f(s)}^{\text{IS}} = \frac{1}{N} \sum_{j=1}^{N} \mathbb{I}_{D_f(s)}(\mathbf{x}^{(j)}) \frac{f_X(\mathbf{x}^{(j)})}{h(\mathbf{x}^{(j)})}, \quad (2)$$

where $\mathbf{x}^{(j)}$ are independent and identically distributed (iid) from $h$. The classical Monte Carlo method [9] is a special case of IS, where $h = f_X$. The main advantage of IS compared to Monte Carlo is an improved convergence speed, particularly for rare event estimation when $h$ is well chosen.

Assuming the gradients $\nabla_X g(\mathbf{s}, \mathbf{x}) \neq 0$ for all $\mathbf{x}$ and $\mathbf{s}$ on the limit state surface $\{g(\mathbf{s}, \mathbf{x}) = 0\}$, the derivatives of $P_f(s)$, with respect to $s_\ell \in \mathbb{R}$ for $\ell \in [1, \dots, p]$, are defined by the following surface integral [10]

$$\frac{\partial P_f(s)}{\partial s_\ell} = - \int_{g(s, x)=0} \frac{1}{\| \nabla_X g(s, x) \|} \frac{\partial g(s, x)}{\partial s_\ell} f_X(\mathbf{x}) d\mathbf{c}(\mathbf{x}). \quad (3)$$

where $d\mathbf{c}(\mathbf{x})$ stands for surface integration over the limit state surface $\{g(s, \mathbf{x}) = 0\}$. The resulting expression of the derivative of $P_f$ in Equation (3) is, thus, a surface integral, which depends on the derivatives of g with respect to $s_\ell$ and its gradient $\nabla_X g$.

### 2.2. Weak Approach: Approximation of the Indicator Function

The failure domain indicator function makes the direct differentiation of the failure probability integral not possible. By replacing the failure domain indicator function with a smoother function, the differentiation can be performed and results in a domain inte-
gral [11]; it is the Weak approach. Several smooth approximations of the indicator function have been derived in the literature [12], which are typically cumulative distribution functions (cdfs) of continuous univariate variables [8]. Here, we derive the approximation chosen in [4,11], which comes from the following limit

$$ \forall \mathbf{x} \in \mathbb{R}^d \quad \mathbb{I}_{D_f}(\mathbf{x}) = \lim_{\sigma \to 0} \Psi\left(-\frac{g(\mathbf{s}, \mathbf{x})}{\sigma}\right), \quad (4) $$

where $\Psi$ is the standard normal univariate cdf, and $\sigma > 0$.

In the sense of distributions, the derivative of $P_f$, with respect to $s_\ell$, can then be defined as in [8]

$$ \begin{aligned} \frac{\partial P_f(\mathbf{s})}{\partial s_\ell} &= \lim_{\sigma \to 0} \int_{\mathbb{R}^d} \frac{\partial \Psi(-g(\mathbf{s}, \mathbf{x})/\sigma)}{\partial s_\ell} f_X(\mathbf{x}) \mathrm{d}\mathbf{x} \quad (5) \\ &= -\lim_{\sigma \to 0} \int_{\mathbb{R}^d} \frac{1}{\sigma} \frac{\partial g(\mathbf{s}, \mathbf{x})}{\partial s_\ell} \phi\left(-\frac{g(\mathbf{s}, \mathbf{x})}{\sigma}\right) f_X(\mathbf{x}) \mathrm{d}\mathbf{x}, \end{aligned} $$

where $\phi$ is the univariate standard normal probability density function (pdf). Consequently, an approximation of the failure probability sensitivity can be obtained with the following domain integral

$$ \frac{\partial P_f(\mathbf{s})}{\partial s_\ell} \approx - \int_{\mathbb{R}^d} \frac{1}{\tilde{\sigma}} \frac{\partial g(\mathbf{s}, \mathbf{x})}{\partial s_\ell} \phi\left(-\frac{g(\mathbf{s}, \mathbf{x})}{\tilde{\sigma}}\right) f_X(\mathbf{x}) \mathrm{d}\mathbf{x} = \frac{\partial \tilde{P}_f(\mathbf{s}, \tilde{\sigma})}{\partial s_\ell} \quad (6) $$

$$ \text{where} \quad \tilde{P}_f(\mathbf{s}, \tilde{\sigma}) = \int_{\mathbb{R}^d} \Psi\left(-\frac{g(\mathbf{s}, \mathbf{x})}{\tilde{\sigma}}\right) f_X(\mathbf{x}) \mathrm{d}\mathbf{x}, \quad (7) $$

for $\tilde{\sigma}$ a fixed positive value. The domain integral $\partial \tilde{P}_f(\mathbf{s}, \tilde{\sigma}) / \partial s_\ell$ can then be estimated with IS with the following estimate

$$ \frac{\partial \widetilde{P}_f(\mathbf{s}, \tilde{\sigma})}{\partial s_\ell} = - \frac{1}{N} \sum_{j=1}^{N} \frac{1}{\tilde{\sigma}} \frac{\partial g(\mathbf{s}, \mathbf{X}^{(j)})}{\partial s_\ell} \phi\left(-\frac{g(\mathbf{s}, \mathbf{X}^{(j)})}{\tilde{\sigma}}\right) \frac{f_X(\mathbf{X}^{(j)})}{h(\mathbf{X}^{(j)})}, \quad (8) $$

where $\mathbf{X}^{(j)}$ are iid from $h$. The same sample can be reused to compute both the failure probability estimate with IS and this sensitivity estimate. However, for each observation $\mathbf{X}^{(j)}$, the derivative of the lsf, with respect to $s_\ell$, has to be evaluated, which may increase the simulation budget.

This estimate is biased since $\tilde{\sigma} \neq 0$. The choice of $\tilde{\sigma}$ is crucial and has influence on both the bias and the variance of the estimate [8]. Decreasing the parameter $\tilde{\sigma}$ greatly reduces the bias. Nevertheless, past a certain point, the variance of the estimate increases for smaller $\tilde{\sigma}$ values, if $N$ is kept constant.

The Weak approach has been applied with various advanced sampling techniques. It has notably been associated with the sequential importance sampling (SIS) framework in [4] and the subset sampling framework in [12], to obtain more efficient sensitivity estimates.

### 3. Sensitivity with Respect to Deterministic Inputs through Taylor Series Expansion

The proposed method is greatly inspired by the Weak approach, as an approximation of the failure domain indicator function is also employed. The main idea of the proposed approach is to use sampling methods to compute the derivative of the failure probability as a byproduct of a heteroscedastic polynomial regression. In doing so, the bias included in the Weak approach can be controlled.

The first step is, thus, to approximate the failure domain indicator function with a smoother function. Then, a random variable change in the image measure is introduced. The resulting integral is differentiable, and its derivative, with respect to $s_\ell$ with $\ell = 1, \dots, p$,
is expressed as an expected value. Next, the Taylor series expansion of the expected value is derived and the sensitivity with respect to $s_\ell$ is finally identified amongst the Taylor series coefficients. These different steps are detailed in the following sections.

### 3.1. Approximation of the Failure Indicator Function

As the indicator function is not differentiable, several smoother functions have been used in the literature as surrogates, as mentioned before [8,11,12]. Here, we focus on approximations which are continuous cumulative distribution functions $\Xi_\sigma$, defined with a parameter $\sigma > 0$, which verify the following property

$$ \forall \mathbf{x} \in \mathbb{R}^d \quad \mathbb{I}_{D_f(\mathbf{s})}(\mathbf{x}) = \mathbb{I}_{y \le 0}(g(\mathbf{s}, \mathbf{x})) = \lim_{\sigma \to 0} \Xi_\sigma(-g(\mathbf{s}, \mathbf{x})), \quad (9) $$

where 0 is contained in the interior of the support of $\Xi_\sigma$. A typical example is presented in Equation (4) with the standard normal cdf. The failure probability function $\overline{P}_f(\mathbf{s}, \cdot)$ is then defined on $\mathbb{R}^+ \setminus \{0\}$, such as

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \overline{P_f(\mathbf{s}, \sigma)} = \int_{\mathbb{R}^d} \Xi_\sigma(-g(\mathbf{s}, \mathbf{x})) f_X(\mathbf{x}) d\mathbf{x} = \mathbb{E}_{f_X}[\Xi_\sigma(-g(\mathbf{s}, \mathbf{X}))], \quad (10) $$

and we further assume the cdf $\Xi_\sigma$ to be regular enough to verify the following properties:

$$ P_f(\mathbf{s}) = \lim_{\sigma \to 0} \overline{P_f(\mathbf{s}, \sigma)} \quad \text{and} \quad \frac{\partial P_f(\mathbf{s})}{\partial s_\ell} = \lim_{\sigma \to 0} \frac{\partial \overline{P_f(\mathbf{s}, \sigma)}}{\partial s_\ell}. \quad (11) $$

In the Weak approach framework presented in Section 2.2, the sensitivity of $P_f$ is obtained by differentiating Equation (10) with respect to $s_\ell$, $\ell = 1, \dots, p$, for a fixed value of $\sigma$, denoted as $\tilde{\sigma}$. Here, we introduce a modification to the approach with a random variable change to derive another expression of the failure probability function of Equation (10).

### 3.2. Change of Variable

Throughout the rest of this paper, it is assumed that the random response of the system $g(\mathbf{s}, \mathbf{X})$ is an absolutely continuous univariate random variable. Considering the following random variable change in the image measure: $H_s = g(\mathbf{s}, \mathbf{X})$, the failure probability function Equation (10) becomes

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \overline{P_f(\mathbf{s}, \sigma)} = \mathbb{E}_{f_{H_s}}[\Xi_\sigma(-H_s)] = \int_{\mathbb{R}} \Xi_\sigma(-h) f_{H_s}(h) dh, \quad (12) $$

where $f_{H_s}$ is the unknown density of the univariate random variable $H_s$. It should be noted that the failure probability function is, thus, written as an integral defined over $\mathbb{R}$ rather than $\mathbb{R}^d$. Denoting $\tilde{\xi}_\sigma$, the pdf associated with the $\Xi_\sigma$ distribution, the failure probability function of Equation (12) can then be rewritten in the following way

$$ \begin{aligned} \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \overline{P_f(\mathbf{s}, \sigma)} &= \int_{\mathbb{R}} \Xi_\sigma(-h) f_{H_s}(h) dh = \int_{\mathbb{R}} \left( \int_{-\infty}^{-h} \tilde{\xi}_\sigma(z) dz \right) f_{H_s}(h) dh \quad (13) \\ &= \int_{\mathbb{R}} \left( \int_{\mathbb{R}} \mathbb{I}_{\{z \le -h\}} \tilde{\xi}_\sigma(z) dz \right) f_{H_s}(h) dh. \end{aligned} $$

Applying the theorem of Fubini to integral Equation (13) results in

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \overline{P_f(\mathbf{s}, \sigma)} = \int_{\mathbb{R}} \left( \int_{\mathbb{R}} \mathbb{I}_{\{h \le -z\}} f_{H_s}(h) dh \right) \tilde{\xi}_\sigma(z) dz = \int_{\mathbb{R}} F_{H_s}(-z) \tilde{\xi}_\sigma(z) dz, \quad (14) $$
where $F_{H_s}$ is the unknown cdf of the univariate random variable $H_s$. Employing this new expression of the failure probability function $\overline{P_f(s, \sigma)}$ to compute the derivative with regard to $s_\ell$ with $\ell = 1, \dots, p$ gives

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \frac{\partial \overline{P_f(s, \sigma)}}{\partial s_\ell} = \int_{\mathbb{R}} \frac{\partial F_{H_s}(-z)}{\partial s_\ell} \zeta_\sigma(z) dz = \mathbb{E}_{\zeta_\sigma} \left[ \frac{\partial F_{H_s}(-Z)}{\partial s_\ell} \right], \quad (15) $$

where $Z$ is a univariate random variable of pdf $\zeta_\sigma$. Equation (15) is a new expression of the failure probability derivative function in a modified Weak approach context, in which the cdf nature of the approximation function $\Xi_\sigma$ is taken advantage of. Instead of trying to evaluate this expression of the derivative function at a specific value of $\sigma$, we use a Taylor series expansion to remove the dependence in $\sigma$, as detailed in the next section.

### 3.3. Taylor Series Expansion

The Taylor series expansion is a powerful tool to derive the expressions of expected values [13–15]. In our specific context, we derive the Taylor series expansion in the neighborhood of 0 of function $T$, defined as $T(z) = \partial F_{H_s}(-z)/\partial s_\ell$ for all $z \in \mathbb{R}$. Assuming $F_{H_s}$ is $C^{n+2}$ in 0, the Taylor series expansion leads to

$$ T(Z) = T(0) + \frac{Z}{1!} T'(0) + \frac{Z^2}{2!} T''(0) + \dots + \frac{Z^n}{n!} T^{(n)}(0) + R_n(Z), \quad (16) $$

where the remainder is expressed in the integral form with [16]

$$ R_n(Z) = \int_0^Z \frac{(Z-t)^n}{n!} T^{(n+1)}(t) dt. \quad (17) $$

Applying the expectation to both sides of the Equation (16) results in the following equation

$$ \mathbb{E}_{\zeta_{\sigma}}[T(Z)] = T(0) + \frac{\mathbb{E}_{\zeta_{\sigma}}[Z]}{1!} T'(0) + \frac{\mathbb{E}_{\zeta_{\sigma}}[Z^2]}{2!} T''(0) + \dots + \frac{\mathbb{E}_{\zeta_{\sigma}}[Z^n]}{n!} T^{(n)}(0) + \mathbb{E}_{\zeta_{\sigma}}[R_n(Z)]. \quad (18) $$

The first term $T(0)$ on the right side of the equality is equal to the derivative of the function $F_{H_s}$, with respect to $s_\ell$, evaluated in $z = 0$. This term is equal to $\partial P_f(s)/\partial s_\ell$ since $F_{H_s}(0) = P_f(s)$. Therefore, the following expression of $\partial \overline{P_f(s, \sigma)} / \partial s_\ell$ is obtained

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \frac{\partial \overline{P_f(s, \sigma)}}{\partial s_\ell} = \frac{\partial P_f(s)}{\partial s_\ell} + \frac{\mathbb{E}_{\zeta_\sigma}[Z]}{1!} T'(0) + \frac{\mathbb{E}_{\zeta_\sigma}[Z^2]}{2!} T''(0) + \dots + \frac{\mathbb{E}_{\zeta_\sigma}[Z^n]}{n!} T^{(n)}(0) + \mathbb{E}_{\zeta_\sigma}[R_n(Z)]. \quad (19) $$

From Equation (19), the derivative of $P_f$, with respect to $s_\ell$, appears to be the constant term of a polynomial expression of the moments of $Z$, where $Z$ is a random univariate variable of cdf $\Xi_\sigma$ and pdf $\zeta_\sigma$.

### 3.4. Combining Sampling Methods and Polynomial Regression to Derive the Failure Probability Sensitivity

In this section, and throughout the rest of this article, we denote $P_n$ as the polynomial of order $n$, such as

$$ \forall \sigma \in \mathbb{R}^{+} \quad P_n(\sigma) = \frac{\partial P_f(s)}{\partial s_\ell} + \frac{\mathbb{E}_{\zeta_\sigma}[Z]}{1!} T'(0) + \frac{\mathbb{E}_{\zeta_\sigma}[Z^2]}{2!} T''(0) + \dots + \frac{\mathbb{E}_{\zeta_\sigma}[Z^n]}{n!} T^{(n)}(0). \quad (20) $$
Consequently, one has

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \frac{\partial \overline{P_f(s, \sigma)}}{\partial s_{\ell}} - P_n(\sigma) = \mathbb{E}_{\xi_{\sigma}}[R_n(Z)]. $$

In order to illustrate the results derived from Equation (19), Figure 1 shows the evolution of the functions $\partial \overline{P_f(s, \cdot)} / \partial s_\ell$, $P_2$ and $P_4$, with a simple toy example in dimension one, where $g(s, X) = aX + b = H_s$, $X \sim \mathcal{N}(0, 1)$ and $s = [a, b]$. The continuous cdf $\Xi_\sigma$ selected here is the univariate centered normal cdf of variance $\sigma^2$.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustration of the polynomials $P_2$ and $P_4$ for both of the derivatives of $P_f$. The parameters of the toy example are set as follows: $a = 2$ and $b = 5$. In the figure, (a) displays the derivative with respect to $a$ while (b) displays the derivative with respect to $b$; (c,d) display the corresponding absolute difference $\frac{\partial \overline{P_f(s, \sigma)}}{\partial s_\ell} - P_n(\sigma)$. The failure probability is equal to $6.21 \times 10^{-3}$, the derivative with respect to $a$ is equal to $2.19 \times 10^{-2}$, and the derivative with respect to $b$ is equal to $8.76 \times 10^{-3}$.</figcaption></figure>

From the graphs, it is noticeable that for small values of $\sigma$, $P_n \approx \partial \overline{P_f(s, \cdot)} / \partial s_\ell$. The smaller the $\sigma$, the lower the value of $n$ needed to reach an accurate equivalence between $P_n$ and $\partial \overline{P_f(s, \cdot)} / \partial s_\ell$.

In order to evaluate the coefficients of the polynomial $P_n$, one must have at least $n + 1$ evaluations of the function $\sigma \mapsto P_n(\sigma)$. The idea of the proposed approach is, thus, to find the coefficient of $P_n$ of order zero, which is the probability sensitivity, by performing a
polynomial regression with $m$ evaluations of $(\overline{\partial P_f(s, \sigma_l) / \partial s_l})_{l=1,...,m}$ with $m \ge n + 1$ since the two functions are equivalent for small $\sigma$. Recalling that

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \frac{\overline{\partial P_f(s, \sigma)}}{\partial s_l} = \int_{\mathbb{R}^d} \frac{\partial g(s, x)}{\partial s_l} \zeta_\sigma(-g(s, x)) f_x(x) dx, \quad (21) $$

is a domain integral, it can be estimated $\forall \sigma \in \mathbb{R}^+$ with Monte Carlo methods, and in our specific case, with importance sampling, as in the Weak approach framework presented in Section 2.2. It should be noted that since $P_n \ne \partial P_f(s, \cdot) / \partial s_l$, due to the remainder of Taylor series expansion, a small global bias is introduced in the proposed approach.

The failure probability derivative function $\sigma \mapsto \partial P_f(s, \sigma) / \partial s_l$ of Equation (21) can be computed with the following IS estimate

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \frac{\overline{\partial P_f(s, \sigma)}^{\text{IS}}}{\partial s_l} = -\frac{1}{N} \sum_{j=1}^{N} \frac{\partial g(s, \mathbf{X}^{(j)})}{\partial s_l} \zeta_{\sigma}(-g(s, \mathbf{X}^{(j)})) \frac{f_{\mathbf{X}}(\mathbf{X}^{(j)})}{h(\mathbf{X}^{(j)})}, \quad (22) $$

where the observations $\mathbf{X}^{(j)}$ are iid from $h$, the IS auxiliary density. These observations are reused from the estimation of the failure probability with an IS estimate. Therefore, the only additional simulation budget concerns the evaluation of the limit state function derivative with respect to $s_l$ for each observation. The IS estimate is unbiased and its variance is given by

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \mathrm{Var} \left( \frac{\overline{\partial P_f(s, \sigma)}^{\mathrm{IS}}}{\partial s_l} \right) = \frac{1}{N} \mathrm{Var}_h \left( \frac{\partial g(s, \mathbf{X})}{\partial s_l} \zeta_{\sigma}(-g(s, \mathbf{X})) \frac{f_{\mathbf{X}}(\mathbf{X})}{h(\mathbf{X})} \right). \quad (23) $$

Depending on the pdf $\zeta_\sigma$, the expression of this variance can be further detailed. For instance, if $\zeta_\sigma$ is the pdf of a centered normal variable of variance $\sigma^2$, then $\forall y \in \mathbb{R} \ \zeta_\sigma(y) = (1/\sigma)\phi(y/\sigma)$, and the variance becomes

$$ \forall \sigma \in \mathbb{R}^{+} \setminus \{0\} \quad \mathrm{Var} \left( \frac{\overline{\partial P_f(s, \sigma)}^{\mathrm{IS}}}{\partial s_l} \right) = \frac{1}{N\sigma^2} \mathrm{Var}_h \left( \frac{\partial g(s, \mathbf{X})}{\partial s_l} \phi \left( \frac{-g(s, \mathbf{X})}{\sigma} \right) \frac{f_{\mathbf{X}}(\mathbf{X})}{h(\mathbf{X})} \right). \quad (24) $$

Therefore, for a fixed $N$, each estimate $\overline{\partial P_f(s, \sigma_l) / \partial s_l}$, for $l = 1, ..., m$, has a different noise. In order to identify the coefficients of the polynomial $P_n$, a heteroscedastic polynomial regression must then be performed. This heteroscedastic polynomial regression is further detailed in the next section.

## 4. Heteroscedastic Polynomial Regression

The main parameters of the heteroscedastic polynomial regression are the selected degree $n$ of the polynomial, the number $m$ of estimates $(\overline{\partial P_f(s, \sigma_l) / \partial s_l})_{l=1,...,m}$ with $m \ge n + 1$, and the interval of regression $[\sigma_{\min}, \sigma_{\max}]$, given that $n$ and the interval are dependent. However, the correlation between the different estimates $(\overline{\partial P_f(s, \sigma_l) / \partial s_l})_{l=1,...,m}$ is first discussed, as it directly influences the regression framework.

We denote $\mathbf{V}$ as the $m$-length vector of the estimates, given by $\mathbf{V} = (\overline{\partial P_f(s, \sigma_l) / \partial s_l})_{l=1,...,m}$, $\mathbf{S}$ as the Vandermonde matrix of $\sigma$ of size $m \times (n + 1)$, where $\mathbf{S}_{l,i} = \sigma_l^i$, for $l = 1, ..., m$ and
$i = 0, \dots, n$, and $\alpha$ denotes the polynomial coefficients, such that the regression is written as follows

$$V = S\alpha + \epsilon, \qquad (25)$$

where $\epsilon$ is the $m$-length vector of random errors $\epsilon_l$ with expected value $\mathbb{E}_{f_X}[\epsilon_l] = 0$, as the estimates $V_l$ are unbiased. The error's variance $\mathbb{V}ar(\epsilon_l)$ depends on the $\sigma_l$ cf Equation (23). The $\sigma$-vector is denoted in ascending order: $\sigma_1 > \sigma_2 > \dots > \sigma_m$. This linear regression is addressed in this article with linear least square methods [17,18].

### 4.1. Linear Least Squares Method in Our Specific Context

As previously mentioned, the observations $(\mathbf{X}^{(j)})_{j=1,...,N}$ needed for each estimate $V_l$ are reused from the failure probability estimation procedure to minimize the additional simulation budget. However, reusing the exact same sample for each $V_l$ results in a highly correlated database $(\sigma_l, V_l)_{l=1,...,m}$. The generalized least squares (GLS) framework extends the generalization of the ordinary least squares (OLS) framework to situations where the variance in the random error, denoted as $\epsilon_l$, is not constant (heteroscedasticity) and the estimate vector, $V$, is not uncorrelated. We denote $\Sigma_\epsilon$ as the covariance matrix of the random errors $\epsilon$ and $\hat{\Sigma}_\epsilon$ as its corresponding empirical estimate. The GLS estimate of the coefficient $\alpha$ is then obtained with the following equations

$$\hat{\alpha} = (\mathbf{S}^{\top}\Sigma_{\epsilon}^{-1}\mathbf{S})^{-1}(\mathbf{S}^{\top}\Sigma_{\epsilon}^{-1}\mathbf{V}), \qquad (26)$$

and it is the best linear unbiased estimate [19]. It is possible to derive the covariance matrix estimate of $\hat{\alpha}$ with the following formula:

$$\widehat{\mathbb{V}ar}(\hat{\alpha}) = (\mathbf{S}^{\top}\Sigma_{\epsilon}^{-1}\mathbf{S})^{-1}. \qquad (27)$$

In our specific context, the covariance matrix $\Sigma_\epsilon$ is not analytically known and has to be estimated with the sample $(\mathbf{X}^{(j)})_{j=1,...,N}$. Therefore, in Equations (26) and (27), replacing $\Sigma_\epsilon$ with $\hat{\Sigma}_\epsilon$ gives the final framework suited for the heteroscedastic polynomial regression needed in our approach, called the feasible GLS (FGLS).

From Equations (26) and (27), one can notice that the inverse of the covariance matrix is required to compute both quantities. In our specific case, this covariance matrix is an estimation of the real covariance matrix. Therefore, each coefficient already comes with an estimation error. All those errors combined make it very difficult to accurately compute $\hat{\Sigma}_\epsilon^{-1}$ as the conditioning number of the matrix $\hat{\Sigma}_\epsilon$ is generally very high, making its numerical inversion challenging. Consequently, although the FGLS framework is theoretically best suited for the proposed method, it is practically inapplicable. For this reason, a regression framework simplification must be considered, which is detailed next.

The simplification proposed here is to decrease the correlation of the vector of estimates $V$ with bootstrap [20]. The logic of bootstrapping is to learn an empirical discrete cdf from a vector of iid observations in order to generate a new sample. The new sample shares the same property as the original sample but they are independent. Here, we apply bootstrap to the iid observations $(\mathbf{X}^{(j)})_{j=1,...,N}$ already available from the failure probability estimation. The first estimate $V_1$ is computed with the original sample, and for each $\sigma_l$ with $l > 1$, bootstrap is employed. The resulting vector $V$ is, thus, independent. Bootstrap does not require any additional call to the limit state function or its derivatives.

The resulting covariance matrix $\hat{\Sigma}_\epsilon$ is, thus, assumed diagonal. The regression framework is then referred to as weighted least squares (WLS). Denoting $\mathbf{W}$ as the diagonal matrix, where $W_{l,l} = 1/\widehat{\mathbb{V}ar}(V_l)$, and where $\widehat{\mathbb{V}ar}(V_l)$ is given by Equation (23), the WLS estimate is then written

$$\hat{\alpha}_W = (\mathbf{S}^{\top}\mathbf{W}\mathbf{S})^{-1}(\mathbf{S}^{\top}\mathbf{W}\mathbf{V}), \qquad (28)$$
and it is the best linear unbiased estimate [19]. It is possible to derive the covariance matrix estimate of $\hat{\alpha}_W$ with the following formula

$$ \widehat{\mathrm{Var}}(\hat{\alpha}_W) = (\mathbf{S}^{\top} \mathbf{W} \mathbf{S})^{-1}. \qquad (29) $$

Therefore, with only one simulation run, an estimation of the variance of the sensitivity estimate is also available. The sensitivity estimate is equal to the coefficient of order zero of the polynomial, which is the first component of the vector $\hat{\alpha}_W$, and its theoretical variance estimate is the component $\widehat{\mathrm{Var}}(\hat{\alpha}_W)_{1,1}$.

## 4.2. Settings of the Regression Parameters

The settings presented here are the results of several tests. Before presenting the settings, the scaling of the limit state function is first addressed, as it greatly influences the evolution of the failure probability derivative functions $\sigma \mapsto \partial P_f(\mathbf{s}, \sigma) / \partial s_\ell$.

### 4.2.1. Scaling of the Limit State Function

Depending on the limit state function, the order of magnitude of $\mathbf{x} \mapsto g(\mathbf{s}, \mathbf{x})$ and its derivatives, with respect to $\mathbf{s}$, can significantly vary, especially in the vicinity of the failure surface, as underlined in [4,12]. As a result, the behavior of the failure probability derivative functions $\sigma \mapsto \partial P_f(\mathbf{s}, \sigma) / \partial s_\ell$ can considerably vary as well. Here, general regression settings are presented, which aim to be applied to various limit state functions. Consequently, the following scaling is first performed.

Let $(Y^{(j)} = g(\mathbf{s}, \mathbf{X}^{(j)}) \quad j=1,...,N)$ be the lsf values vector obtained from the estimation of the failure probability. The scaling proposed here is to divide the vector $\mathbf{Y}$ as well as the vector of derivatives $(\partial g(\mathbf{s}, \mathbf{X}^{(j)}) / \partial s_\ell)_{j=1,...,N}$ by the standard deviation $\delta$ obtained from the negative lsf values only of $\mathbf{Y}$

$$ \delta^2 = \sum_{j=1}^{N} \bar{w}^{(j)} (Y^{(j)} - \mu)^2 \quad \text{with} \quad \mu = \sum_{j=1}^{N} \bar{w}^{(j)} Y^{(j)}, \qquad (30) $$

where $\bar{w}^{(j)} = w^{(j)} / \sum_{i=1}^{N} w^{(i)}$ and $w^{(j)} = \mathbb{I}_{Y<0} (Y^{(j)})$. Therefore, $\delta$ represents the order of magnitude of the lsf in the failure domain.

This division neither affects the value of the failure probability nor the failure probability sensitivity. Indeed, provided that $\delta > 0$, we have $P_f(\mathbf{s}) = \mathbb{P}(g(\mathbf{s}, \mathbf{X}) < 0) = \mathbb{P}(g(\mathbf{s}, \mathbf{X}) / \delta < 0)$, and the derivatives do not change. However, this division influences the value of the coefficients of the order superior to zero in the Taylor series expansion of $T$. Since these coefficients are not of interest in the proposed approach, the consequences are negligible. Thanks to this scaling, the order of magnitude of the lsf is expected to have little influence on the polynomial regression settings presented in the following paragraphs.

### 4.2.2. Choice of the Regression Interval and the Polynomial Degree

**Definition of the regression interval** $[\sigma_{\min}, \sigma_{\max}]$.

The selection of the interval is crucial in the regression process, as it greatly influences the quality of the sensitivity estimate. If the lower bound of the interval is too close to zero, then the estimates $(V_l)_{l=1,...,m}$ in the vicinity of the lower bound will be very noisy and their variance might be inaccurately estimated, as $N$ is fixed. Consequently, the WLS framework might lead to erroneous results. Moreover, if the upper bound of the regression interval is too far from zero, a higher degree polynomial is needed to correctly approximate the Taylor series expansion. The polynomial regression is then harder to achieve and leads to a probability sensitivity of higher variance. Therefore, there is a trade-off between the accuracy of the estimates $(V_l)_{l=1,...,m}$ and $(\widehat{\mathrm{Var}}(V_l))_{l=1,...,m}$, and the variance of the
sensitivity estimate obtained with the polynomial regression. In the proposed approach, it was decided to define the interval bounds inside the interval $\sigma \in [0.01, 1]$. Thanks to the scaling process mentioned above, restricting the search for the regression bounds in this interval has proven to be efficient for the various limit state functions tested.

Within this interval, the regression bounds are set, depending on the theoretical coefficient of variation (CV) estimate of the IS estimate $V_1$. The CV of $V_1$ is equal to the ratio between the square root of its variance, given by Equation (23), and its value, given by Equation (22). This CV is, thus, a function of $\sigma$. Choosing a CV criterion to set the interval bounds is inspired by [4], where the optimal $\tilde{\sigma}$ of the Weak approach framework is selected with a target CV technique. Here, the values of $\sigma_{min}$ and $\sigma_{max}$ are, thus, respectively, the lowest and the highest values of $\sigma$, such as $CV(\sigma) < CV_{target}$ for all $\sigma \in [\sigma_{min}, \sigma_{max}]$. This CV criterion ensures the regression is performed with estimates $(V_l)_{l=1,...,m}$ and $(\tilde{\text{Var}}(V_l))_{l=1,...,m}$, which have reasonable noise. In this article, the $CV_{target}$ is dependently set on the minimum theoretical CV estimate of the $V_1$ estimate with

$$CV_{target} = \min_{\sigma \in [0.01,1]} CV_{V_1}(\sigma) + 5\%. \qquad (31)$$

The additional 5% ensures a sufficiently large interval and is an arbitrary value. It has proven to be efficient after several tests. This threshold is independent of the quality of the estimation of $P_f$, as it has been underlined that $P_f$ and its derivatives do not depend on the same quantities [8]. Furthermore, this threshold changes for each deterministic input $s_\ell$ and allows the regression interval to be specifically suited for the evolution of $\sigma \mapsto \partial P_f(s, \sigma)/\partial s_\ell$.

Once the bounds are computed, the values of the vector $(\sigma_l)_{l=1,...,m}$ are set uniformly in the interval $[\sigma_{min}, \sigma_{max}]$.

**Selection of the polynomial degree $n$ and number $m$ of estimations.** The choice of the polynomial degree $n$ influences the global bias of the proposed approach. The higher the degree $n$, the lower the bias for a fixed regression interval, as illustrated in Figure 1. However, a high degree $n$ induces a more intricate polynomial regression, as the polynomial $P_n$ is then more complex, resulting in a sensitivity estimate with higher variance. Therefore, there is a trade-off between the value of the theoretical bias of the sensitivity estimate and its variance.

The choice of $n$ affects the number $m$ of estimations $V_l$ required for the regression, as there are $n+1$ coefficients that have to be estimated with the regression. For simplification purposes, the number $m$ of estimations of the failure probability derivatives $V_l = (\partial \widehat{P_f(s, \sigma_l)}/\partial s_\ell)_{l=1,...,m}$ is set to $m = n+2$, as several tests showed that an increase in this number did not improve the quality of the sensitivity estimates.

These $m$ estimations are obtained with bootstrap; therefore, they do not increase the simulation budget. Consequently, several tests can be performed to assess the minimal degree $n$ needed for a correct sensitivity estimate, without affecting the simulation budget. If—with a higher degree polynomial—the resulting sensitivity estimate value no longer changes, we can assume that the estimate has reached the correct value and the bias is controlled since it has become negligible. This control of the bias depends naturally on the level of accuracy wanted for the estimation of the sensitivity.

In the numerical evaluations of the following section, different values of $n$ are, thus, set to identify the minimal degree needed to obtain a sensitivity estimate with controlled bias. However, different model selection methods, like AIC (Akaike information criterion) and BIC (Bayesian information criterion) [21], can be used to choose the best statistical model from a set of candidate models. These criteria help balance the trade-off between model complexity and goodness of fit and can, thus, be an improvement to the proposed algorithm.
## 5. Numerical Investigation

The performance of the proposed approach is investigated with two numerical applications, taken from the failure probability sensitivity literature. These examples focus on applications in a rather low-dimensional space and feature a singular failure: the cantilever beam [4,5,8,11] and the roof truss [4,5,22,23]. As they are very common examples in the sensitivity literature, they are relevant cases to present the proposed approach. Function g is known analytically in both use cases but this is not a requirement to apply the proposed algorithm as the method only calls on g and $\partial g / \partial s_i$ in a point-wise manner.

In the first example, the deterministic inputs s are design parameters, while in the second example, they are the distribution parameters of the original inputs, denoted as Z. An isoprobabilistic transformation allows transforming the original inputs Z into the standard normal inputs X. The limit state function in the original space is denoted as gZ while the transformed limit state function is denoted as g.

The proposed approach is combined with various IS algorithms for the numerical investigation. The non-parametric adaptive IS (NAIS) algorithm [24] for the cantilever beam and the iCE-AIS algorithm [25] for the roof truss are considered in the following to show the robustness of the proposed approach to sampling algorithms. For both algorithms, the last generated sample of the adaptive procedure is used to compute the sensitivity with the new approach.

For all applications, the indicator approximation function $\Xi_{\sigma}$ selected is the cdf of a centered normal random variable of variance $\sigma^2$. As previously mentioned, other cdfs of continuous univariate variables could be relevant [8]. However, the Gaussian approximation, being the most considered in practice [4,11,12], is the one selected here. The impact of the selection of another approximation $\Xi_{\sigma}$ is not studied in this paper and is left to future work. This particular choice results in even polynomials $P_n$; therefore, $n = 2k$ and the number of $V_i$ estimates is fixed to $m = k + 2$. The proposed approach is studied for three different polynomial degrees $n = 2, n = 4$, and $n = 6$. It is compared to the Weak approach, when $\tilde{\sigma} = \sigma_{min}$, the lower bound of the regression interval. For comparison purposes, 500 independent simulation runs are performed to calculate the statistics of the probability estimates and the other quantities of interest.

### 5.1. Cantilever Beam

#### 5.1.1. Presentation of the Application

The first example is a cantilever beam subject to biaxial bending, as illustrated in Figure 2. This example is quite popular in sensitivity analysis [4,5,8,11] and was first studied in [26] in an RBDO context.

<figure><img src="image_5.png" /><figcaption>Figure 2. Illustration of a cantilever beam subject to biaxial bending.</figcaption></figure>

Two different limit state functions are considered for this example, both defined in the original space Z of dimension 4. The first one $g_Z^{(1)}$ represents yielding at the fixed end of the beam, with s = [w, t], respectively, denoting the width and the height of the cross-section beam

$$g_Z^{(1)}(\mathbf{s}, \mathbf{Z}) = Z_3 - \left( \frac{600}{wt^2} Z_1 + \frac{600}{w^2t} Z_2 \right). \quad (32)$$
The second limit state function $g_Z^{(2)}$ restricts the maximum allowed displacement at the tip of the beam to the value $d_0$. Therefore, one has $\mathbf{s} = [w, t, d_0]$

$$g_Z^{(2)}(\mathbf{s}, \mathbf{Z}) = d_0 - \frac{4L^3}{Z_4 w t} Z_3 \sqrt{\left(\frac{Z_1}{t^2}\right)^2 + \left(\frac{Z_2}{w^2}\right)^2} \quad (33)$$

where $L = 100$ m. The random variables $Z_1$ and $Z_2$ represent the loads, $Z_3$ is the yield strength of the beam, and $Z_4$ is the Young's modulus. We assume the vector $\mathbf{Z}$ to be an independent normal vector. For comparison purposes, the distribution of each random variable is the same as in [26]; therefore, denoting $\mu_{Z_i}$ as the mean value and $\delta_{Z_i}$ as the standard deviation, one has

$$(\mu_{Z_1}, \delta_{Z_1}) = (1000, 100) \qquad (34)$$

$$(\mu_{Z_2}, \delta_{Z_2}) = (500, 100) \qquad (35)$$

$$(\mu_{Z_3}, \delta_{Z_3}) = (40000, 2000) \qquad (36)$$

$$(\mu_{Z_4}, \delta_{Z_4}) = (29 \times 10^6, 1.45 \times 10^6). \qquad (37)$$

The isoprobabilistic transformation is quite simple in this case, as the inputs are already independent. Therefore, $X_i = (Z_i - \mu_{Z_i})/\delta_{Z_i}$ for $i = 1, \dots, 4$; this transformation is linear.

The design parameters are fixed as $[w, t, d_0] = [2.4, 3.9, 2.5]$, which corresponds to the optimal reliability-based design [4].

### 5.1.2. Sensitivity Analysis for the First Failure of the System

As the first limit state function $g_Z^{(1)}$ is linear and the transformation from $\mathbf{Z}$ to $\mathbf{X}$ is linear as well, the transformed lsf $g^{(1)}$ is also linear in the random variables $\mathbf{X}$. Consequently, the failure probability as well as its sensitivities can be exactly determined with a FORM analysis [27]. This FORM analysis gives the following reference values [4]: $P_f^{(1)} = 3.03 \times 10^{-3}$, $\partial P_f^{(1)} / \partial w = -5.76 \times 10^{-2}$, and $\partial P_f^{(1)} / \partial t = -3.53 \times 10^{-2}$.

For this first application, the NAIS algorithm is combined with the proposed method. The mean simulation budget required for the probability estimation is near 8000. The results are presented in Table 1, with the empirical CVs given, as well as the theoretical CV estimates of the proposed method in parentheses.

<table><caption>**Table 1.** Comparison of the results of the polynomial regression with 3 different degrees, for the first lsf of the cantilever beam with NAIS. The failure probability is equal to $2.97 \times 10^{-3}$ with an empirical CV of 8.4%. The reference values for the sensitivity are $\partial P_f^{(1)} / \partial w = -5.76 \times 10^{-2}$ and $\partial P_f^{(1)} / \partial t = -3.53 \times 10^{-2}$. Empirical CVs of the proposed method are given in %, and the theoretical CV estimates are in parentheses.</caption><thead><tr><th colspan="8">Cantilever Beam 1</th></tr><tr><th></th><th>Regression<br>2k = 2</th><th>Regression<br>2k = 4</th><th>Regression<br>2k = 6</th><th colspan="3">Weak Approach<br>σ̃ = σ<sub>min</sub></th></tr></thead><tbody><tr><td>∂P<sub>f</sub>/∂w</td><td>-5.76 × 10<sup>-2</sup><br>CV</td><td>-5.74 × 10<sup>-2</sup><br>CV</td><td>-5.76 × 10<sup>-2</sup><br>CV</td><td>-5.76 × 10<sup>-2</sup><br>CV</td><td>-5.76 × 10<sup>-2</sup><br>CV</td><td>7.0%</td></tr><tr><td>∂P<sub>f</sub>/∂t</td><td>-3.53 × 10<sup>-2</sup><br>CV</td><td>-3.52 × 10<sup>-2</sup><br>CV</td><td>-3.53 × 10<sup>-2</sup><br>CV</td><td>-3.53 × 10<sup>-2</sup><br>CV</td><td>-3.53 × 10<sup>-3</sup><br>CV</td><td>7.0%</td></tr></tbody></table>

From Table 1, the sensitivity estimates of both the proposed approach and the Weak approach are very close to the reference values of the failure probability sensitivities.
For the three degrees selected and the two derivatives, the value of the estimate is quite constant. Therefore, a polynomial of degree 2 is sufficient for correctly approximating the coefficient of order zero of the Taylor series expansion of the function $\sigma \mapsto \partial \overline{P_f(s, \sigma)} / \partial s_l$. The sensitivity estimate obtained with $2k = 2$ has a smaller CV than the one obtained with the Weak approach. Consequently, this application illustrates how the proposed method is an improvement of the Weak approach; it results in a more precise estimate, without any additional simulation budget. However, the theoretical CV estimates are slightly underrated. They have large variations (above 45%, not displayed in Table 1); therefore, they cannot be considered accurately estimated.

### 5.1.3. Sensitivity Analysis for the Second Failure of the System

The limit state function $g_Z^{(2)}$ is nonlinear in all the random variables; therefore, FORM cannot be applied to accurately compute the failure probability and its sensitivities. The reference values used here are, thus, the ones given in [4], which are the outcomes of a large line sampling simulation. The line sampling simulation gives the following values: $P_f^{(2)} = 2.54 \times 10^{-4}$, $\partial P_f^{(2)} / \partial w = -8.84 \times 10^{-3}$, $\partial P_f^{(2)} / \partial t = -2.95 \times 10^{-3}$ and $\partial P_f^{(2)} / \partial d_0 = -3.27 \times 10^{-3}$.

For this second application, the proposed approach is combined with the iCE-AIS presented in [25], with a single Gaussian density as the auxiliary density. The mean simulation budget required for the probability estimation is near 4000. The results are presented in Table 2, with the empirical CVs given, as well as the theoretical CV estimates of the proposed method in parentheses.

<table><caption>**Table 2.** Comparison of the results of the polynomial regression with 3 different degrees, for the second lsf of the cantilever beam with iCE-SG. The failure probability is equal to $2.51 \times 10^{-3}$ with an empirical CV of 4.5%. The reference values for the sensitivity are $\partial P_f^{(2)} / \partial w = -8.84 \times 10^{-4}$, $\partial P_f^{(2)} / \partial t = -2.95 \times 10^{-3}$ and $\partial P_f^{(2)} / \partial d_0 = -3.27 \times 10^{-3}$. Empirical CVs of the proposed method are given in %, as well as the theoretical CV estimates in parentheses.</caption><thead><tr><th colspan="7">Cantilever Beam 2</th></tr><tr><th></th><th>Regression<br>2k = 2</th><th>Regression<br>2k = 4</th><th>Regression<br>2k = 6</th><th colspan="2">Weak Approach</th></tr><tr><th></th><th></th><th></th><th></th><th colspan="2">&sigma;&#770; = &sigma;<sub>min</sub></th></tr></thead><tbody><tr><td>&Delta;P&#770;<sub>f</sub>/&Delta;w</td><td>-8.78 &times; 10<sup>-3</sup><br>CV</td><td>-8.82 &times; 10<sup>-3</sup><br>CV</td><td>-8.83 &times; 10<sup>-3</sup><br>CV</td><td colspan="2">-8.93 &times; 10<sup>-3</sup><br>CV</td></tr><tr><td></td><td></td><td>7.3%<br>(5.0%)</td><td>8.6%<br>(6.8%)</td><td colspan="2">9.9%<br>(8.8%)</td></tr><tr><td>&Delta;P&#770;<sub>f</sub>/&Delta;t</td><td>-2.92 &times; 10<sup>-3</sup><br>CV</td><td>-2.94 &times; 10<sup>-3</sup><br>CV</td><td>-2.94 &times; 10<sup>-3</sup><br>CV</td><td colspan="2">-2.98 &times; 10<sup>-3</sup><br>CV</td></tr><tr><td></td><td></td><td>7.4%<br>(5.0%)</td><td>8.6%<br>(6.8%)</td><td colspan="2">9.9%<br>(8.8%)</td></tr><tr><td>&Delta;P&#770;<sub>f</sub>/&Delta;d<sub>0</sub></td><td>-3.24 &times; 10<sup>-3</sup><br>CV</td><td>-3.26 &times; 10<sup>-3</sup><br>CV</td><td>-3.26 &times; 10<sup>-3</sup><br>CV</td><td colspan="2">-3.29 &times; 10<sup>-3</sup><br>CV</td></tr><tr><td></td><td></td><td>7.3%<br>(5.0%)</td><td>8.6%<br>(6.8%)</td><td colspan="2">9.9%<br>(8.8%)</td></tr></tbody></table>

From Table 2, it is noticeable that the polynomial of degree $2k = 2$ leads to sensitivities whose values are slightly different from those obtained with higher degrees, for the three parameters. Therefore, it appears that a polynomial of degree 2 is not sufficient for accurately estimating the coefficient of order zero of the Taylor series expansion of the functions $\sigma \mapsto \partial \overline{P_f(s, \sigma)} / \partial s_l$.

The estimate obtained with $2k = 2$ shares the same properties as the Weak approach estimate. Their biases are quite similar along with their CVs. For $2k \ge 4$, the estimates obtained with the proposed approach have a smaller bias than the estimates obtained with the Weak approach, but they have a slightly higher CV. The theoretical CV estimates are still slightly underrated. They have a moderate variation for $2k \le 4$, with an empirical
CV that is close to 25% (not displayed in Table 2). For $2k = 6$, the variation is quite high, i.e., 39%.

Consequently, for this application, the proposed approach brings improvement to the Weak approach framework, as estimates with smaller biases can be obtained and the biases are globally controlled. For $2k \ge 4$—since the values of the estimates no longer vary—it can be assumed that the estimates have reached accurate values. Such analysis cannot be performed with the Weak approach, where the bias is not properly managed.

## 5.2. Roof Truss

### 5.2.1. Presentation of the Application

The second example is a roof truss subject to random loading as illustrated in Figure 3. This example is also very commonly used in sensitivity analysis [4,5,22,23], and we keep the same framework as presented in [23]. The top boom and the compression bars are reinforced by concrete; the bottom boom and the tension bars are made of steel.

<figure><img src="image_6.png" /><figcaption>Figure 3. Illustration of a roof truss subject to random loading.</figcaption></figure>

The perpendicular deflection of the peak of the structure must not exceed 3 cm [23]. Consequently, the limit state function is defined in the original space Z of dimension 6 with the following equation

$$g_{\mathbf{Z}}^{(3)}(\mathbf{s}, \mathbf{Z}) = 0.03 - \frac{Z_1 Z_2^2}{2} \left( \frac{3.81}{Z_4 Z_6} + \frac{1.13}{Z_3 Z_5} \right) \qquad (38)$$

where $\mathbf{s} = [\mu_{Z_1}, \delta_{Z_1}, \dots, \mu_{Z_6}, \delta_{Z_6}]$, with $\delta_{Z_i}$ denoting the standard deviation of $Z_i$. The random variable $Z_1$ represents the uniformly distributed load applied on the roof truss and $Z_2$ is the roof span. $Z_3$ is the cross-section area of the bottom boom and the tension bars made of steel, whose Young's modulus is $Z_5$. $Z_4$ is the cross-section area of the top boom and compression bars reinforced in concrete, whose Young's modulus is $Z_6$. We assume that vector Z is an independent normal vector. For comparison purposes, the distribution of each random variable is the same as in [23]. Therefore, we have

$$(\mu_{Z_1}, \delta_{Z_1}) = (20000, 1400) \qquad (39)$$

$$(\mu_{Z_2}, \delta_{Z_2}) = (12, 0.12) \qquad (40)$$

$$(\mu_{Z_3}, \delta_{Z_3}) = (9.82 \times 10^{-4}, 5.892 \times 10^{-5}) \qquad (41)$$

$$(\mu_{Z_4}, \delta_{Z_4}) = (0.04, 0.0048) \qquad (42)$$

$$(\mu_{Z_5}, \delta_{Z_5}) = (1 \times 10^{11}, 6 \times 10^9) \qquad (43)$$

$$(\mu_{Z_6}, \delta_{Z_6}) = (2 \times 10^{10}, 1.2 \times 10^9) \qquad (44)$$

As for the previous application, it should be noted that the normal distribution is strictly not an appropriate choice for modeling $Z_2, Z_3, Z_4, Z_5$, and $Z_6$, as they represent
physical variables of positive support. Once more, the isoprobabilistic transformation is simple in this case, as the inputs are already independent, resulting in $X_i = (Z_i - \mu_{Z_i})/\delta_{Z_i}$ for $i = 1, \dots, 6$. This transformation is linear.

### 5.2.2. Sensitivity Analysis of the System

The reference values of the failure probability and its sensitivities are the same as in [4], obtained with the score function method, and combined with IS [28]. The resulting failure probability is equal to $9.38 \times 10^{-3}$ and the sensitivities are equal to

$$ (\partial P_f / \partial \mu_{Z_1}, \partial P_f / \partial \delta_{Z_1}) = (1.11 \times 10^{-5}, 1.59 \times 10^{-5}) \quad (45) $$

$$ (\partial P_f / \partial \mu_{Z_2}, \partial P_f / \partial \delta_{Z_2}) = (4.03 \times 10^{-2}, 1.80 \times 10^{-2}) \quad (46) $$

$$ (\partial P_f / \partial \mu_{Z_3}, \partial P_f / \partial \delta_{Z_3}) = (-1.86 \times 10^2, 2.05 \times 10^2) \quad (47) $$

$$ (\partial P_f / \partial \mu_{Z_4}, \partial P_f / \partial \delta_{Z_4}) = (-2.14, 2.56) \quad (48) $$

$$ (\partial P_f / \partial \mu_{Z_5}, \partial P_f / \partial \delta_{Z_5}) = (-1.83 \times 10^{-12}, 2.00 \times 10^{-12}) \quad (49) $$

$$ (\partial P_f / \partial \mu_{Z_6}, \partial P_f / \partial \delta_{Z_6}) = (-3.77 \times 10^{-12}, 2.03 \times 10^{-12}) \quad (50) $$

For this application, the NAIS algorithm is combined with the proposed method. The mean simulation budget required for the probability estimation is near 6000. The results are presented in Tables 3 and 4, with the empirical CVs given; the theoretical CV estimates of the proposed method are in parentheses.

<table><caption>**Table 3.** Comparison of the results of the polynomial regression with 3 different degrees for the roof truss with NAIS. The failure probability is equal to $9.28 \times 10^{-3}$ with an empirical CV of 12.5%. The reference values for the sensitivity are $(\partial P_f / \partial \mu_{Z_1}, \partial P_f / \partial \delta_{Z_1}) = (1.11 \times 10^{-5}, 1.59 \times 10^{-5})$, $(\partial P_f / \partial \mu_{Z_2}, \partial P_f / \partial \delta_{Z_2}) = (4.03 \times 10^{-2}, 1.80 \times 10^{-2})$ and $(\partial P_f / \partial \mu_{Z_3}, \partial P_f / \partial \delta_{Z_3}) = (-1.86 \times 10^2, 2.05 \times 10^2)$. Empirical CVs of the proposed method are given in %, as well as the theoretical CV estimates in parentheses.</caption><thead><tr><th colspan="8">Roof Truss, Part 1</th></tr><tr><th></th><th>Regression 2k = 2</th><th>Regression 2k = 4</th><th>Regression 2k = 6</th><th>Weak Approach</th><th></th><th></th><th></th></tr></thead><tbody><tr><td>$$ \frac{\partial P_f}{\partial \mu_{Z_1}} $$</td><td>1.10 x 10<sup>-5</sup></td><td>1.10 x 10<sup>-5</sup></td><td>1.10 x 10<sup>-5</sup></td><td>1.11 x 10<sup>-5</sup></td><td></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>7.4%<br>(5.2%)</td><td>CV</td><td>9.0%<br>(7.2%)</td><td>CV</td><td>14.4%<br>(9.2%)</td><td>CV</td><td>8.2%</td></tr><tr><td>$$ \frac{\partial P_f}{\partial \delta_{Z_1}} $$</td><td>1.57 x 10<sup>-5</sup></td><td>1.57 x 10<sup>-5</sup></td><td>1.58 x 10<sup>-5</sup></td><td>1.58 x 10<sup>-5</sup></td><td></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>8.2%<br>(5.8%)</td><td>CV</td><td>10.4%<br>(8.2%)</td><td>CV</td><td>15.6%<br>(11.2%)</td><td>CV</td><td>9.2%</td></tr><tr><td>$$ \frac{\partial P_f}{\partial \mu_{Z_2}} $$</td><td>4.01 x 10<sup>-2</sup></td><td>4.03 x 10<sup>-2</sup></td><td>4.00 x 10<sup>-2</sup></td><td>4.05 x 10<sup>-2</sup></td><td></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>7.4%<br>(5.2%)</td><td>CV</td><td>9.0%<br>(7.2%)</td><td>CV</td><td>16.0%<br>(9.2%)</td><td>CV</td><td>8.3%</td></tr><tr><td>$$ \frac{\partial P_f}{\partial \delta_{Z_2}} $$</td><td>1.79 x 10<sup>-2</sup></td><td>1.88 x 10<sup>-2</sup></td><td>2.04 x 10<sup>-2</sup></td><td>1.85 x 10<sup>-2</sup></td><td></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>18.7%<br>(15.1%)</td><td>CV</td><td>36.1%<br>(26.0%)</td><td>CV</td><td>375%<br>(12.8%)</td><td>CV</td><td>14.4%</td></tr><tr><td>$$ \frac{\partial P_f}{\partial \mu_{Z_3}} $$</td><td>-1.85 x 10<sup>2</sup></td><td>-1.85 x 10<sup>2</sup></td><td>-1.84 x 10<sup>2</sup></td><td>-1.86 x 10<sup>2</sup></td><td></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>7.4%<br>(5.2%)</td><td>CV</td><td>9.1%<br>(7.2%)</td><td>CV</td><td>14.8%<br>(9.2%)</td><td>CV</td><td>8.3%</td></tr><tr><td>$$ \frac{\partial P_f}{\partial \delta_{Z_3}} $$</td><td>2.02 x 10<sup>2</sup></td><td>2.03 x 10<sup>2</sup></td><td>2.05 x 10<sup>2</sup></td><td>2.03 x 10<sup>2</sup></td><td></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>10.0%<br>(7.3%)</td><td>CV</td><td>13.6%<br>(9.1%)</td><td>CV</td><td>26.1%<br>(14.0%)</td><td>CV</td><td>10.1%</td></tr></tbody></table>
<table><caption>**Table 4.** Comparison of the results of the polynomial regression with 3 different degrees for the roof truss with NAIS. The failure probability is equal to $9.28 \times 10^{-3}$ with an empirical CV of 12.5%. The reference values for the sensitivity are $(\partial P_f / \partial \mu_{Z_4}, \partial P_f / \delta_{Z_4}) = (-2.14, 2.56)$, $(\partial P_f / \partial \mu_{Z_5}, \partial P_f / \delta_{Z_5}) = (-1.83 \times 10^{-12}, 2.00 \times 10^{-12})$ and $(\partial P_f / \partial \mu_{Z_6}, \partial P_f / \delta_{Z_6}) = (-3.77 \times 10^{-12}, 2.03 \times 10^{-12})$. Empirical CVs of the proposed method are given in %, as well as the theoretical CV estimates in parentheses.</caption><thead><tr><th colspan="7">Roof Truss, Part 2</th></tr><tr><th></th><th>Regression 2k = 2</th><th>Regression 2k = 4</th><th>Regression 2k = 6</th><th>Weak approach</th><th></th><th></th></tr><tr><th></th><th></th><th></th><th></th><th>σ̄ = σ<sub>min</sub></th><th></th><th></th></tr></thead><tbody><tr><td>∂P<sub>f</sub>/∂μ<sub>Z<sub>4</sub></sub></td><td>-2.11</td><td>-2.11</td><td>-2.11</td><td>-2.13</td><td></td><td></td></tr><tr><td></td><td>CV</td><td>7.5%<br>(5.4%)</td><td>CV</td><td>9.0%<br>(7.4%)</td><td>CV</td><td>10.6%<br>(9.6%)</td></tr><tr><td>∂P<sub>f</sub>/∂δ<sub>Z<sub>4</sub></sub></td><td>2.49</td><td>2.48</td><td>2.50</td><td>2.50</td><td></td><td></td></tr><tr><td></td><td>CV</td><td>10.7%<br>(8.1%)</td><td>CV</td><td>13.3%<br>(11.6%)</td><td>CV</td><td>19.1%<br>(18.2%)</td></tr><tr><td>∂P<sub>f</sub>/∂μ<sub>Z<sub>5</sub></sub></td><td>-1.82 × 10<sup>-12</sup></td><td>-1.82 × 10<sup>-12</sup></td><td>-1.81 × 10<sup>-12</sup></td><td>-1.83 × 10<sup>-12</sup></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>7.5%<br>(5.2%)</td><td>CV</td><td>9.1%<br>(7.2%)</td><td>CV</td><td>16.5%<br>(9.6%)</td></tr><tr><td>∂P<sub>f</sub>/∂δ<sub>Z<sub>5</sub></sub></td><td>2.01 × 10<sup>-12</sup></td><td>2.02 × 10<sup>-12</sup></td><td>1.98 × 10<sup>-12</sup></td><td>2.02 × 10<sup>-12</sup></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>10.2%<br>(7.3%)</td><td>CV</td><td>13.2%<br>(10.5%)</td><td>CV</td><td>48.5%<br>(15.3%)</td></tr><tr><td>∂P<sub>f</sub>/∂μ<sub>Z<sub>6</sub></sub></td><td>-3.74 × 10<sup>-12</sup></td><td>-3.75 × 10<sup>-12</sup></td><td>-3.74 × 10<sup>-12</sup></td><td>-3.77 × 10<sup>-12</sup></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>7.4%<br>(5.3%)</td><td>CV</td><td>9.3%<br>(7.3%)</td><td>CV</td><td>16.4%<br>(9.2%)</td></tr><tr><td>∂P<sub>f</sub>/∂δ<sub>Z<sub>6</sub></sub></td><td>1.97 × 10<sup>-12</sup></td><td>2.00 × 10<sup>-12</sup></td><td>2.32 × 10<sup>-12</sup></td><td>2.02 × 10<sup>-12</sup></td><td></td><td></td></tr><tr><td></td><td>CV</td><td>15.2%<br>(13.2%)</td><td>CV</td><td>45.1%<br>(20.9%)</td><td>CV</td><td>320%<br>(41.2%)</td></tr></tbody></table>

From Tables 3 and 4, the sensitivity estimates of both the proposed approach (omitting the estimates with respect to $\delta_{Z_2}$ and $\delta_{Z_6}$ in red) and the Weak approach are very close to the reference values. For the three degrees selected, the value of the estimate is quite constant, except for the derivatives, with respect to $\delta_{Z_2}$ and $\delta_{Z_6}$. Therefore, a polynomial of degree 2 is sufficient for correctly approximating the coefficient of order zero of the Taylor series expansion of the functions $\sigma \mapsto \partial P_f(\mathbf{s}, \sigma) / \partial s_\ell$ for $s_\ell \neq \delta_{Z_2}$ and $s_\ell \neq \delta_{Z_6}$. The sensitivity estimates obtained with $2k = 2$ have a smaller CV than the ones obtained with the Weak approach for an equivalent bias. Consequently, for this application, the proposed method is an improvement of the Weak approach for 10 out of 12 derivatives.

### 5.2.3. Focus on the Derivatives with Respect to $\delta_{Z_2}$ and $\delta_{Z_6}$

The sensitivity estimates with respect to the distribution parameters $\delta_{Z_2}$ and $\delta_{Z_6}$ obtained with the proposed approach have a larger CV than those obtained with the Weak approach, no matter the polynomial degree. For $2k \le 4$, the estimates have a small bias and a moderate CV. However, when $2k = 6$, the bias is higher and the CV considerably increases, reaching values above 300%. It should also be noted that the theoretical CV estimates are meaningless, as they all have an empirical CV above 100% (not displayed here) for the three degrees.

These poor results can be explained by the poor quality of the estimates $V_l$ and their variance estimates, but it is mostly due to the phenomenon of polynomial overfitting. As previously underlined, polynomial regressions of higher degrees are harder to perform and the resulting polynomial is more flexible than with smaller degrees. For these two distribution parameters, the noise of each $V_l$ is quite high and the regression interval can be particularly narrow, as illustrated in Figure 4.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Illustration of the proposed method with $2k = 6$ for the roof truss application; (**a–c**) illustrate the method with respect to $s_ℓ = δ_{Z_2}$, while (**d–f**) illustrate the method with respect to $s_ℓ = δ_{Z_6}$. The magenta vertical dashed lines represent the regression interval bounds; (**a,d**) represent the evolution of the IS estimate according to $σ$ (in navy blue); (**b,e**) show the evolution of the theoretical CV estimate of the IS estimate and the value of the CV$_{target}$; (**c,f**) represent the estimates $(V_l)_{l=1,...,5}$ (displayed with navy blue crosses) obtained with bootstrap. The resulting polynomial $\hat{P}_6$ is displayed in green; the failure probability sensitivity estimate value plotted is displayed in a dotted green line.</figcaption></figure>

The resulting polynomials shown in Figure 4c,f have much more fluctuations in the regression interval than the original IS estimates shown in Figure 4a,d. For these parameters, the polynomial regression does not smooth the noise of the IS estimates $V_l$. The polynomials are too flexible and they excessively take into account the noise of each IS estimate; it is the phenomenon of overfitting. This phenomenon does not occur for lower degrees, as the polynomials are then less flexible. Consequently, when $2k = 6$, the polynomial regression can lead to sensitivity estimates that are very biased, which explains the large CVs.

We should note that in other studies that focused on this roof truss application, these two distribution parameters always led to sensitivity estimates with significantly higher CVs compared to the other distribution parameters [4,22,23].

**6. Conclusions**

In this article, we present a new method to compute the local sensitivity of a failure probability, with respect to design parameters or distribution parameters, based on a heteroscedastic polynomial regression. This approach is inspired by the recent Weak approach framework and is presented as an improvement of the latter. The main innovation of the proposed approach is to express the sensitivity estimate as the constant coefficient of a Taylor series expansion, which can be recovered with a polynomial regression. The proposed approach can be applied to various simulation methods and is presented here with IS. Moreover, this approach is independent of the dimension of the system, the distribution of the inputs, along with the shape of the failure domain. Indeed, after a variable change in the integral of interest, only the scalar response of the system matters.
One of the main outlooks of the proposed method is to improve the stability of the resulting sensitivity estimate. Indeed, in several applications, the mean value of the estimates is very close to the reference value but the CV can be slightly too large for the method to be a definitive upgrade over the Weak approach. Another interesting outlook of the proposed approach is to obtain a better estimation of the theoretical variance of the sensitivity estimate. Indeed, with the Weak approach, the theoretical variance is already available using the formula in Equation (23) for IS. The availability of an accurate estimation of the theoretical variance is of great interest.

**Author Contributions:** Conceptualization, M.C., J.M. and S.D.; Methodology, M.C., J.M. and S.D.; Software, M.C.; Validation, M.C.; Formal analysis, M.C., J.M. and S.D.; Investigation, M.C.; Writing—original draft, M.C.; Writing—review & editing, M.C., J.M. and S.D.; Supervision, J.M. and S.D.; Funding acquisition, J.M. and S.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Chabridon, V. Reliability-Oriented Sensitivity Analysis under Probabilistic Model Uncertainty—Application to Aerospace Systems. Ph.D. Thesis, Université Clermont Auvergne, Clermont-Ferrand, France, 2018.

2. Millwater, H.R.; Wieland, D.H. Probabilistic sensitivity-based ranking of damage tolerance analysis elements. *J. Aircr.* **2010**, *47*, 161–171. [CrossRef]

3. Moustapha, M.; Sudret, B. Surrogate-assisted reliability-based design optimization: A survey and a unified modular framework. *Struct. Multidiscip. Optim.* **2019**, *60*, 2157–2176. [CrossRef]

4. Papaioannou, I.; Breitung, K.; Straub, D. Reliability sensitivity estimation with sequential importance sampling. *Struct. Saf.* **2018**, *75*, 24–34. [CrossRef]

5. Torii, A.J.; Novotny, A.A. A priori error estimates for local reliability-based sensitivity analysis with Monte Carlo Simulation. *Reliab. Eng. Syst. Saf.* **2021**, *213*, 107749. [CrossRef]

6. Morio, J.; Balesdent, M. *Estimation of Rare Event Probabilities in Complex Aerospace and Other Systems: A Practical Approach*; Woodhead Publishing: Sawston, UK, 2015.

7. Bucklew, J.A. *Introduction to Rare Event Simulation*; Springer: Berlin/Heidelberg, Germany, 2004; Volume 5.

8. Torii, A.J. On sampling-based schemes for probability of failure sensitivity analysis. *Probabilistic Eng. Mech.* **2020**, *62*, 103099. [CrossRef]

9. Robert, C.P.; Casella, G. *Monte Carlo Statistical Methods*; Springer: Berlin/Heidelberg, Germany, 1999; Volume 2.

10. Breitung, K.W. *Asymptotic Approximations for Probability Integrals*; Springer: Berlin/Heidelberg, Germany, 2006.

11. Papaioannou, I.; Breitung, K.; Straub, D. Reliability sensitivity analysis with Monte Carlo methods. In *Proceedings of the ICOSSAR 2013, New York, NY, USA, 16–20 June 2013*.

12. Lacaze, S.; Brevault, L.; Missoum, S.; Balesdent, M. Probability of failure sensitivity with respect to decision variables. *Struct. Multidiscip. Optim.* **2015**, *52*, 375–381. [CrossRef]

13. Ang, A.H.S.; Tang, W.H. *Probability Concepts in Engineering Planning and Design: Emphasis on Application to Civil and Environmental Engineering*; Wiley: Hoboken, NJ, USA, 2007.

14. Khan, R.A. Approximation for the expectation of a function of the sample mean. *Statistics* **2004**, *38*, 117–122. [CrossRef]

15. Biau, G.; Mason, D.M. High-dimensional *p*-norms. In *Mathematical Statistics and Limit Theorems: Festschrift in Honour of Paul Deheuvels*; Springer: Cham, Switzerland, 2015; pp. 21–40.

16. Apostol, T. *Calculus*; John WiJey & Sons, Inc.: Hoboken, NJ, USA, 1967.

17. Watson, G.S. Linear least squares regression. *Ann. Math. Stat.* **1967**, *38*, 1679–1699. [CrossRef]

18. Lai, T.L.; Robbins, H.; Wei, C.Z. Strong consistency of least squares estimates in multiple regression II. *J. Multivar. Anal.* **1979**, *9*, 343–361. [CrossRef]

19. Taboga, M. Generalized least squares. In *Lectures on Probability Theory and Mathematical Statistics*; John Wiley & Sons: Hoboken, NJ, USA, 2021.

20. Horowitz, J.L. The bootstrap. In *Handbook of Econometrics*; Elsevier: Amsterdam, The Netherlands, 2001; Volume 5, pp. 3159–3228.

21. Stoica, P.; Selen, Y. Model-order selection: A review of information criterion rules. *IEEE Signal Process. Mag.* **2004**, *21*, 36–47. [CrossRef]

22. Proppe, C. Local reliability based sensitivity analysis with the moving particles method. *Reliab. Eng. Syst. Saf.* **2021**, *207*, 107269. [CrossRef]
23. Song, S.; Lu, Z.; Qiao, H. Subset simulation for structural reliability sensitivity analysis. *Reliab. Eng. Syst. Saf.* **2009**, *94*, 658–665.
[CrossRef]

24. Morio, J. Extreme quantile estimation with nonparametric adaptive importance sampling. *Simul. Model. Pract. Theory* **2012**, *27*, 76–89.
[CrossRef]

25. Papaioannou, I.; Geyer, S.; Straub, D. Improved cross entropy-based importance sampling with a flexible mixture model. *Reliab. Eng. Syst. Saf.* **2019**, *191*, 106564.
[CrossRef]

26. Yang, R.; Gu, L. Experience with approximate reliability-based optimization methods. *Struct. Multidiscip. Optim.* **2004**, *26*, 152–159.
[CrossRef]

27. Madsen, H.O.; Krenk, S.; Lind, N.C. *Methods of Structural Safety*; Courier Corporation: Chelmsford, MA, USA, 2006.

28. Rubinstein, R.Y. The score function approach for sensitivity analysis of computer simulation models. *Math. Comput. Simul.* **1986**, *28*, 351–379.
[CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.