# Distributions of Outputs Given Subsets of Inputs and Dependent Generalized Sensitivity Indices

Matieyendou Lamboni

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

## To cite this version:

Matieyendou Lamboni. Distributions of Outputs Given Subsets of Inputs and Dependent Generalized Sensitivity Indices. Mathematics, 2025, 13 (5), pp.766. 10.3390/math13050766. hal-05083368

HAL Id: hal-05083368

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

Submitted on 24 May 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Distributions of Outputs Given Subsets of Inputs and Dependent Generalized Sensitivity Indices

Matieyendou Lamboni <fnref n="1" />,<fnref n="2" />

**Abstract:** Better understanding mathematical and numerical models often requires investigating the impacts of inputs on the model outputs, as well as interactions. Quantifying such effects for models with non-independent input variables (NIVs) relies on conditional distributions of the outputs given every subset of inputs. In this paper, by firstly providing additional dependency models of NIVs, functional outputs are composed by dependency models (yielding equivalent representations of outputs) to derive distributions of outputs conditional on inputs. We then provide an algorithm for selecting the necessary and sufficient equivalent representations that allow for obtaining all the conditional distributions of outputs given every subset of inputs, and for assessing the main, total, and interaction effects (i.e., indices) of every subset of NIVs. Unbiased estimators of covariances of sensitivity functionals and consistent estimators of such indices are derived by distinguishing the case of the multivariate and/or functional outputs, including dynamic models. Finally, analytical results and numerical results are provided, including an illustration based on a dynamic model.

**Keywords:** computer models; copulas; dependency models; multivariate sensitivity analysis; non-independent and discrete variables

MSC: 49Q12; 60E05; 62F12

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

check for

updates

Academic Editor: Manuel Alberto M. Ferreira

Received: 31 December 2024

Revised: 20 February 2025

Accepted: 23 February 2025

Published: 26 February 2025

**Citation:** Lamboni, M. Distributions of Outputs Given Subsets of Inputs and Dependent Generalized Sensitivity Indices. Mathematics **2025**, *13*, 766. https://doi.org/10.3390/math13050766

**Copyright:** © 2025 by the author. 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

Mathematical and numerical models and/or simulators are being increasingly developed and used for supporting decision-making in different scientific fields, such as engineering, environment, agronomy, and biology. To better understand such complex models or to design emulators, it is worth investigating the impacts of inputs on the model outputs, including the interactions among inputs. Uncertainty quantification and global sensitivity analysis are dedicated to addressing such issues. While variance-based methods (see, e.g., Sobol' indices [1] for real-valued functions and generalized sensitivity indices [2–7] for dynamic and multivariate outputs) are well established for independent inputs, three are cases where the input variables are correlated or dependent.

Non-independent variables (NIVs) arise when two or more variables do not vary freely and are widely encountered in different scientific fields, such as in data analysis, quantitative risk analysis, and inverse problems. Performing uncertainty quantification and variance-based sensitivity analysis of computer and/or mathematical models in the presence of dependent and/or correlated input variables (i.e., NIVs) still remains a challenge when one is interested in assessing the contributions of any subset of input variables

<footnote n="1">¹ Department DFR-ST, University of Guyane, 97346 Cayenne, France; matieyendou.lamboni@gmail.com or matieyendou.lamboni@univ-guyane.fr</footnote>

<footnote n="2">² 228-UMR Espace-Dev, University of Guyane, University of Réunion, IRD, University of Montpellier, 34090 Montpellier, France</footnote>
and their interactions over all the outputs. Indeed, the dependency structures inferred by
the constraints imposed on model outputs and/or inputs, and the dependency structures
of the initial distributions of inputs may have significant impacts on the results of Sobol'
indices and generalized sensitivity indices. A number of papers report inconsistent ranking
of inputs using the above indices in the presence of NIVs (e.g., [8–11]).

In the presence of NIVs, most existing variance-based methodologies provide the same first-order index of one input or one group of dependent inputs (e.g., [8–13]). Some of them provide the total indices of inputs, and there can be cases in which the first-order index is greater than the total index. To be able to rank input variables, the works in [14,15] introduced a methodology that ensures that the first-order index is always less than the total index for every single input and for real-valued functions. In the same sense, the recent works in [16–18] provided in-depth approaches for quantifying the effects of each single input and some subsets of inputs (but not all the subsets) by making use of dependency models of NIVs. Note that dependency models provide functions that model the dependency structures of NIVs ([16,19,20]), and such recent approaches rely on different DMs of the same NIVs for defining and computing dependent generalized sensitivity indices (dGSIs) of some subsets of inputs.

In this paper, we propose a new methodology for (i) deriving all the conditional distributions of outputs given every subset of NIVs and (ii) assessing the impacts of every subset of uncertain NIVs and their interactions on the model outputs by making use of the necessary and sufficient equivalent representations of the model of interest. Equivalent representations are obtained by composing the model outputs by the necessary and sufficient DMs, leading to reducing the number of model runs for computing such quantities of interest.

This paper is organized as follows: in Section 2, we provide additional and generic DMs of NIVs so as to cover more distribution functions, and we investigate the DM transforms, which avoids searching dependency functions for some distributions by making use of known DMs already derived in [16,19,20]. We also provide computational DMs when analytical distributions of the input variables are not available, such as the resultant distribution of inputs associated with complex mathematical models under constraints. Section 3 provides the famous algorithm for selecting the necessary and sufficient DMs. It also deals with different representations of the model outputs associated with subsets of NIVs that are equivalent in distributions using the necessary and sufficient DMs. Such representations are useful for assessing the main, total, and interaction effects of any subset of NIVs in Section 4 by providing dGSIs of every subset of inputs and for the multivariate and/or functional outputs, including spatiotemporal models and dynamic models. Section 5 aims at constructing unbiased estimators of the cross-covariances of sensitivity functionals and the consistent estimators of dGSIs. We also provide the asymptotic distributions of dGSIs. Analytical results and numerical results are provided, including an illustration based on a dynamic model (see Section 6), and we conclude this work in Section 7.

General Notation

For an integer $d > 0$, we use $\mathbf{X}: = (X_1, \dots, X_d)$ for a random vector of NIVs. Given $u \subseteq \{1, \dots, d\}$, we use $\mathbf{X}_u := (X_j, \forall j \in u)$, $\mathbf{X}_{\sim u} := (X_j, \forall j \in \{1, \dots, d\} \setminus u)$ and $|u|$ for its cardinality, leading to the partition $\mathbf{X} = (\mathbf{X}_u, \mathbf{X}_{\sim u})$. We also use $\mathbf{Z} \stackrel{d}{=} \mathbf{X}$ to say that $\mathbf{Z}$ and $\mathbf{X}$ have the same cumulative distribution function (CDF). For $\mathbf{a} \in \mathbb{R}^n$, we use $||\mathbf{a}||_{L^2}$ for the Euclidean norm. For a matrix $\mathbf{\Sigma} \in \mathbb{R}^{n \times n}$, we use $\text{Tr}(\mathbf{\Sigma})$ for the trace of $\mathbf{\Sigma}$, and $||\mathbf{\Sigma}||_F := \sqrt{\text{Tr}(\mathbf{\Sigma}\mathbf{\Sigma}^T)}$ for the Frobenius norm of $\mathbf{\Sigma}$. We use $\mathbb{E}[\cdot]$ for the expectation operator and $\mathbb{V}[\cdot]$ for the variance-covariance operator.
## 2. Dependency Functions of Non-Independent Random Variables

This section provides generic DMs of NIVs and some transformations of such DMs. Formally, the inputs X have F as the CDF and C as the copula, that is, $F(x) = C(F_1(x_1), \dots, F_d(x_d))$ with $F_{x_j}$ or $F_j$ the marginal CDF of $X_j$, $j = 1, \dots, d$ and x a sample value of X. We use $F_j^{\leftarrow}$ for the generalized inverse of $F_j$, and $F_{j|k}$ for the distribution of $X_j$ conditional on $X_k$ for all $j, k \in \{1, \dots, d\}$ and $j \neq k$. For a given discrete variable $X_i$, let us consider the distribution transform of $X_i$ given by $\tau_{F_i}(x_i, \lambda_i) = P(X_i < x_i) + \lambda_i P(X_i = x_i)$ with $\lambda_i \in [0, 1]$. Such a distribution transform ensures that for all $U_i \sim U(0, 1)$ [21–24]

$$V_i := \tau_{F_i}(X_i, U_i) \sim U(0, 1), \quad X_i \stackrel{d}{=} F_i^{\leftarrow}(V_i) \quad \text{a.s.} \quad (1)$$

For continuous distributions, the former term of Equation (1) comes down to the Rosenblatt transform [25], and the latter term is equivalent to the inverse of the Rosenblatt transform. Denote with $(w_1, \dots, w_{d-1})$ an arbitrary permutation of $\{1, \dots, d\} \setminus \{j\}$, $\mathbf{X}_{\sim j} := (X_{w_1}, \dots, X_{w_{d-1}})$ and $\mathbf{Z} \sim U(0, 1)^{d-1}$. A generic DM of X is given by [16,20]

$$\mathbf{X}_{\sim j} \stackrel{d}{=} r_j(\mathbf{X}_j, \mathbf{Z}) = (r_{w_1}(\mathbf{X}_j, Z_{w_1}), r_{w_2}(\mathbf{X}_j, Z_{w_1}, Z_{w_2}), \dots, r_{w_{d-1}}(\mathbf{X}_j, Z_{w_1}, \dots, Z_{w_p})), \quad (2)$$

where $\mathbf{Z}$ is a random vector of independent variables, and it is independent of $X_j$ as well. It is worth noting that the dependency function $r_j := (r_{w_1}, \dots, r_{w_{d-1}}) : \mathbb{R}^d \to \mathbb{R}^{d-1}$ is not unique. Indeed, one may replace, in Equation, (2) $Z_{w_i}$ with $F_{w_i}(V_{w_i})$ for any continuous variable $V_{w_i} \sim F_{w_i}$ or with $\tau_{G_{w_i}}(W_{w_i}, U_{w_i})$ for any discrete variable $W_{w_i} \sim G_{w_i}$, provided that such variables are independent of $X_j$. But, DMs are uniquely defined once the marginal CDFs of $\mathbf{Z}$ are prescribed.

**Proposition 1** ([18]). Let $r_j$ be a dependency function; $0 < p \le d-1$; $v := (j, w_1, \dots, w_p)$ be a vector and $\mathbf{Z}_{\sim v} \sim U(0, 1)^{d-p-1}$. Then,

$$\mathbf{X}_v := (X_j, X_{w_1}, \dots, X_{w_p}) \stackrel{d}{=} (X_j, r_{w_1}(X_j, Z_{w_1}), \dots, r_{w_p}(X_j, Z_{w_1}, \dots, Z_{w_p})). \quad (3)$$

Moreover, there exists $r_v : \mathbb{R}^d \to \mathbb{R}^{d-p-1}$ such that

$$\mathbf{X}_{\sim v} := (X_{w_{p+1}}, \dots, X_{w_{d-1}}) \stackrel{d}{=} r_v(\mathbf{X}_v, \mathbf{Z}_{\sim u}). \quad (4)$$

### 2.1. Distribution-Based and Copula-Based Expressions of Dependency Models

The distribution-based dependency model of X is given by [16,19]

$$\begin{align}
X_{w_1} &= r_{w_1}(X_j, Z_{w_1}) = F_{w_1|j}^{\leftarrow}(Z_{w_1} | X_j) \nonumber \\
X_{w_{d-1}} &= F_{w_{d-1}|j,w_1,\dots,w_{d-2}}^{\leftarrow}(Z_{w_{d-1}} | X_j, r_{w_1}(X_j, Z_{w_1}), \dots, r_{w_{d-2}}(X_j, Z_{w_1}, \dots, Z_{w_{d-2}})) \tag{5}
\end{align}$$

Likewise, the copula-based dependency function is of interest to master all joint distributions having the same copula C, regardless of their marginal CDFs ([23,26]). It allows one to use the same DMs for the class of distributions sharing the same copula. For independent variables $U \sim U(0, 1)^d$, $\mathbf{Z}$, and $X_j$, the copula-based DM is given by

$$\begin{align}
X_{w_1} &= r_{w_1}(X_j, Z_{w_1}, U_j) = F_{w_1|j}^{\leftarrow}(C_{w_1|j}^{\leftarrow}(Z_{w_1} | \tau_{F_j}(X_j, U_j))) \nonumber \\
X_{w_{d-1}} &= F_{w_{d-1}|j,(\sim w_{d-1})}^{\leftarrow}(C_{w_{d-1}|j,(\sim w_{d-1})}^{\leftarrow}(Z_{w_{d-1}} | \tau_{F_j}(X_j, U_j), \tau_{F_{w_1}}(r_{w_1}(X_j, Z_{w_1}, U_j), U_{w_1}), \dots)) \tag{6}
\end{align}$$

by making use of the conditional sampling algorithm based on copulas [26,27]. Note that for Gaussian copulas and continuous marginal CDFs of inputs, DMs in (6) have a particular
form provided in [16]. Lemma 1 extends such DMs to cope with discrete CDFs as well. To that end, denote with $C^{\text{Gauss}}(U_1, \dots, U_d, \mathcal{R})$ the Gauss copula having $\mathcal{R}$ as the correlation matrix; $\mathcal{L}$ the Cholesky factor of $\mathcal{R}$ (i.e., $\mathcal{R} = \mathcal{LL}^T$); $\Phi$ the CDF of the standard Gaussian variable; and $\mathcal{I}$ the identity matrix.

**Lemma 1.** Let $X_j, Z \sim \mathcal{N}_{d-1}(0, \mathcal{I})$ and $U_j \sim \mathcal{U}(0, 1)$ be independent variables.

If $(X_j, \mathbf{X}_{\sim j})$ has the copula $C^{\text{Gauss}}(U_j, U_{w_1}, \dots, U_{w_{d-1}}, \mathcal{R})$, then

$$ \mathbf{X}_{\sim j} \stackrel{d}{=} r_j(\mathbf{X}_j, \mathbf{Z}, U_j) = \begin{bmatrix} X_{w_1} = F_{w_1}^{\leftarrow}[\Phi(Y_{w_1})] \\ \vdots \\ X_{w_{d-1}} = F_{w_{d-1}}^{\leftarrow}[\Phi(Y_{w_{d-1}})] \end{bmatrix}, \quad (7) $$

with

$$ \begin{bmatrix} Y_j \\ \mathbf{Y}_{\sim j} \end{bmatrix} := \mathcal{L} \begin{bmatrix} \Phi^{-1}(\tau_{F_j}(X_j, U_j)) \\ \mathbf{Z} \end{bmatrix}. $$

**Proof.** See Appendix A. □

Likewise, to provide DMs for the Student copulas (Lemma 2), we use $t(\nu, 0, 1)$ for the standard $t$-distribution with $\nu$ degrees of freedom and $T_\nu$ for its CDF.

**Lemma 2.** Let $X_j, \{Z_{w_i} \sim t(\nu + i, 0, 1)\}_{i=1}^{d-1}$ and $U_j$ be independent variables.

If $(X_j, \mathbf{X}_{\sim j})$ has the Student copula $C^{\text{St}}(U_j, U_{w_1}, \dots, U_{w_{d-1}}, \nu, \mathcal{R})$, then

$$ \mathbf{X}_{\sim j} \stackrel{d}{=} r_j(\mathbf{X}_j, \mathbf{Z}, U_j) = \begin{bmatrix} X_{w_1} = F_{w_1}^{\leftarrow}[T_\nu(Y_{w_1})] \\ \vdots \\ X_{w_{d-1}} = F_{w_{d-1}}^{\leftarrow}[T_\nu(Y_{w_{d-1}})] \end{bmatrix}, \quad (8) $$

with

$$ \begin{bmatrix} Y_j \\ \mathbf{Y}_{\sim j} \end{bmatrix} := \mathcal{L} \begin{bmatrix} \frac{T_\nu^{-1}(\tau_{F_j}(X_j, U_j))}{\sqrt{\frac{\nu + (T_\nu^{-1}(\tau_{F_j}(X_j, U_j)))^2}{\nu+1}}} Z_{w_1} \\ \vdots \\ \frac{\sqrt{\frac{(\nu + (T_\nu^{-1}(\tau_{F_j}(X_j, U_j)))^2}{\nu+1} \prod_{k=1}^{d-2} (\nu+k+Z_{w_k}^2)} \prod_{k=1}^{d-1} (\nu+k) Z_{w_{d-1}}}{\prod_{k=1}^{d-1} (\nu+k)} \end{bmatrix}. $$

**Proof.** See Appendix B. □

Note that for every continuous variable $X_j$, we have to replace $\tau_{F_j}(X_j, U_j)$ with $F_j(X_j)$ in Lemmas 1 and 2. For discrete variables, we have to include additional and independent uniformly distributed variables $U$ in copula-based DMs (see Equation (6)).

## 2.2. Empirical and Computational Dependency Models

This section deals with the derivation of DMs for unknown distributions of inputs, such as distributions obtained by imposing constraints on the initial inputs or outputs. Formally, given a function $c: \mathbb{R}^d \to \mathbb{R}^n$ and a domain of interest $D$, we are interested in deriving a dependency function of a random vector defined by

$$ \mathbf{X}^c \stackrel{d}{=} \{\mathbf{X} \sim F: c(\mathbf{X}) \in D\}. \quad (9) $$
While we are able to derive the analytical distribution of $X^c$ and its dependency function for some distributions and constraints (see [16,20]), we have to estimate such dependency functions in general. Using Equation (9), we can generate a sample of $X^c$, that is, $X_1^c, \dots, X_m^c$ and a pseudo-sample from the copula $C$ of $X^c$, that is, $\hat{F}_1(X_{i,1}^c), \dots, \hat{F}_d(X_{i,d}^c)$ with $i = 1, \dots, m$ for continuous variables, where $\hat{F}_j$ is an estimator of $F_j$. In general, we use $\tau_{\hat{F}_1}(X_{i,1}^c, U_1), \dots, \tau_{\hat{F}_d}(X_{i,d}^c, U_d)$.

With such samples, we consider two main ways to derive the empirical dependency functions. Firstly, we fit a distribution to such observations and then derive dependency functions using results from Section 2. To that end, there are numerous papers about fitting a distribution to data. For instance, direct methods for estimating densities and distributions can be found in [28–31], and the copula-based methods for modeling distributions are provided in [23,27,32–35].

Secondly, we derive empirical dependency functions using the estimators of the conditional quantile functions [36–43]. Formally, consider the loss function of Koenker and Bassett [36] given by $L(x, u) = x(u - I_{\{x<0\}})$ with $u \in [0, 1]$ and $I_{\{x<0\}}$ the indicator function. A dependency function can be written as follows:

$$r_j(X_j^c, Z_{w_1}) := \arg \min_{f \in \mathcal{F}} \mathbb{E} [L(X_{w_1}^c - f(X_j^c), Z_{w_1}) | X_j^c, Z_{w_1}],$$

where $\mathcal{F}$ is a class of smooth functions, and $Z_{w_1} \sim U(0, 1)$ is independent of $(X_j^c, X_{w_1}^c)$. Using the sample of $(X_j^c, X_{w_1}^c)$, the M-estimator of a dependency function is given by ([43], Lemma 3)

$$\hat{r}_j(X_j^c, Z_{w_1}) := \arg \min_{f \in \mathcal{H}} \sum_{i=1}^{m} L(X_{i,w_1}^c - f(X_{i,j}^c), Z_{w_1}) + \frac{\lambda}{2} \|f(X_j^c) - b\|_{\mathcal{H}}^2, \quad (10)$$

where $\lambda \in \mathbb{R}_+$ is a bandwidth, $\|\cdot\|_{\mathcal{H}}$ is a given norm, and $b \in \mathbb{R}$.

### 3. Equivalent Representations of Functional Outputs

This section formalizes different representations of complex models with NIVs $X$. Formally, given $\Theta \subseteq \mathbb{R}$, $n \in \mathbb{N}^*$, consider a function $f: \mathbb{R}^d \times \Theta \to \mathbb{R}^n$ with random evaluations, that is, $f(X, \theta) \in \mathbb{R}^n$ with $\theta \in \Theta$. It may represent any multivariate and functional outputs. When $\Theta = \{\theta_0\}$ with $\theta_0 \in \mathbb{R}$, we obtain a class of vector-valued functions. In what follows, $X$ is organized as follows:

(O): $X := (X_1, \dots, X_d)$ consists of $K$ independent random vector(s), that is, $X = (X_{\pi_1}, \dots, X_{\pi_K})$ where $\pi_1, \dots, \pi_K$ are the sets that form a partition of $\{1, \dots, d\}$; $X_{\pi_{k_1}} := (X_j, \forall j \in \pi_{k_1})$ is independent of $X_{\pi_{k_2}} := (X_j, \forall j \in \pi_{k_2})$ for all $k_1, k_2 \in \{1, \dots, K\}$ and $k_1 \neq k_2$. Without loss of generality, we use $X_{\pi_1}$ for a random vector of $d_1 \ge 0$ independent variable(s); $X_{\pi_k}$ with $k \ge 2$ for a random vector of $d_k \ge 2$ NIVs. Note that $|\pi_k| = d_k$.

Denote with $w_k := (w_{1,k}, \dots, w_{d_k,k})$ an arbitrary permutation of $\pi_k$ and $w_{\pi_{1,k}} := (w_{2,k}, \dots, w_{d_k,k})$. Thus, $w_{i,k} \in \pi_k, \forall i \in \{1, \dots, d_k\}$. If we use $s := \{w_{1,2}, \dots, w_{1,K}\}$, then $X_s := (X_{w_{1,2}}, \dots, X_{w_{1,K}})$ contains $K-1$ independent variables and $X_{\pi_{\cdot s}} := (X_{w_{2,k}}, \dots, X_{w_{d_k,k}})$ with $X_{\pi_{\cdot s}} := (X_{w_{2,k}}, \dots, X_{w_{d_k,k}})$, $k = 2, \dots, K$.

Using the DMs of Section 2, we can write

$$X_{\pi_{\cdot s}} = (X_{w_{1,2}} \stackrel{d}{=} r_{w_{1,2}}(X_{w_{1,2}}, Z_{w_{1,2}}), \dots, X_{w_{\pi_{\cdot K}} \stackrel{d}{=} r_{w_{\pi_{\cdot K}}}(X_{w_{\pi_{\cdot K}}}, Z_{w_{\pi_{\cdot K}}})) ; \quad (11)$$

where $Z_{w_{\pi_{\cdot K}} \stackrel{d}{=} (Z_{w_{2,k}}, \dots, Z_{w_{d_k,k}})}$ is a vector of independent variables, $X_{w_{1,k}}$ is independent of $Z_{w_{\pi_{\cdot K}} \stackrel{d}{=} (Z_{w_{1,2}}, \dots, Z_{w_{1,K}})}$, $k = 2, \dots, K$. By using $Z_{w_{\pi_{\cdot K}} \stackrel{d}{=} (Z_{w_{2,k}}, \dots, Z_{w_{d_k,k}})}$, we can see that
$(X_{\pi_1}, X_s, Z_{w_{\sim 1}})$ contains only independent variables. Compiling the $K-1$ DMs given by (11) in one function, that is, $r_s: \mathbb{R}^{d-d_1} \to \mathbb{R}^{d-d_1-K+1}$ yields

$$X_{\sim s} \stackrel{d}{=} r_s(X_s, Z_{w_{\sim 1}}) = (r_{w_{1,2}}(X_{w_{1,2}}, Z_{w_{\sim 1,2}}), \dots, r_{w_{1,K}}(X_{w_{1,K}}, Z_{w_{\sim 1,K}})),$$

and we have the following partition:

$$(X_2, \dots, X_K) = (X_s, X_{\sim s}), \quad X \stackrel{d}{=} (X_{\pi_1}, X_s, r_s(X_s, Z_{w_{\sim 1}})). \quad (12)$$

Composing $f(\cdot, \theta)$ by (12) yields

$$g(X_{\pi_1}, X_s, Z_{w_{\sim 1}}, \theta) := f(X_{\pi_1}, X_s, r_s(X_s, Z_{w_{\sim 1}}), \theta), \quad (13)$$

and Lemma 3 provides useful properties of $g$ linked to $f$. For a given integer $0 \le p_k \le d_k$, consider the vectors $v_k := (w_{2,k}, \dots, w_{p_k,k})$ and $u_k := (w_{1,k}, v_k)$ with $k = 2, \dots, K$. By definition, when $p_k = 0$, we have $v_k = u_k = \emptyset$, and when $p_k = 1$, $v_k = \emptyset$ and $u_k = w_{1,k}$.

**Lemma 3.** Let $u_0 \subseteq \pi_1, \{k_1, \dots, k_m\} \subseteq \{2, \dots, K\}, \{w_{1,k_1}, \dots, w_{1,k_m}\} \subseteq s$. Then,

$$f(\mathbf{X}, \theta) | \mathbf{X}_{u_0}, \mathbf{X}_{u_{k_1}}, \dots, \mathbf{X}_{u_{k_m}} \stackrel{d}{=} g(\mathbf{X}_{\pi_1}, \mathbf{X}_s, \mathbf{Z}_{w_{\sim 1}}, \theta) | \mathbf{X}_{u_0}, \mathbf{X}_{w_{1,k_1}}, \mathbf{Z}_{v_{k_1}}, \dots, \mathbf{X}_{w_{1,k_m}}, \mathbf{Z}_{v_{k_m}}.$$

**Proof.** See Appendix C. □

It comes out from Lemma 3 that the distribution of $f(\mathbf{X}, \theta)$ given the inputs $\mathbf{X}_u := (\mathbf{X}_{u_0}, \mathbf{X}_{u_{k_1}}, \dots, \mathbf{X}_{u_{k_m}})$ is equal (in distribution) to the distribution of $g(\mathbf{X}_{\pi_1}, \mathbf{X}_s, \mathbf{Z}_{w_{\sim 1}}, \theta)$ conditional on $(\mathbf{X}_{u_0}, \mathbf{X}_{w_{1,k_1}}, \mathbf{Z}_{v_{k_1}}, \dots, \mathbf{X}_{w_{1,k_m}}, \mathbf{Z}_{v_{k_m}})$. Thus, we are able to assess the effect of $\mathbf{X}_u$ on $f(\mathbf{X}, \theta)$ using $g(\mathbf{X}_{\pi_1}, \mathbf{X}_s, \mathbf{Z}_{w_{\sim 1}}, \theta)$ and $(\mathbf{X}_{u_0}, \mathbf{X}_{w_{1,k_1}}, \mathbf{Z}_{v_{k_1}}, \dots, \mathbf{X}_{w_{1,k_m}}, \mathbf{Z}_{v_{k_m}})$, leading to the following definition.

**Definition 1.** Consider $u \subseteq \{1, \dots, d\}$, and $f$ and $g$ are given by Equation (13).

A function $g$ is said to be an equivalent representation of $f$ regarding the input(s) $\mathbf{X}_u$ if the distribution of $f(\mathbf{X}, \theta) | \mathbf{X}_u$ can be determined using $g$ and some of its inputs.

Different equivalent representations (ERs) of $f$ (i.e., $g$) are necessary for assessing the effects of $\mathbf{X}_u$ for all $u \subseteq \{1, \dots, d\}$. For instance, a representation in Lemma 3 for a given set $s$ can be used to assess the effects of some subsets of inputs given by

$$\left\{ (\mathbf{X}_{u_0}, \mathbf{X}_{u_{k_1}}, \dots, \mathbf{X}_{u_{k_m}}) : \begin{array}{l} \forall u_0 \subseteq \pi_1, \forall \{k_1, \dots, k_m\} \subseteq \{2, \dots, K\}, \\ \forall p_{k_i} \in \{0, \dots, d_{k_i}\}, i = 1, \dots, m \end{array} \right\}.$$

Modifying $s$ and $v_k, k = 2, \dots, K$ leads to another representation of $f$, which allows for assessing other inputs' effects, such as the effects of $X_s$ with $s \notin s$. Permutations of $\pi_k$s give such modifications. Obviously, we have $\prod_{i=2}^K d_i!$ ERs of $f$ that allow for deriving the effects of all the subsets of inputs with the possibility to have some replications.

**Definition 2.** Let $u \subseteq \{1, \dots, d\}$ and $g_1 \neq g_2$ be two ERs of $f$.

The representations $g_1, g_2$ are said to be replicated representations of $f$ regarding $\mathbf{X}_u$ if $g_1, g_2$ allow for determining the distribution of $f(\mathbf{X}, \theta) | \mathbf{X}_u$.
To avoid unnecessary replicated representations and to be able to recover all the subsets of {1, ..., d} using permutations, we use Algorithm 1 for selecting the necessary and sufficient permutations (see Lemma 4). Formally, consider integers

$$j_{0,k} := \begin{cases} \frac{d_k}{2} & \text{if } d_k \text{ is even}, \\ \frac{d_k+1}{2} & \text{otherwise} \end{cases}, k = 2, \dots, K, \quad (14)$$

and the super-sets $\mathcal{A}_{j_{0,k}}$ given by

$$\mathcal{A}_{j_{0,k}} = \{u \subseteq \pi_k : |u| = j_{0,k}\}, \quad k = 2, \dots, K.$$

The set $\mathcal{A}_{j_{0,k}}$ consists of all the subsets of $\pi_k$ that contain exactly $j_{0,k}$ elements, and its cardinality is $|\mathcal{A}_{j_{0,k}}| = \binom{d_k}{j_{0,k}}$. Basically, the algorithm takes $\mathcal{A}_{j_{0,k}}$ as the input and provides the super-sets $\mathcal{B}_k, \mathcal{P}_k, \mathcal{E}_k$, which are initially empty. The first step of Algorithm 1, corresponding to $e_0 = 1$, focuses on selecting $d_k = \binom{d_k}{e_0=1}$ permutations, and on bringing different sets of the form $\{w_{1,k}\}, \dots, \{w_{1,k}, \dots, w_{j_{0,k},k}\}$ in one hand and $\{w_{d_k-j_{0,k}+1,k'}, \dots, w_{d_k,k}\}, \dots, \{w_{d_k,k}\}$ in the other hand. We repeat that process by increasing $e_0$ and bringing the sets of the forms:

$$\{w_{1,k}, \dots, w_{e_0,k}\}, \dots, \{w_{1,k}, \dots, w_{j_{0,k},k}\}; \quad \{w_{d_k-j_{0,k}+1,k'}, \dots, w_{d_k,k}\}, \dots, \{w_{d_k-e_0+1,k'}, \dots, w_{d_k,k}\};$$

until we are able to derive all the subsets of $\pi_k$, that is, until $|\mathcal{A}_{j_{0,k}}| = 0$. The formal algorithm is given as follows:

**Algorithm 1:** Construction of the sets $\mathcal{B}_k$ and $\mathcal{P}_k$ for all $k \in \{2, \dots, K\}$.

initialization: $\mathcal{B}_k \leftarrow \mathcal{P}_k \leftarrow \mathcal{E}_k \leftarrow \emptyset; i \leftarrow e_0 \leftarrow 1$;
$\mathcal{A}_{j_{0,k}} \leftarrow \{u_k \subseteq \pi_k : |u_k| = j_{0,k}\}$;
**while** $|\mathcal{A}_{j_{0,k}}| > 0$ **do**
Find a permutation $w_k = (w_{1,k}, \dots, w_{d_k,k})$ of $\pi_k$ such that
$[\{w_{1,k}, \dots, w_{j,k}\} \notin \mathcal{B}_k, j = e_0, \dots, j_{0,k}]$ and
$[\{w_{d_k-j+1,k'}, \dots, w_{d_k,k'}\} \notin \mathcal{E}, j = e_0, \dots, j_{0,k}]$ and $\{w_{1,k}, \dots, w_{j_{0,k},k}\} \in \mathcal{A}_{j_{0,k}}$;
$\mathcal{A}_{j_{0,k}} \leftarrow \mathcal{A}_{j_{0,k}} \setminus \{w_{1,k}, \dots, w_{j_{0,k},k}\}$;
$\mathcal{B}_k \leftarrow \mathcal{B}_k \cup \{\{w_{1,k}, \dots, w_{j,k}\}, j = e_0, \dots, d_k - e_0 + 1\}$;
$\mathcal{E}_k \leftarrow \mathcal{E}_k \cup \{\{w_{j,k}, \dots, w_{d_k,k}\}, j = j_{0,k} + 1, \dots, d_k - e_0 + 1\}$;
$\mathcal{P}_k \leftarrow \mathcal{P}_k \cup \mathbf{w}_k$;
$i \leftarrow i + 1$;
**if** $\binom{d_k}{e_0} < i \le \binom{d_k}{e_0+1}$ **then**
$e_0 \leftarrow e_0 + 1$;
**end**
**end**

Algorithm 1 provides the permutations selected (i.e., $\mathcal{P}_k$) that are used for deriving different DMs. For instance, if $d_2 = 3$, $\pi_2 = \{1, 2, 3\}$, the permutations (1, 2, 3) and (3, 2, 1) lead to the following DMs:

$$(X_2, X_3) \stackrel{d}{=} r_1(X_1, Z_2, Z_3); \quad (X_2, X_1) \stackrel{d}{=} r_3(X_3, Z_2, Z_1),$$
and the associated ERs, that is,

$$f(\mathbf{X}) \stackrel{d}{=} f(X_1, r_1(X_1, Z_2, Z_3)); \quad f(\mathbf{X}) \stackrel{d}{=} f(X_3, r_3(X_3, Z_2, Z_1)).$$

The set $\mathcal{P}_k$ from Algorithm 1 contains $\binom{d_k}{j_{0,k}}$ permutations. The set $\mathcal{B}_k$ is built using $\mathcal{P}_k$, and it consists of sets containing the first $t$ elements of $\mathbf{w}_k$ for any $t \in \{1, \dots, d_k\}$ and $\mathbf{w}_k \in \mathcal{P}_k$. Lemma 4 shows the properties of such sets.

**Lemma 4.** Consider $j_{0,k}$ given by (14) and $\mathcal{B}_k, \mathcal{P}_k$ given by Algorithm 1. Then,

$$\mathcal{B}_k = \{u \subseteq \pi_k : \forall |u| > 0\}; \qquad (15)$$

$$\mathcal{B}_k = \{\{w_{1,k}, \dots, w_{t,k}\}, t = 1, \dots, d_k : \forall \mathbf{w}_k \in \mathcal{P}_k\}. \qquad (16)$$

**Proof.** See Appendix D. □

Based on Lemma 4, we quantify the necessary and sufficient number of ERs of $f$ given $\mathbf{X}_u$ for all $u \subseteq \{1, \dots, d\}$ in Theorem 1. For $\mathbf{w}_k \in \mathcal{P}_k$, recall that $\mathbf{w}_{\sim 1,k} := (w_{2,k}, \dots, w_{d_k,k})$ and $|\mathcal{P}_k| = \binom{d_k}{j_{0,k}}$ with $k = 2, \dots, K$.

**Theorem 1.** Consider integers $p_2 \le d_2, \dots, p_K \le d_K$. Then, the minimum number of ERs of $f$ given $\mathbf{X}_u$ for all $u \subseteq \{1, \dots, d\}$ is

$$R_{\min} := \prod_{k=2}^{K} \binom{d_k}{j_{0,k}}. \qquad (17)$$

Such necessary and sufficient ERs are given by

$$f(\mathbf{X}, \theta) \stackrel{d}{=} g_{\ell}(\mathbf{X}_{\pi_1}, \mathbf{X}_{w_{1,2}}, \mathbf{Z}_{\mathbf{w}_{\sim 1,2}}, \dots, \mathbf{X}_{w_{1,K}}, \mathbf{Z}_{\mathbf{w}_{\sim 1,K}}, \theta), \qquad (18)$$

where $\ell := (\mathbf{w}_2, \dots, \mathbf{w}_K)$ for all $\mathbf{w}_k \in \mathcal{P}_k$ and $k = 2, \dots, K$.

**Proof.** See Appendix E. □

It is worth noting that each ER of Equation (18) shares the same distribution of $f(\mathbf{X}, \theta)$, and two different ERs should be considered independent to avoid misleading dependencies. When a function includes only independent variables, we see that $R_{\min} = 1$. Reducing $R_{\min}$ will depend on the analysis of interest one is going to perform. For instance, one ER is sufficient to determine the distribution of $f(\mathbf{X}, \theta)$ conditional on $\mathbf{X}_u$ for all

$$u \in \{\{u_0, w_{1,k}, \dots, w_{t,k}\} : \forall u_0 \subseteq \pi_1, \forall t \in \{0, \dots, d_k\}, k = 2, \dots, K\}.$$

Likewise, $R_1 = \max(d_2, \dots, d_K)$ ERs are used for assessing the main and total effects of $X_j$ for any $j \in \{1, \dots, d\}$ in [16]. It is also worth noting that the above ERs can lead to assessing the effects of other inputs or groups of inputs.

## Discussions About High-Dimensional Cases

Implementation of Algorithm 1 can be made fast for moderate values of each $d_k$, that is $d_k \le 10$ with $k = 2, \dots, K$. Since Algorithm 1 must be run separately for each $\pi_k$ with $k = 2, \dots, K$ and $d = \sum_{k=1}^K d_k$, implementation of Algorithm 1 remains fast in high-dimensional settings, provided that each $d_k \le 10$ with $k = 2, \dots, K$. When there is $k_0 \in \{2, \dots, K\}$ such that $d_{k_0} > 10$, Algorithm 1 can still be used with a substantial requirement of time to obtain all the necessary and sufficient permutations of $\pi_{k_0}$. But,
as the selected permutations are not going to change for a given $d_k$ with $k \ge 2$, a table of selected permutations for different values of $d_k$ will avoid tuning this algorithm all the time.

Regarding the $R_{\min}$ ERs, we can check that $2^K = 2^{\frac{d-d_1}{2}} \le R_{\min} \le 2^{d-d_1}$ because $\binom{d_k}{j_0,k} = \binom{\frac{d_k}{2}}{\frac{d_k}{2}} \le 2^{d_k}$ with $\left[\frac{d_k}{2}\right]$ the largest integer that is less than $\frac{d_k}{2}$. Thus, computing all the effects of inputs will require at least $2^{\frac{d-d_1}{2}}$ ERs and at most $2^{d-d_1}$ ERs, which grow exponentially with respect to $d(1-\alpha)$ with $\alpha := \frac{d_1}{d} \in [0, 1]$. If the effects of $\mathbf{X}_{\pi_{k_0}}$ are not significant, then $\alpha = \frac{d_1+d_{k_0}}{d}$. Based on the values of $\alpha$, different conclusions can be drawn.

## 4. Dependent Multivariate Sensitivity Analysis

This section extends dependent generalized sensitivity indices (dGSIs) for models with NIVVs introduced in [16] by defining the dGSIs for every subset of inputs and for functional outputs. Since the ER given by Equation (18) includes only independent variables, we define dGSIs by relying on the multivariate sensitivity analysis ([2–4,44,45]). To ensure that the proposed dGSIs in this section are well defined, assume that

$$ (A1): 0 < \int_{\Theta} \mathbb{E}[||f(\mathbf{X}, \theta)||_{L^2}^2] d\theta < +\infty. $$

Definitions of GSIs and dGSIs are based on sensitivity functionals (SFs), which contain information about the single and overall contributions of input variables over the whole functional outputs [2,6,16,46]. To define SFs, recall that $0 \le p_k \le d_k$, $\mathbf{v}_k = (w_{2,k},...,w_{p_k,k})$ and $\mathbf{u}_k = (w_{1,k},\mathbf{v}_k)$ with $k = 2,...,K$, and define $\bar{\mathbf{u}}_k := (w_{p_k+1,k},...,w_{d_k,k})$. According to Lemmas 3 and 4, for all $u \subseteq \{1,...,d\}$, there exists $u_0 \subseteq \pi_1$, a vector $\mathbf{w}_k \in \mathcal{P}_k$ and $p_k$ with $k = 2,...,K$ such that $\mathbf{X}_u = (\mathbf{X}_{u_0}, \mathbf{X}_{u_2},...,\mathbf{X}_{u_K})$ and

$$ f(\mathbf{X}, \theta) | \mathbf{X}_u \stackrel{d}{=} g_{\ell}(\mathbf{X}_{\pi_1}, X_{w_{1,2}}, \mathbf{Z}_{w_{\sim 1,2}}, \dots, X_{w_{1,K}}, \mathbf{Z}_{w_{\sim 1,K}}, \theta) | \mathbf{X}_{u_0}, \mathbf{R}_{u_2}, \dots, \mathbf{R}_{u_K}, $$

where $\mathbf{R}_{u_k} := (X_{w_{1,k}}, \mathbf{Z}_{w_k})$, $k = 2,...,K$. Thus, the effects of inputs $\mathbf{X}_u$ are equal to the effects of $\mathbf{X}_{u_0}, \mathbf{R}_{u_2},...,\mathbf{R}_{u_K}$ using $g_{\ell}$. For concise notations, we use

$$ \begin{gathered} \mathbf{X}_{\{\pi_1, s\}} := (X_j, \forall j \in \{\pi_1, s\}); \quad \mathbf{Z}_{w_{\sim 1}} = (\mathbf{Z}_{w_{\sim 1,2}}, \dots, \mathbf{Z}_{w_{\sim 1,K}}), \\ \mathbf{R}_u := (\mathbf{X}_{u_0}, \mathbf{R}_{u_2}, \dots, \mathbf{R}_{u_K}); \quad \mathbf{R}_{\sim u} := (\mathbf{X}_{\pi_1 \setminus u_0}, \mathbf{Z}_{\bar{\mathbf{u}}_2}, \dots, \mathbf{Z}_{\bar{\mathbf{u}}_K}). \end{gathered} $$

Note that $\mathbf{R} := (\mathbf{R}_u, \mathbf{R}_{\sim u})$ is a partition of $(\mathbf{X}_{\{\pi_1,s\}}, \mathbf{Z}_{w_{\sim 1}})$.

The first-order SF of $\mathbf{X}_u$ with $u = \{u_0, \mathbf{u}_2, ..., \mathbf{u}_K\}$ is given by

$$ f_{u}^{fo}(\mathbf{R}_{u}, \theta) := \mathbb{E}[g_{\ell}(\mathbf{X}_{\{\pi_{1}, s\}}, \mathbf{Z}_{w_{\sim 1}}, \theta) | \mathbf{R}_{u}] - \mathbb{E}[g_{\ell}(\mathbf{X}_{\{\pi_{1}, s\}}, \mathbf{Z}_{w_{\sim 1}}, \theta)], \quad (19) $$

and the total SF, which contains the overall information about $\mathbf{X}_u$, is given by

$$ f_{u}^{tot}(\mathbf{X}_{\{\pi_{1}, s\}}, \mathbf{Z}_{w_{\sim 1}}, \theta) := g_{\ell}(\mathbf{X}_{\{\pi_{1}, s\}}, \mathbf{Z}_{w_{\sim 1}}, \theta) - \mathbb{E}_{\mathbf{R}_{u}}[g_{\ell}(\mathbf{X}_{\{\pi_{1}, s\}}, \mathbf{Z}_{w_{\sim 1}}, \theta)], \quad (20) $$

where $\mathbb{E}_{\mathbf{R}_u}$ means that the expectation is taken with respect to $\mathbf{R}_u$. The SFs given by (19) and (20) are random processes, and their components may be correlated and/or dependent. Using the variance–covariance as an importance measure, the definitions of dGSIs rely on the cross-covariances of SFs. For $\theta_1, \theta_2 \in \Theta$, the first-order cross-covariance or the cross-covariance of $f_{u}^{fo}$ is given by

$$ \Sigma_u(\theta_1, \theta_2) := \mathbb{E}[f_{u}^{fo}(\mathbf{R}_u, \theta_1) f_{u}^{fo}(\mathbf{R}_u, \theta_2)^T]. $$
Furthermore, the cross-covariances of $f_u^{tot}$ and $f$ are given as follows:

$$ \Sigma_u^{tot}(\theta_1, \theta_2) := \mathbb{E} \left[ f_u^{tot}(\mathbf{X}_{\{\pi_1, s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta_1) f_u^{tot}(\mathbf{X}_{\{\pi_1, s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta_2)^T \right], $$

$$ \Sigma(\theta_1, \theta_2) := \mathbb{E} \left[ g_{\ell}(\mathbf{X}_{\{\pi_1, s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta_1) g_{\ell}(\mathbf{X}_{\{\pi_1, s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta_2)^T \right] \\ - \mathbb{E} \left[ g_{\ell}(\mathbf{X}_{\{\pi_1, s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta_1) \right] \mathbb{E} \left[ g_{\ell}(\mathbf{X}_{\{\pi_1, s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta_2)^T \right]. $$

To account for the different properties of the aforementioned SFs, we distinguish three different types of dGSIs.

**Definition 3.** Consider the cross-covariances of SFs, and assume that (A1) holds.

(i) For the first-type dGSIs, the first-order and total indices are given by

$$ dGSI_{u,1} := \frac{\int_{\Theta} \mathrm{Tr}(\Sigma_u(\theta, \theta)) d\theta}{\int_{\Theta} \mathrm{Tr}(\Sigma(\theta, \theta)) d\theta}, \qquad dGSI_{T_{u,1}} := \frac{\int_{\Theta} \mathrm{Tr}(\Sigma_u^{tot}(\theta, \theta)) d\theta}{\int_{\Theta} \mathrm{Tr}(\Sigma(\theta, \theta)) d\theta}. \quad (21) $$

(ii) The second-type dGSIs are defined as follows:

$$ dGSI_{u,2} := \left( \frac{\int_{\Theta} ||\Sigma_u(\theta, \theta)||_F^2 d\theta}{\int_{\Theta} ||\Sigma(\theta, \theta)||_F^2 d\theta} \right)^{1/2}, \qquad dGSI_{T_{u,2}} := \left( \frac{\int_{\Theta} ||\Sigma_u^{tot}(\theta, \theta)||_F^2 d\theta}{\int_{\Theta} ||\Sigma(\theta, \theta)||_F^2 d\theta} \right)^{1/2}. \quad (22) $$

(iii) The third-type dGSIs are given as follows:

$$ dGSI_{u,3} := \left( \frac{\int_{\Theta^2} ||\Sigma_u(\theta_1, \theta_2)||_F^2 d\theta_1 d\theta_2}{\int_{\Theta^2} ||\Sigma(\theta_1, \theta_2)||_F^2 d\theta_1 d\theta_2} \right)^{1/2}, \qquad dGSI_{T_{u,3}} := \left( \frac{\int_{\Theta^2} ||\Sigma_u^{tot}(\theta_1, \theta_2)||_F^2 d\theta_1 d\theta_2}{\int_{\Theta^2} ||\Sigma(\theta_1, \theta_2)||_F^2 d\theta_1 d\theta_2} \right)^{1/2}. \quad (23) $$

The first-type and the second-type dGSIs treat independently the outputs $f(\mathbf{X}, \theta_1)$ and $f(\mathbf{X}, \theta_2)$ when $\theta_1 \neq \theta_2$, but the second-type dGSIs account for the correlations among the components of SFs from the same output $f(\mathbf{X}, \theta)$. Furthermore, the third-type dGSIs account for the correlations among the cross-components of SFs.

### 4.1. Properties of Dependent Generalized Sensitivity Indices

The two types of dGSIs share the same properties as those proposed in [16] for $0 \le p_k \le 1$, only. Proposition 2 extends such properties for all $p_k \in \{0, \dots, d_k\}$.

**Proposition 2.** Under (A1), consider the dGSIs of Definition 3. Then,

$$ 0 \le dGSI_{u,1} \le dGSI_{T_{u,1}} \le 1; \qquad 0 \le dGSI_{u,2} \le dGSI_{T_{u,2}} \le 1. $$

Moreover, if the cross-covariances are positive semi-definite, then we have

$$ 0 \le dGSI_{u,3} \le dGSI_{T_{u,3}} \le 1. $$

**Proof.** See Appendix F. □

When the total dGSI of $\mathbf{X}_u = (\mathbf{X}_{u0}, \mathbf{X}_{u1}, \dots, \mathbf{X}_{uK})$ is zero or almost zero (i.e., $dGSI_{T_{u,\bullet}} \approx 0$), we have to fix $\mathbf{X}_u$ using the DM associated with each $\mathbf{X}_{\pi_k}$, $k = 2, \dots, K$. Indeed, for the DM $\mathbf{X}_{u_k} = r_{\sim u_k}(\mathbf{X}_{\sim u_k}, \mathbf{Z}_{u_k})$, fixing $\mathbf{X}_{u_K}$ comes down to fix $\mathbf{Z}_{u_k}$ to its nominal values. Since we can compute the total dGSI of each block of NIVs, that is, $\mathbf{X}_{\pi_k}$ using any ER of $f$, it
becomes possible to quickly identify the non-influential block of NIVs, and then put our computational efforts on the most important groups of NIVs.

**Remark 1.** Given $j_1, j_2 \in \{1, \dots, d\}$, the same rankings of the inputs $X_{j_1}$ and $X_{j_2}$ using either $dGSI_{T_{j_1},1}$ or $dGSI_{T_{j_1},2}$ with $j \in \{j_1, j_2\}$ are obtained under the assumptions $\Sigma_{j_1}^{\text{tot}}(\theta, \theta) \le \Sigma_{j_2}^{\text{tot}}(\theta, \theta)$ or $\Sigma_{j_2}^{\text{tot}}(\theta, \theta) \le \Sigma_{j_1}^{\text{tot}}(\theta, \theta)$. Note that $A_1 \preceq A_2$ means that $A_2 - A_1$ is a positive semi-definite matrix, a.k.a. the Loewner partial ordering between matrices (see also Section 6.3 in [16]).

**Remark 2.** Case of the multivariate dynamic function

Consider a model that is evaluated at $\mathbf{X}$ and provides $n$ dynamic(s), such as a spatiotemporal model. Such a model is a particular case of multivariate and functional outputs. Indeed, the multivariate dynamic model given by $f: \mathbb{R}^d \times [0, T] \to \mathbb{R}^n$ and $f(\mathbf{X}, t) \in \mathbb{R}^n$ with $t \in [0, T]$ is mathematically identical to $f(\mathbf{X}, \theta)$ using $\Theta = [0, T]$ and $\theta = t$.

## 4.2. Case of the Multivariate Response Models

When $\Theta = \{\theta_0\}$, the multivariate and functional outputs come down to $f(\mathbf{X}, \theta_0) =: h(\mathbf{X})$ with $h: \mathbb{R}^d \to \mathbb{R}^n$. Thus, the dGSIs of Definition 3 can be adapted for quantifying the effect of inputs. It is worth noting that the third-type dGSIs are equal to the second-type dGSIs, and both types of dGSIs boil down to the second-type dGSIs provided in Definition 4. Moreover, the cross-covariances become the covariances, and we use $\Sigma_u := \Sigma_u(\theta_0, \theta_0)$, $\Sigma_{u1}^{\text{tot}} := \Sigma_{u1}^{\text{tot}}(\theta_0, \theta_0)$ and $\Sigma := \Sigma(\theta_0, \theta_0)$.

**Definition 4.** Consider the above covariances of SFs and assume (A1) holds.

(i) The first-type dGSIs for a given multivariate response function are

$$dGSI_{u,1}^{M} := \frac{\mathrm{Tr}(\Sigma_{u})}{\mathrm{Tr}(\Sigma)}, \qquad dGSI_{T_{u},1}^{M} := \frac{\mathrm{Tr}(\Sigma_{u}^{\mathrm{tot}})}{\mathrm{Tr}(\Sigma)}. \qquad (24)$$

(ii) The second-type dGSIs are defined as follows:

$$dGSI_{u,2}^{M} := \frac{\|\Sigma_{u}\|_{F}}{\|\Sigma\|_{F}}, \qquad dGSI_{T_{u},2}^{M} := \frac{\|\Sigma_{u}^{\text{tot}}\|_{F}}{\|\Sigma\|_{F}}. \qquad (25)$$

Also, in the case where $n = 1$, the two types of dGSIs of Definition 4 are equal and boil down to dependent sensitivity indices (dSIs) for single-response models (see Definition 5).

**Definition 5.** For a function $f: \mathbb{R}^d \to \mathbb{R}$ ($n = 1$), the dSIs of $\mathbf{X}_u$ are given by

$$dS_u := \frac{\Sigma_u}{\Sigma}, \qquad dS_{T_u} := \frac{\Sigma_{u}^{\text{tot}}}{\Sigma}. \qquad (26)$$

# 5. Estimators of Dependent Generalized Sensitivity Indices

In this section, we provide unbiased estimators of the cross-covariances and covariances of SFs, consistent estimators of dGSIs, and their asymptotic distributions. For the sake of simplicity, we provide estimators of covariances and dGSIs using the functions $g_\ell$ that include only independent variables, that is, $\mathbf{X}_{\{\pi_1,s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}$ Note that the Supplementary Materials provides such estimators using $f$ directly thanks to the relation

$$g_{\ell}(\mathbf{X}_{\{\pi_1,s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta) = f(\mathbf{X}_{\{\pi_1,s\}}, r_{\mathbf{s}}(\mathbf{X}_{\mathbf{s}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}}), \theta) \\ = f(\mathbf{X}_{\{\pi_1,s\}}, r_{w_{1,2}}(X_{w_{1,2}}, \mathbf{Z}_{\mathbf{w}_{\sim 1,2}}), \dots, r_{w_{1,K}}(X_{w_{1,K}}, \mathbf{Z}_{\mathbf{w}_{\sim 1,K}}), \theta).$$
To derive the estimators of cross-covariances that are useful for computing dGSIs, we are given two i.i.d. samples from $(\mathbf{X}_{\{\pi_1,s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}})$, that is, $\{(\mathbf{X}_{i,\{\pi_1,s\}}', \mathbf{Z}_{i,\mathbf{w}_{\sim 1}}^{(1)})\}_{i=1}^m$ and $\{(\mathbf{X}_{i,\{\pi_1,s\}}^{(2)}, \mathbf{Z}_{i,\mathbf{w}_{\sim 1}}^{(2)})\}_{i=1}^m$, and we use $\mathbf{X}_i^{(k)} := (\mathbf{X}_{i,\{\pi_1,s\}}^{(k)}, r_s(\mathbf{X}_{i,s}^{(k)}, \mathbf{Z}_{i,\mathbf{w}_{\sim 1}}^{(k)}))$, $k = 1, 2$, $i = 1, \dots, m$. Since $\mathbf{R} := (\mathbf{R}_u, \mathbf{R}_{\sim u})$ is a partition of $(\mathbf{X}_{\{\pi_1,s\}}, \mathbf{Z}_{\mathbf{w}_{\sim 1}})$, we can deduce the following two i.i.d. samples: $\{\mathbf{R}_i^{(1)} := (\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(1)})\}_{i=1}^m$ and $\{\mathbf{R}_i^{(2)} := (\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(2)})\}_{i=1}^m$.

**Theorem 2.** Assume (A1) holds. Then, unbiased and consistent estimators of $\Sigma_u(\theta_1, \theta_2)$, $\Sigma_u^{tot}(\theta_1, \theta_2)$ and $\Sigma(\theta_1, \theta_2)$ are, respectively, given by

$$ \Sigma_u(\theta_1, \theta_2) := \frac{1}{2m} \sum_{i=1}^{m} [g_{\ell}(\mathbf{R}_i^{(1)}, \theta_1) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)}, \theta_1)] [g_{\ell}(\mathbf{R}_i^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}, \theta_2) - g_{\ell}(\mathbf{R}_i^{(2)}, \theta_2)]^T; $$

$$ \Sigma_{iu}^{tot}(\theta_1, \theta_2) := \frac{1}{2m} \sum_{i=1}^{m} [g_{\ell}(\mathbf{R}_i^{(1)}, \theta_1) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)}, \theta_1)] [g_{\ell}(\mathbf{R}_i^{(1)}, \mathbf{R}_{i,\sim u}^{(1)}, \theta_2) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)}, \theta_2)]^T; $$

$$ \Sigma(\theta_1, \theta_2) := \frac{1}{2m} \sum_{i=1}^{m} [g_{\ell}(\mathbf{R}_i^{(1)}, \theta_1) - g_{\ell}(\mathbf{R}_i^{(2)}, \theta_1)] [g_{\ell}(\mathbf{R}_i^{(1)}, \theta_2) - g_{\ell}(\mathbf{R}_i^{(2)}, \theta_2)]^T \\ = \frac{1}{2m} \sum_{i=1}^{m} [f(\mathbf{X}_i^{(1)}, \theta_1) - f(\mathbf{X}_i^{(2)}, \theta_1)] [f(\mathbf{X}_i^{(1)}, \theta_2) - f(\mathbf{X}_i^{(2)}, \theta_2)]^T. $$

**Proof.** See Appendix G. □

Note that when $\theta_1 = \theta_2 = \theta$, the estimators provided in Theorem 2 are minimum-variance unbiased estimators (MVUEs). Using Theorem 2, we deduce the estimators of the variance-covariances of SFs for the multivariate response models and single-response models. To provide such results for the vector-valued function $f: \mathbb{R}^d \to \mathbb{R}^m$ in Corollary 1, let us consider the following symmetric kernels:

$$ K_u^{fo}(\mathbf{R}_i^{(1)}, \mathbf{R}_i^{(2)}) := [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})] [g_{\ell}(\mathbf{R}_i^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_i^{(2)})]^T \\ + [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_i^{(2)})] [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^T; $$

$$ K_u^{tot}(\mathbf{R}_i^{(1)}, \mathbf{R}_i^{(2)}) := [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})] [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^T \\ + [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_i^{(2)})] [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_i^{(2)})]^T; $$

$$ K(\mathbf{R}_i^{(1)}, \mathbf{R}_i^{(2)}) := [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})] [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^T \\ + [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_i^{(2)})] [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_i^{(2)})]^T. $$
**Corollary 1.** Assume that f has finite fourth moments (i.e., (A2)) and (A1) hold. Then, the MVUEs and consistent estimators of the covariance matrices $\Sigma_u, \Sigma_u^{tot}$, and $\Sigma$ are, respectively, given by

$$ \hat{\Sigma}_u := \frac{1}{4m} \sum_{i=1}^{m} K_u^{fo}(\mathbf{R}_i^{(1)}, \mathbf{R}_i^{(2)}), \qquad \hat{\Sigma}_u^{tot} := \frac{1}{4m} \sum_{i=1}^{m} K_u^{tot}(\mathbf{R}_i^{(1)}, \mathbf{R}_i^{(2)}); \\ \hat{\Sigma} := \frac{1}{4m} \sum_{i=1}^{m} K(\mathbf{R}_i^{(1)}, \mathbf{R}_i^{(2)}). $$

**Proof.** See Appendix H. $\square$

For single-response models (i.e., $n = 1$), the MVUEs of the covariances of SFs in Corollary 1 have simple expressions, given below.

**Corollary 2.** Under the assumptions (A1)-(A2) and $n = 1$, the MVUEs $\hat{\Sigma}_u, \hat{\Sigma}_u^{tot}$, and $\hat{\Sigma}$ become, respectively,

$$ \hat{\sigma}_u := \frac{1}{2m} \sum_{i=1}^{m} [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})] [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_i^{(2)})]; \\ \hat{\sigma}_u^{tot} := \frac{1}{4m} \sum_{i=1}^{m} \{[g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^2 + [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_i^{(2)})]^2\}; \\ \hat{\sigma} := \frac{1}{4m} \sum_{i=1}^{m} \{[g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_i^{(2)})]^2 + [g_{\ell}(\mathbf{R}_{i,u}^{(1)}, \mathbf{R}_{i,\sim u}^{(2)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^2\}. $$

If we are only interested in the total effects, we should use the expressions of $\hat{\Sigma}_u^{tot}$ and $\hat{\sigma}_u^{tot}$ given, respectively, by

$$ \frac{1}{2m} \sum_{i=1}^{m} [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})] [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^T; \\ \frac{1}{2m} \sum_{i=1}^{m} [g_{\ell}(\mathbf{R}_i^{(1)}) - g_{\ell}(\mathbf{R}_{i,u}^{(2)}, \mathbf{R}_{i,\sim u}^{(1)})]^2. $$

Using the results from Theorem 2, Corollary 1, and Corollary 2, we derive the estimators of dGSIs and dSIs in Corollary 3, Theorem 3, and Corollary 4, respectively.

**Corollary 3.** Assume that (A1) and (A2) hold. If we observe the model outputs at $\theta_{\ell} \in \Theta$ with $\ell = 1, \dots, L$, then

(i) the consistent estimators of the first-type dGSIs are given as follows:

$$ \widehat{dGSI}_{u,1} := \frac{\sum_{\ell=1}^{L} \mathrm{Tr}(\widehat{\Sigma}_u(\widehat{\theta}_{\ell}, \widehat{\theta}_{\ell}))}{\sum_{\ell=1}^{L} \mathrm{Tr}(\widehat{\Sigma}(\widehat{\theta}_{\ell}, \widehat{\theta}_{\ell}))} \xrightarrow{P} dGSI_{u,1}, $$

when $L \to \infty$, $m \to \infty$, where $\xrightarrow{P}$ denotes the convergence in probability.

$$ \widehat{dGSI}_{T_u,1} := \frac{\sum_{\ell=1}^{L} \mathrm{Tr}(\widehat{\Sigma}_u^{tot}(\widehat{\theta}_{\ell}, \widehat{\theta}_{\ell}))}{\sum_{\ell=1}^{L} \mathrm{Tr}(\widehat{\Sigma}(\widehat{\theta}_{\ell}, \widehat{\theta}_{\ell}))} \xrightarrow{P} dGSI_{T_u,1}. $$
(ii) The estimators of the second-type dGSIs are given as follows:

$$ \widehat{dGSI}_{u,2} := \left( \frac{\sum_{\ell=1}^{L} \left| \left| \Sigma_u(\widehat{\theta}_\ell, \theta_\ell) \right|_F \right|^2}{\sum_{\ell=1}^{L} \left| \left| \Sigma(\widehat{\theta}_\ell, \theta_\ell) \right|_F \right|^2} \right)^{1/2} \xrightarrow{P} dGSI_{u,2}, $$

$$ d\widehat{GSI}_{T_u,2} := \left( \frac{\sum_{\ell=1}^{L} \left| \left| \Sigma_u^{tot}(\widehat{\theta}_\ell, \theta_\ell) \right|_F \right|^2}{\sum_{\ell=1}^{L} \left| \left| \Sigma(\widehat{\theta}_\ell, \theta_\ell) \right|_F \right|^2} \right)^{1/2} \xrightarrow{P} dGSI_{T_u,2}, $$

(iii) The estimators of the third-type dGSIs are given as follows:

$$ \widehat{dGSI}_{u,3} := \left( \frac{\sum_{\ell_1=1}^{L} \sum_{\ell_2=1}^{L} \left| \left| \Sigma_u(\widehat{\theta}_{\ell_1}, \theta_{\ell_2}) \right|_F \right|^2}{\sum_{\ell_1=1}^{L} \sum_{\ell_2=1}^{L} \left| \left| \Sigma(\widehat{\theta}_{\ell_1}, \theta_{\ell_2}) \right|_F \right|^2} \right)^{1/2} \xrightarrow{P} dGSI_{u,3}, $$

$$ d\widehat{GSI}_{T_u,3} := \left( \frac{\sum_{\ell_1=1}^{L} \sum_{\ell_2=1}^{L} \left| \left| \Sigma_u^{tot}(\widehat{\theta}_{\ell_1}, \theta_{\ell_2}) \right|_F \right|^2}{\sum_{\ell_1=1}^{L} \sum_{\ell_2=1}^{L} \left| \left| \Sigma(\widehat{\theta}_{\ell_1}, \theta_{\ell_2}) \right|_F \right|^2} \right)^{1/2} \xrightarrow{P} dGSI_{T_u,3}. $$

**Proof.** Using Theorem 2, the results hold by applying the Slutsky theorem. □

For computing the first-order and total-effect covariances of each $X_j$, $2m(d+1)$ model runs are needed. Some of such runs can be combined for computing the model outputs' covariance. From now on, $M$ model runs are used for computing the model outputs' covariance. The operator $\text{Vec}(\cdot)$ transforms a matrix $\Sigma \in \mathbb{R}^{n \times n}$ into a vector, that is, $\text{Vec}(\Sigma) \in \mathbb{R}^{n^2}$ and $O \in \mathbb{R}^{n \times n}$ denotes the null matrix.

**Theorem 3.** Assume that (A1) and (A2) hold, $m \to +\infty$, $M \to +\infty$ and $m/M \to 0$.

(i) The estimators of the first-type dGSIs are given as follows:

$$ \widehat{dGSI}_{u,1}^{M} := \frac{\text{Tr}(\widehat{\Sigma}_{u})}{\text{Tr}(\widehat{\Sigma})} \xrightarrow{P} dGSI_{u,1}^{M}; \qquad \widehat{dGSI}_{T_u,1}^{M} := \frac{\text{Tr}(\widehat{\Sigma}_{u}^{\text{tot}})}{\text{Tr}(\widehat{\Sigma})} \xrightarrow{P} dGSI_{T_u,1}^{M}, \quad (27) $$

with the following asymptotic distributions:

$$ \sqrt{m} \left( d\widehat{GSI}_{u,1}^{M} - dGSI_{u,1}^{M} \right) \xrightarrow{D} \mathcal{N} \left( 0, \frac{\mathbb{V}[\text{Tr}\{K_u^{fo}(\mathbf{R}_1^{(1)}, \mathbf{R}_1^{(2)})\}]}{(\text{Tr}(\Sigma))^2} \right); $$

$$ \sqrt{m} \left( d\widehat{GSI}_{T_u,1}^{M} - dGSI_{T_u,1}^{M} \right) \xrightarrow{D} \mathcal{N} \left( 0, \frac{\mathbb{V}[\text{Tr}\{K_u^{tot}(\mathbf{R}_1^{(1)}, \mathbf{R}_1^{(2)})\}]}{(\text{Tr}(\Sigma))^2} \right). $$

(ii) For the second-type dGSIs, we have

$$ d\widehat{GSI}_{u,2}^{M} := \frac{\left| \left| \widehat{\Sigma}_{u} \right|_{F}}{\left| \left| \widehat{\Sigma} \right|_{F}} \xrightarrow{P} dGSI_{u,2}^{M}; \qquad \widehat{dGSI}_{T_u,2}^{M} := \frac{\left| \left| \widehat{\Sigma}_{u}^{\text{tot}} \right|_{F}}{\left| \left| \widehat{\Sigma} \right|_{F}} \xrightarrow{P} dGSI_{T_u,2}^{M}, \quad (28) $$

$$ \sqrt{m} \left( d\widehat{GSI}_{u,2}^{M} - dGSI_{u,2}^{M} \right) \xrightarrow{D} \mathcal{N} \left( 0, \frac{\text{Vec}(\Sigma_u)^T \mathbb{V} [ \text{Vec} \{ K_u^{fo}(\mathbf{R}_1^{(1)}, \mathbf{R}_1^{(2)}) \} ] \text{Vec}(\Sigma_u)}{\|\Sigma_u\|_F^2 \|\Sigma\|_F^2} \right). $$
$$
\sqrt{m} \left( d\widehat{GSI}_{T_u, 2}^{M} - dGSI_{T_u, 2}^{M} \right) \xrightarrow{D} \mathcal{N} \left( 0, \frac{\mathrm{Vec}(\Sigma_u^{\mathrm{tot}})^T \mathbb{V} \left[ \mathrm{Vec}\left\{ K_u^{\mathrm{tot}} \left( \mathbf{R}_1^{(1)}, \mathbf{R}_1^{(2)} \right) \right\} \right] \mathrm{Vec}(\Sigma_u^{\mathrm{tot}})}{||\Sigma_u^{\mathrm{tot}}||_F^2 ||\Sigma||_F^2} \right),
$$
provided that $\Sigma_u \neq 0$ and $\Sigma_u^{\mathrm{tot}} \neq 0$.

**Proof.** See Appendix I. □

Using Theorem 3, we give the estimators of dSIs (see Definition 5) for real-valued functions in Corollary 4.

**Corollary 4.** Let $n = 1$ and $\sigma^2 := \mathbb{V}[f(\mathbf{X})]$. Assume that (A1)-(A2) hold, $m \to +\infty$, $M \to +\infty$ and $m/M \to 0$. The estimators of dSIs are given as follows:

$$
\widehat{dS}_u := \frac{\hat{\sigma}_u}{\hat{\sigma}} \xrightarrow{P} dS_u; \quad \sqrt{m}(\widehat{dS}_u - dS_u) \xrightarrow{D} \mathcal{N} \left( 0, \frac{\mathbb{V}[K_u^{fo}(\mathbf{R}_1^{(1)}, \mathbf{R}_1^{(2)})]}{\sigma^4} \right); \quad (29)
$$

$$
\widehat{dS}_{T_u} := \frac{\widehat{\sigma}_{u}^{\text{tot}}}{\widehat{\sigma}} \xrightarrow{P} dS_{T_u}; \quad \sqrt{m}(\widehat{dS}_{T_u} - dS_{T_u}) \xrightarrow{D} \mathcal{N} \left( 0, \frac{\mathbb{V}[K_{u}^{\text{tot}}(\mathbf{R}_{1}^{(1)}, \mathbf{R}_{1}^{(2)})]}{\sigma^4} \right). \quad (30)
$$

The computation of the dGSIs or dSIs of $X_u$ for all $u \subseteq \{1, \dots, d\}$ using the above estimators will require $R_{\min}$ ERs of $f$. When we are only interested in $u \subseteq \{1, \dots, d\}$ with $|u| = 1$, $R_{\min,1} := \max(d_2, \dots, d_K)$, the ERs of $f$ are necessary and sufficient, and we need $2 \times m(d+1)$ model evaluations to estimate the first-order and total dGSIs or dSIs of $X_j$ for all $j \in \{1, \dots, d\}$. Note that with the same $R_{\min,1}$ ERs and additional model runs, we are able to compute the effects of other subsets of inputs and interactions.

**6. Analytical and Numerical Results**

In this section, we illustrate our approach by means of analytical test cases, including
a dynamic model so as to highlight some theoretical properties of the new indices.

6.1. Linear Function Without Explicit Interaction ($d = 3, n = 1$)

$$
\text{We consider } f(\mathbf{X}) = X_1 + X_2 + X_3 \text{ with } \mathbf{X} \sim \mathcal{N} \left( \mathbf{0}, \left[ \begin{array}{cccc} \sigma_1^2 & \rho_{12}\sigma_1\sigma_2 & \rho_{13}\sigma_1\sigma_3 & \\ \rho_{12}\sigma_1\sigma_2 & \sigma_2^2 & \rho_{23}\sigma_2\sigma_3 & \\ \rho_{13}\sigma_1\sigma_3 & \rho_{23}\sigma_2\sigma_3 & \sigma_3^2 & \end{array} \right] \right).
$$

A dependency function of $\mathbf{X}$ is given by $(X_2, X_3) = r_1(X_1, Z_2, Z_3)$, where

$$
\left\{
\begin{aligned}
X_2 &= \frac{\rho_{12}\sigma_2}{\sigma_1} X_1 + \sqrt{1 - \rho_{12}^2} Z_2 \\
X_3 &= \frac{\rho_{13}\sigma_3}{\sigma_1} X_1 + \frac{\sigma_3(\rho_{23} - \rho_{12}\rho_{13})}{\sigma_2\sqrt{1-\rho_{12}^2}} Z_2 + \sqrt{\frac{1-\rho_{12}^2-\rho_{13}^2-\rho_{23}^2+2\rho_{12}\rho_{13}\rho_{23}}{1-\rho_{12}^2}} Z_3,
\end{aligned}
\right.
$$

$Z_i \sim \mathcal{N}(0, \sigma_i^2), i = 2, 3,$ and an ER of $f$ is given by

$$
g_1(X_1, Z_2, Z_3) = \left(1 + \frac{\rho_{12}\sigma_2}{\sigma_1} + \frac{\rho_{13}\sigma_3}{\sigma_1}\right) X_1 + \left(\sqrt{1 - \rho_{12}^2} + \frac{\sigma_3(\rho_{23} - \rho_{12}\rho_{13})}{\sigma_2\sqrt{1 - \rho_{12}^2}}\right) Z_2 \\
+ \sqrt{\frac{1 - \rho_{12}^2 - \rho_{13}^2 - \rho_{23}^2 + 2\rho_{12}\rho_{13}\rho_{23}}{1 - \rho_{12}^2}} Z_3.
$$

Using such ERs of $f$, we have the following dSIs of $X_1$ and $(X_1, X_2)$:
$$dS_1 = dS_{T_1} = \frac{(\sigma_1 + \rho_{12}\sigma_2 + \rho_{13}\sigma_3)^2}{\sum_{j=1}^{3} \sigma_j^2 + 2\rho_{12}\sigma_1\sigma_2 + 2\rho_{13}\sigma_1\sigma_3 + 2\rho_{23}\sigma_2\sigma_3},$$

$$dS_{12} = dS_{T_{12}} = \frac{(1 - \rho_{12}^2)(\sigma_1 + \rho_{12}\sigma_2 + \rho_{13}\sigma_3)^2 + (\sigma_2(1 - \rho_{12}^2) + \sigma_3(\rho_{23} - \rho_{12}\rho_{13}))^2}{(1 - \rho_{12}^2)(\sum_{j=1}^{3} \sigma_j^2 + 2\rho_{12}\sigma_1\sigma_2 + 2\rho_{13}\sigma_1\sigma_3 + 2\rho_{23}\sigma_2\sigma_3)},$$

respectively. Using the same reasoning and knowing that $R_{min} = 3$, the two extra ERs of f lead to the remaining results. When using $g_2(X_2, Z_3, Z_1)$, we have

$$dS_2 = dS_{T_2} = \frac{(\sigma_2 + \rho_{12}\sigma_1 + \rho_{23}\sigma_3)^2}{\sum_{j=1}^{3} \sigma_j^2 + 2\rho_{12}\sigma_1\sigma_2 + 2\rho_{13}\sigma_1\sigma_3 + 2\rho_{23}\sigma_2\sigma_3},$$

$$dS_{23} = dS_{T_{23}} = \frac{(1 - \rho_{23}^2)(\sigma_2 + \rho_{12}\sigma_1 + \rho_{23}\sigma_3)^2 + (\sigma_3(1 - \rho_{23}^2) + \sigma_1(\rho_{13} - \rho_{12}\rho_{23}))^2}{(1 - \rho_{23}^2)(\sum_{j=1}^{3} \sigma_j^2 + 2\rho_{12}\sigma_1\sigma_2 + 2\rho_{13}\sigma_1\sigma_3 + 2\rho_{23}\sigma_2\sigma_3)}.$$

Likewise, using $g_3(X_3, Z_1, Z_2)$, we obtain

$$dS_3 = dS_{T_3} = \frac{(\sigma_3 + \rho_{13}\sigma_1 + \rho_{23}\sigma_2)^2}{\sum_{j=1}^{3} \sigma_j^2 + 2\rho_{12}\sigma_1\sigma_2 + 2\rho_{13}\sigma_1\sigma_3 + 2\rho_{23}\sigma_2\sigma_3},$$

$$dS_{13} = dS_{T_{13}} = \frac{(1 - \rho_{13}^2)(\sigma_3 + \rho_{13}\sigma_1 + \rho_{23}\sigma_2)^2 + (\sigma_1(1 - \rho_{13}^2) + \sigma_2(\rho_{12} - \rho_{13}\rho_{23}))^2}{(1 - \rho_{13}^2)(\sum_{j=1}^{3} \sigma_j^2 + 2\rho_{12}\sigma_1\sigma_2 + 2\rho_{13}\sigma_1\sigma_3 + 2\rho_{23}\sigma_2\sigma_3)}.$$

## 6.2. Functional Outputs: Dynamic Model ($d = 2, n = 1$)

The following dynamic model includes two inputs $X \sim N(0, \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix})$:

$$f(\mathbf{X}, t) = (X_1 + X_2 + aX_1X_2)(2 - \alpha(t)) + (X_1^2 + \sqrt{2}X_2)(-1 + \alpha(t)),$$

where $a \in \mathbb{R}$, $t \in [0, 365]$ and $\alpha(t) \in \mathbb{R}$. When $a = 0$, there is no explicit interaction between $X_1$ and $X_2$. To illustrate the difference between both types of dGSIs, we suppose that we have observed the model outputs at $t \in \{t_1, t_2\}$ with $\alpha(t_1) = 1$ and $\alpha(t_2) = 2$. Using the dependency functions $X_2 = r_1(X_1, Z_2) = \rho X_1 + \sqrt{1 - \rho^2}Z_2$ and $X_1 = r_2(X_2, Z_1) = \rho X_2 + \sqrt{1 - \rho^2}Z_1$ with $Z_i \sim N(0, 1)$ and $i = 1, 2$, Figure 1 shows the two types of dGSIs.

The first figure (top-left panel) depicts the first-type and second-type dGSIs, and we see that both types of dGSIs give the same ranking of inputs for negative values of the correlation. In the absence of correlation (i.e., $\rho = 0$), $X_1$ and $X_2$ have the same effect according to the first-type dGSIs, while the second-type dGSIs identify $X_2$ as the most influential input. When $\rho > 0.5$, the first-type dGSIs show that $X_2$ is the most important input, while the second-type dGSIs suggest that both inputs have the same total effects. The second figure (top-right panel) compares the prime second-type and the second-type total dGSIs, and the results are similar to those of the first figure. The third figure (bottom-left panel) compares the first-order dGSIs of the first and second types, and it comes out that such dGSIs give different main effects of inputs for positive correlations. In the last figure, the first-type and second-type dGSIs give the same ranking of inputs except for $\rho = 0$. As different rankings of inputs can happen using both types of dGSIs, and knowing that the second-type dGSIs account for the correlations among SFs, we should prefer such indices.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** First-type, prime second-type, and second-type dGSIs for different values of the correlation between the two inputs and for $a = -2, 0$.</figcaption></figure>

### 6.3. Sobol's g-Function ($d = 10, K = 3$)

Here, we consider $f(\mathbf{x}) := \prod_{j=1}^{d=10} \frac{|4x_j - 2| + a_j}{1+a_j}$ with $\mathbf{a} := (35, 35, \dots, 35) \in \mathbb{R}^d$. The $d = 10$ inputs are organized into $K = 3$ blocks as follows:

* $(X_j \sim U(0, 1), j = 4, \dots, 8)$ are independent variables, that is, $\pi_1 = \{4, \dots, 8\}$;

* $\pi_2 = \{1, \dots, 3\}$ and $(X_1, X_2, X_3)$ have a Gaussian copula with $\rho_{12} = 0$, $\rho_{13} = 0.01, \rho_{23} = 0.85$ as the correlation values, and $X_j \sim U(0, 1), j = 1, 2, 3$;

* $\pi_3 = \{9, 10\}$, where $X_9 \sim U(0, 1), X_{10} \sim U(0, 1)$ with $X_9 + X_{10} \le 1$.

The dependency functions of $X_{\pi_3}$ are given by $X_{10}^c = U_{10}(1 - X_9^c)$ and $X_9^c = U_9(1 - X_{10}^c)$ where $U_{10} \sim U(0, 1), U_9 \sim U(0, 1)$; and $U_9, U_{10}$ are independent of $X_9^c \sim \text{Beta}(1, 2), X_{10}^c \sim \text{Beta}(1, 2)$ [20].

Based on the $R_{\min} = 6$ ERs of $f(\mathbf{x})$ (see Appendix J), we computed the dSIs using Sobol's sequences and the sample size $m = 10,000$ (see Table 1). We can see that $X_2, X_3, X_9, X_{10}$ are the most influential inputs. For fixing $X_1$, we need the copula-based dependency model of $X_{\pi_2}$, that is, $X_1 = r_1(X_2, Z_3, Z_1)$. Since the block of inputs $X_{\pi_1}$ is not important, we have also computed the dSIs of the pairs of variables selected out of $(X_{\pi_2}, X_{\pi_3})$ (see Table 1). As expected, the total indices are always greater than the first-order indices.

<table><caption>**Table 1.** Estimates of dSIs for Sobol's g-function with NIVs.</caption><thead><tr><th rowspan="2"></th><th colspan="2">dSI</th><th rowspan="2">dSI</th><th rowspan="2">Total</th></tr><tr><th>Main</th><th>Total</th></tr></thead><tbody><tr><td>X1</td><td>0.089</td><td>0.090</td><td>X1:X2</td><td>0.321</td><td>0.372</td></tr><tr><td>X2</td><td>0.241</td><td>0.294</td><td>X1:X3</td><td>0.321</td><td>0.372</td></tr><tr><td>X3</td><td>0.231</td><td>0.282</td><td>X1:X9</td><td>0.192</td><td>0.216</td></tr><tr><td>X4</td><td>0.093</td><td>0.093</td><td>X1:X10</td><td>0.193</td><td>0.217</td></tr><tr><td>X5</td><td>0.093</td><td>0.093</td><td>X2:X3</td><td>0.298</td><td>0.299</td></tr><tr><td>X6</td><td>0.093</td><td>0.093</td><td>X2:X9</td><td>0.349</td><td>0.426</td></tr><tr><td>X7</td><td>0.093</td><td>0.093</td><td>X2:X10</td><td>0.348</td><td>0.425</td></tr><tr><td>X8</td><td>0.093</td><td>0.093</td><td>X3:X9</td><td>0.334</td><td>0.407</td></tr><tr><td>X9</td><td>0.107</td><td>0.131</td><td>X3:X10</td><td>0.335</td><td>0.409</td></tr><tr><td>X10</td><td>0.108</td><td>0.132</td><td>X9:X10</td><td>0.185</td><td>0.185</td></tr></tbody></table>
# 7. Conclusions

In this paper, we have proposed a new way of deriving the distributions of the model outputs conditional on every subset of inputs using dependency models of random vectors of NIVs. We have provided additional generic dependency models, including empirical or computational dependency functions, of *d*-dimensional random vectors following many distributions, such as copula-based distributions with discrete variables, as well as distributions of inputs of complex mathematical models subjected to constraints. It came out that $\prod_{k=2}^{K} \binom{d_k}{j_0 k}$ different equivalent representations of the model output of interest are necessary and sufficient for recovering all the conditional distributions and for assessing the effects of every subset of NIVs, including their interactions. An algorithm is then provided for selecting such representations or equivalently such dependency models.

Based on such conditional distributions and using variance-covariance as an important measure, we have extended i) dGSIs for the multivariate and/or functional outputs (including spatiotemporal models and dynamic models), and ii) dependent sensitivity indices for single-response models so as to cope with every subset of NIVs. Such indices are also well suited for models with both discrete and continuous variables. Consistent estimators of such indices and their asymptotic distributions are provided. It is worth noting that such conditional distributions are also relevant for commuting the variance-based Shapley effects [47].

Analytical test cases confirmed that the first-order index of any subset of inputs is less than its total index, as expected. In the case of the dynamic model considered, it came out that the second-type dGSIs, which account for the correlations among the components of sensitivity functionals, and the first-type dGSIs can give different rankings of input variables. Therefore, we should prefer the second-type dGSIs in practice. Moreover, it came out that the sum of the main and interaction indices can be greater than one. In the next works, it will be interesting to investigate a new approach for which the main and interactions indices sum up to one. Also, the derivation of MSEs of the estimators of covariances of sensitivity functionals is quite interesting.

Since the computations of the effects of all the subsets of inputs rely on $2^{\frac{d-d_1}{2}} \le R_{\min} 2^{d-d_1}$ ERs, which grow exponentially with respect to $d(1-a)$, it becomes difficult to deploy such efficient approaches in higher dimensions under the following conditions: $a = d_1/d \ll 1$ and all the $K$ groups of NIVs are important, which is not often the case.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/math13050766/s1, Derivation of estimators of covariances and dGSIs.

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

**Data Availability Statement:** No data sets are used in this paper. This paper provides theoretical results with analytical test cases. Simulated data are in the paper.

**Acknowledgments:** We would like to thank the two referees and the associate editor for their comments and suggestions that have helped improving this paper.

**Conflicts of Interest:** The author has no conflicts of interest to declare regarding the content of this paper.

# Appendix A. Proof of Lemma 1

Consider the variable $Y_k = \Phi^{-1}(F_k(X_k))$ if $X_k$ is continuous and $Y_k = \Phi^{-1}(\tau_{F_k}(X_k, U_k))$ otherwise. It is well known that $Y_k$ follows the standard normal distribution with $k = 1, \dots, d$, and $\mathbf{Y} = (Y_j, \mathbf{Y}_{\sim j})$ has the same copula as $\mathbf{X}$ [26], as $\Phi^{-1} \circ F_k$ (resp. $\Phi^{-1} \circ \tau_{F_k}$)
is a strictly increasing transformation on the range of $X_k$. Therefore, $\mathbf{Y} \sim \mathcal{N}_d(\mathbf{0}, \mathcal{R})$. Knowing that the dependency function of $\mathbf{Y}$ is given by $(Y_j, \mathbf{Y}_{\sim j}) = \mathcal{L}[Y_j, \mathbf{Z}^T]^T$ (see [16]), the result follows using the inverse transformation of the form $X_k = F_k^{-1} \circ \Phi(Y_k)$.

## Appendix B. Proof of Lemma 2

Using the same reasoning as in Appendix A, we can see that $\mathbf{Y} \sim \mathbf{t}_d(\mathbf{v}, \mathbf{0}, \mathcal{R})$ with $Y_k = T_v^{-1}(F_k(X_k))$ for continuous variable $X_k$ and $Y_k = T_v^{-1}(\tau_{F_k}(X_k, U_k))$ otherwise. We then have to derive the dependency function of $\mathbf{Y}$ to obtain the result, and it is performed below. As $\mathbf{Y} \sim \mathbf{t}_d(\mathbf{v}, \mathbf{0}, \mathcal{I}) \iff \mathcal{L}\mathbf{Y} \sim \mathbf{t}_d(\mathbf{v}, \mathbf{0}, \mathcal{R})$, with $\mathcal{L}$ being the Cholesky factor of $\mathcal{R}$, the result holds knowing that [20]

$$ \left[ \begin{array}{c} Y_{w_1} = \sqrt{\frac{\nu + Y_j^2}{\nu + 1}} Z_{w_1} \\ \vdots \\ Y_{w_{d-1}} = \sqrt{\frac{(\nu + Y_j^2) \prod_{k=1}^{d-2} (\nu + k + (Z_{w_k})^2)}{\prod_{k=1}^{d-1} (\nu + k)}} Z_{w_{d-1}} \end{array} \right] \sim \mathbf{t}_d(\mathbf{v}, \mathbf{0}, \mathcal{I}). $$

## Appendix C. Proof of Lemma 3

Consider any measurable and integrable function $h: \mathbb{R}^n \to \mathbb{R}$. It is known from Proposition 1 that

$$ \mathbf{X}_{\mathbf{w}_{\sim 1,k}} \stackrel{d}{=} r_{w_{1,k}}(\mathbf{X}_{w_{1,k}}, \mathbf{Z}_{\mathbf{v}_k}, \mathbf{Z}_{\mathbf{w}_{\sim 1,k} \setminus \mathbf{v}_k}) \implies \mathbf{X}_{\sim \mathbf{u}_k} \stackrel{d}{=} r_{\mathbf{u}_k}(\mathbf{X}_{\mathbf{u}_k}, \mathbf{Z}_{\mathbf{w}_{\sim 1,k} \setminus \mathbf{v}_k}). $$

Using $\mathbf{Y} := (\mathbf{X}_{\pi_1 \setminus u_0}, \mathbf{X}_s \setminus \{w_{1,k_1}, \dots, w_{1,k_m}\}, \mathbf{Z}_{\mathbf{w}_{\sim 1,k_1} \setminus \mathbf{v}_{k_1}}, \dots, \mathbf{Z}_{\mathbf{w}_{\sim 1,k_m} \setminus \mathbf{v}_{k_m}}, \mathbf{Z}_{\mathbf{w}_{\sim 1} \setminus \{\mathbf{w}_{\sim 1,k_1}, \dots, \mathbf{w}_{\sim 1,k_m}\}})$ and the fact that the components of $(\mathbf{X}_{\pi_1}, \mathbf{X}_s, \mathbf{Z}_{\mathbf{w}_{\sim 1}})$ are independent, we can write

$$ \begin{aligned} & \mathbb{E}[h(g(\mathbf{X}_{\pi_1}, \mathbf{X}_s, \mathbf{Z}_{\mathbf{w}_{\sim 1}}, \theta)) | \mathbf{X}_{u_0}, \mathbf{X}_{w_{1,k_1}}, \mathbf{Z}_{\mathbf{v}_{k_1}}, \dots, \mathbf{X}_{w_{1,k_m}}, \mathbf{Z}_{\mathbf{v}_{k_m}}] \\ &= \mathbb{E}_{\mathbf{Y}}[h(f(\mathbf{X}_{\pi_1}, \mathbf{X}_{w_{1,2}}, r_{w_{1,2}}(\mathbf{X}_{w_{1,2}}, \mathbf{Z}_{\mathbf{w}_{\sim 1,2}}), \dots, \mathbf{X}_{w_{1,K}}, r_{w_{1,K}}(\mathbf{X}_{w_{1,K}}, \mathbf{Z}_{\mathbf{w}_{\sim 1,K}}), \theta))] \\ &= \mathbb{E}_{\mathbf{Y}}[h(f(\mathbf{X}_{\pi_1}, \mathbf{X}_{w_{1,k}}, r_{w_{1,k}}(\mathbf{X}_{w_{1,k}}, \mathbf{Z}_{\mathbf{v}_k}, \mathbf{Z}_{\mathbf{w}_{\sim 1,k} \setminus \mathbf{v}_k}), k = 2, \dots, K, \theta))] \\ &\stackrel{d}{=} \mathbb{E}_{\mathbf{Y}}[h(f(\mathbf{X}_{\pi_1}, \mathbf{X}_{\mathbf{u}_k}, r_{\mathbf{u}_k}(\mathbf{X}_{\mathbf{u}_k}, \mathbf{Z}_{\mathbf{w}_{\sim 1,k} \setminus \mathbf{v}_k}), k = 2, \dots, K, \theta))] \\ &= \mathbb{E}[h(f(\mathbf{X}, \theta)) | \mathbf{X}_{u_0}, \mathbf{X}_{\mathbf{u}_k}, \dots, \mathbf{X}_{\mathbf{u}_{k_m}}]. \end{aligned} $$

## Appendix D. Proof of Lemma 4

For Equation (15), at the end of the first step (i.e., $e_0 = 1$), $\mathcal{B}_k$ contains super-sets of $\{j_k\}$ for all $j_k \in \pi_k$ of the form $(j_k, t_k)$ with $t_k \subseteq \pi_k \setminus \{j_k\}$. Indeed, for two super-sets $(j_{k_1}, t_{k_1})$ and $(j_{k_2}, t_{k_2})$ of $\mathcal{B}_k$, we must have

$$ \{j_{k_1}, v_{1,k_1}, \dots, v_{j_{k_1}}\} \neq \{j_{k_2}, v_{1,k_2}, \dots, v_{j_{k_2}}\} \text{ for all } j \in \{0, \dots, j_{0,k}\}, \text{ and} $$

$$ \{v_{j,k_1}, \dots, v_{d_k,k_1}\} \neq \{v_{j,k_2}, \dots, v_{d_k,k_2}\} \text{ for all } j \in \{j_{0,k} + 1, \dots, d_k\}, $$

leading to $\{u \subseteq \pi_k : |u| = 1\} \subseteq \mathcal{B}_k$.

Secondly, when $e_0 = 2$ (from iteration $d_k + 1$ to $\frac{d_k(d_k-1)}{2}$), we add the super-sets of $\{j_{k_1}, j_{k_2}\}$ of the form $(j_{k_1}, j_{k_2}, t_{k_1}, t_{k_2})$, which were not in $\mathcal{B}_k$ at the end of the first step $e_0 = 1$. As for the two new super-sets, that is, $\{j_{k_1}, j_{k_2}, t_{k_1}, k_2\}$, $\{j_{k_3}, j_{k_4}, t_{k_3}, k_4\}$, we must have $\{j_{k_1}, j_{k_2}\} \neq \{j_{k_3}, j_{k_4}\}$, $\{j_{k_1}, j_{k_2}\} \neq \{j_{k_1}, v_{1,k_1}\}$, and $\{j_{k_1}, v_{1,k_1}\} \neq \{j_{k_3}, j_{k_4}\}$; the first two steps allow for obtaining $\{u \subseteq \pi_k : 1 \le |u| \le 2\} \subseteq \mathcal{B}_k$.
Thirdly, we repeat that procedure up to $e_0 = j_{0,k} - 1$ to obtain the super-sets of $\{j_{k_1}, \dots, j_{j_{0,k}-1}\}$ and $\{u \subseteq \pi_k : 1 \le |u| \le j_{0,k}-1\} \subseteq \mathcal{B}_k$. These operations are possible because $(\binom{d_k}{e_0}) \le (\binom{d_k}{j_{0,k}})$ for all $e_0 = 1, \dots, j_{0,k}-1$, and we avoid permutations $(\mathbf{w}_k)$ that bring replicated sets in both $\mathcal{B}_k$ and $\mathcal{E}_k$.

Fourthly, the iterations $(\binom{d_k}{j_{0,k}-1}) < i \le (\binom{d_k}{j_{0,k}})$ (when possible) aim to add the remaining subsets of $j_{0,k}$ elements.

Fifthly, we have $\{u \subseteq \pi_k : |u| = j_{0,k} + 1\} \subseteq \mathcal{B}_k$ because for any $\mathbf{v}_1 \subseteq \pi_k$ with $|\mathbf{v}_1| = j_{0,k} + 1$, there exists $\mathbf{w}_k^* \in \mathcal{P}_k$ such that $\{w_{j_{0,k}+2,k}^*, \dots, w_{d_k,k}^*\} \cap \mathbf{v}_1 = \emptyset$. Indeed, $\mathbf{w}_k^*$ was added in $\mathcal{P}_k$ when constructing all the subsets $u \subseteq \pi_k$ with $|u| = d_k - j_{0,k} - 1 < j_{0,k}$ thanks to $\mathcal{E}_k$ and the fact that $(\binom{d_k}{j_{0,k}+1}) = (\binom{d_k}{d_k - j_{0,k}-1})$. Thus, $\mathbf{v}_1 = \{w_{1,k}^*, \dots, w_{j_{0,k}+1,k}^*\}$. Finally, we use the same reasoning to obtain the results.

Equation (16) is obvious by construction (see Algorithm 1).

## Appendix E. Proof of Theorem 1

First, for $\mathbf{u}_k \subseteq \pi_k$ with $|\mathbf{u}_k| > 0$, there exists $\mathbf{w}_k^* \in \mathcal{P}_k$ such that $\mathbf{u}_k = \{w_{1,k}^*, \dots, w_{|\mathbf{u}_k|,k}\}$ according to Lemma 4. Lemma 3 ensures the determination of the distribution of $f$ given $\mathbf{X}_{\mathbf{u}_k}$ using $g$ associated with $\mathbf{w}_k^*$.

Second, for $u := (u_0, \mathbf{u}_2, \dots, \mathbf{u}_K)$, where $u_0 \subseteq \pi_1$, $\mathbf{u}_k \subseteq \pi_k$, and $|\mathbf{u}_k| > 0$ with $k = 2, \dots, K$, there exists only one permutation $\mathbf{w}_k^* \in \mathcal{P}_k$ such that $\mathbf{u}_k = \{w_{1,k}^*, \dots, w_{|\mathbf{u}_k|,k}\}$, $\forall k = 2, \dots, K$. As only one representation of $f$ associated with $\mathbf{w}_k^*$, $k = 2, \dots, K$ allows for determining the distribution of $f$ given $\mathbf{X}_{\mathbf{u}_k}$, and $|\mathcal{P}_k| = (\binom{d_k}{j_{0,k}})$, then $R_{\min} := \prod_{k=2}^K (\binom{d_k}{j_{0,k}})$ different representations of $f$ are needed to obtain the distribution of $f$ given $\mathbf{X}_{\mathbf{u}_k}$ for all $u \subseteq \{1, \dots, d\}$. The result follows because $R_{\min}$ is the highest number of possibilities of $\{\mathbf{u}_k \subseteq \pi_k, k = 2, \dots, K : |\mathbf{u}_k| = j_{0,k}\}$, and other possibilities are in the $R_{\min}$ representations (Lemma 4).

## Appendix F. Proof of Proposition 2

The proofs are straightforward. The results rely on the Hoeffding decomposition of an equivalent representation of $f$ and the fact that for two positive semi-definite matrices $A_1, A_2$, the Loewner partial ordering, that is, $A_1 \preceq A_2$, implies $\text{Tr}(A_1) \le \text{Tr}(A_2)$ and $||A_1||_F \le ||A_2||_F$. See [16] for more details.

## Appendix G. Proof of Theorem 2

Since $\theta_1, \theta_2$ are not random quantities, the proofs of Points (i)–(iii) are straightforward. Indeed, by expanding the above expressions of the estimators, we obtain unbiased estimators, and by applying the law of large numbers, we obtain consistent estimators. Detailed similar proofs can be found in [6,16,44].

## Appendix H. Proof of Corollary 1

Firstly, the proposed estimators are unbiased and consistent by applying Theorem 2 where $\theta_1 = \theta_2 = \theta_0$ is a constant. Secondly, the MVU properties are due to the symmetric properties of the kernels used [48]. Indeed, each kernel remains unchanged when one permutes its arguments. More details can be found in [44] (Theorems 2 and 3).

## Appendix I. Proof of Theorem 3

The results for the consistency are deduced from Corollary 3, as $\Theta = \{\theta_0\}$.

For the asymptotic distribution of Points (i)–(ii), the derivation of the results is similar to the one provided in [45] (Theorem 6), under the condition $m/M \to 0$.
# Appendix J. Equivalent Representations of the Function Used in Section 6.3

The $R_{min} = 6$ equivalent representations of $f(\mathbf{X})$ are given as follows:

$$
\begin{align*}
g_1 &:= f(X_1, r_1(X_1, Z_2, Z_3), \mathbf{X}_{\pi_1}, X_9, r_9(X_9, Z_{10})) \quad \text{for } X_1, X_9, (X_1, X_2), (X_1, X_9), \\
    &\phantom{:=} X_k, (X_1, X_k), \forall k \in \pi_1, \dots; \\
g_2 &:= f(X_2, r_2(X_2, Z_3, Z_1), \mathbf{X}_{\pi_1}, X_9, r_9(X_9, Z_{10})) \quad \text{for } X_2, (X_2, X_3), (X_2, X_9), (X_9, X_{10}), \dots; \\
g_3 &:= f(X_3, r_3(X_3, Z_1, Z_2), \mathbf{X}_{\pi_1}, X_9, r_9(X_9, Z_{10})) \quad \text{for } X_3, (X_3, X_1), (X_3, X_9), \dots; \\
g_4 &:= f(X_1, r_1(X_1, Z_2, Z_3), \mathbf{X}_{\pi_1}, X_{10}, r_{10}(X_{10}, Z_9)) \quad \text{for } X_{10}, (X_1, X_{10}), \dots; \\
g_5 &:= f(X_2, r_2(X_2, Z_3, Z_1), \mathbf{X}_{\pi_1}, X_{10}, r_{10}(X_{10}, Z_9)) \quad \text{for } (X_2, X_{10}), \dots; \\
g_6 &:= f(X_3, r_3(X_3, Z_1, Z_2), \mathbf{X}_{\pi_1}, X_{10}, r_{10}(X_{10}, Z_9)) \quad \text{for } (X_3, X_{10}), \dots,
\end{align*}
$$

where $r_1, r_2, r_3$ are particular cases of DMs provided in Lemma 1.

## References

1. Sobol, I.M. Sensitivity analysis for non-linear mathematical models. *Math. Model. Comput. Exp.* **1993**, *1*, 407–414.

2. Lamboni, M. Multivariate sensitivity analysis: Minimum variance unbiased estimators of the first-order and total-effect covariance matrices. *Reliab. Eng. Syst. Saf.* **2019**, *187*, 67–92. [CrossRef]

3. Lamboni, M.; Monod, H.; Makowski, D. Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. *Reliab. Eng. Syst. Saf.* **2011**, *96*, 450–459. [CrossRef]

4. Gamboa, F.; Janon, A.; Klein, T.; Lagnoux, A. Sensitivity indices for multivariate outputs. *Comptes Rendus Math.* **2013**, *351*, 307–310. [CrossRef]

5. Xiao, S.; Lu, Z.; Xu, L. Multivariate sensitivity analysis based on the direction of eigen space through principal component analysis. *Reliab. Eng. Syst. Saf.* **2017**, *165*, 1–10. [CrossRef]

6. Lamboni, M. Derivative-based generalized sensitivity indices and Sobol' indices. *Math. Comput. Simul.* **2020**, *170*, 236–256. [CrossRef]

7. Perrin, T.; Roustant, O.; Rohmer, J.; Alata, O.; Naulin, J.; Idier, D.; Pedreros, R.; Moncoulon, D.; Tinard, P. Functional principal component analysis for global sensitivity analysis of model with spatial output. *Reliab. Eng. Syst. Saf.* **2021**, *211*, 107522. [CrossRef]

8. Veiga, S.D.; Wahl, F.; Gamboa, F. Local Polynomial Estimation for Sensitivity Analysis on Models With Correlated Inputs. *Technometrics* **2009**, *51*, 452–463. [CrossRef]

9. Mara, T.A.; Tarantola, S. Variance-based sensitivity indices for models with dependent inputs. *Reliab. Eng. Syst. Saf.* **2012**, *107*, 115–121. [CrossRef]

10. Kucherenko, S.; Tarantola, S.; Annoni, P. Estimation of global sensitivity indices for models with dependent variables. *Comput. Phys. Commun.* **2012**, *183*, 937–946. [CrossRef]

11. Hao, W.; Zhenzhou, L.; Wei, P. Uncertainty importance measure for models with correlated normal variables. *Reliab. Eng. Syst. Saf.* **2013**, *112*, 48–58. [CrossRef]

12. Chastaing, G.; Gamboa, F.; Prieur, C. Generalized Hoeffding-Sobol' decomposition for dependent variables—Applications to sensitivity analysis. *Electron. J. Stat.* **2012**, *6*, 2420–2448. [CrossRef]

13. Kucherenko, S.; Klymenko, O.; Shah, N. Sobol' indices for problems defined in non-rectangular domains. *Reliab. Eng. Syst. Saf.* **2017**, *167*, 218–231. [CrossRef]

14. Mara, T.A.; Tarantola, S.; Annoni, P. Non-parametric methods for global sensitivity analysis of model output with dependent inputs. *Environ. Model. Softw.* **2015**, *72*, 173–183. [CrossRef]

15. Tarantola, S.; Mara, T.A. Variance-based sensitivity indices of computer models with dependent inputs: The fourier amplitude sensitivity test. *Int. J. Uncertain. Quantif.* **2017**, *7*, 511–523. [CrossRef]

16. Lamboni, M.; Kucherenko, S. Multivariate sensitivity analysis and derivative-based global sensitivity measures with dependent variables. *Reliab. Eng. Syst. Saf.* **2021**, *212*, 107519. [CrossRef]

17. Lamboni, M. On exact distribution for multivariate weighted distributions and classification. *Methodol. Comput. Appl. Probab.* **2023**, *25*, 1–41. [CrossRef]

18. Lamboni, M. Kernel-based Measures of Association Between Inputs and Outputs Using ANOVA. *Sankhya A* **2024**, *86*, 790–826. [CrossRef]

19. Skorohod, A.V. On a representation of random variables. *Theory Probab. Appl.* **1976**, *21*, 645–648.
20. Lamboni, M. Efficient dependency models: Simulating dependent random variables. *Math. Comput. Simul.* **2021**, *200*, 199–217.
[CrossRef]

21. Ferguson, T. *Mathematical Statistics: A Decision Theoretic Approach*; Academic Press: New York, NY, USA, 1967.

22. Rschendorf, L. Stochastically ordered distributions and monotonicity of the oc-function of sequential probability ratio tests. *Ser. Stat.* **1981**, *12*, 327–338. [CrossRef]

23. Rschendorf, L. Stochastic ordering of risks, influence of dependence and a.s. constructions. In *Advances on Models, Characterizations and Applications*; Balakrishnan, N., Bairamov, I.G., Gebizlioglu, O.L., Eds.; CRC Press: Boca Raton, FL, USA, 2005.

24. Rschendorf, L. On the distributional transform, Sklar's theorem, and the empirical copula process. *J. Stat. Plan. Inference* **2009**, *139*, 3921–3927. [CrossRef]

25. Rosenblatt, M. Remarks on a Multivariate Transformation. *Ann. Math. Statist.* **1952**, *23*, 470–472. [CrossRef]

26. Nelsen, R. *An Introduction to Copulas*; Springer: New York, NY, USA, 2006.

27. McNeil, A.J.; Frey, R.; Embrechts, P. *Quantitative Risk Management*; Princeton University Press: Princeton, NJ, USA; Oxford, UK, 2015.

28. Rosenblatt, M. Remarks on some nonparametric estimates of a density function. *Ann. Math. Stat.* **1956**, *27*, 832–837. [CrossRef]

29. Parzen, E. On estimation of a probability density function and mode. *Ann. Math. Stat.* **1962**, *33*, 1065–1076. [CrossRef]

30. Epanechnikov, V. Nonparametric estimation of a multidimensional probability density. *Theory Probab. Appl.* **1969**, *14*, 153–158.
[CrossRef]

31. Silverman, B. *Density Estimation for Statistics and Data Analysis*; Chapman & Hall: New York, NY, USA, 1986.

32. Clayton, D.G. A Model for Association in Bivariate Life Tables and Its Application in Epidemiological Studies of Familial Tendency in Chronic Disease Incidence. *Biometrika* **1978**, *65*, 141–151. [CrossRef]

33. Joe, H. *Multivariate Models and Dependence Concepts*; Chapman & Hall/CRC: Boca Raton, FL, USA; London, UK; New York, NY, USA, 1997.

34. Smith, M.; Min, A.; Almeida, C.; Czado, C. Modeling Longitudinal Data Using a Pair-Copula Decomposition of Serial Dependence. *J. Am. Statist. Assoc.* **2010**, *105*, 1467–1479. [CrossRef]

35. Durante, F.; Sempi, C. *Principles of Copula Theory*; CRC/Chapman & Hall: London, UK, 2015.

36. Koenker, R.; Bassett, G. Regression quantiles. *Econometrica* **1978**, *46*, 33–50. [CrossRef]

37. Truong, Y.K. Asymptotic Properties of Kernel Estimators Based on Local Medians. *Ann. Stat.* **1989**, *17*, 606–617. [CrossRef]

38. Hendricks, W.; Koenker, R. Hierarchical Spline Models for Conditional Quantiles and the Demand for Electricity. *J. Am. Statist. Assoc.* **1992**, *87*, 58–68. [CrossRef]

39. Koenker, R.; Ng, P.; Portnoy, S. Quantile smoothing splines. *Biometrika* **1994**, *81*, 673–680. [CrossRef]

40. Koenker, R.; Hallock, K.F. Quantile Regression. *J. Econ. Perspect.* **2001**, *15*, 143–156. [CrossRef]

41. Jr., G.B.; Koenker, R. An Empirical Quantile Function for Linear Models with iid Errors. *J. Am. Statist. Assoc.* **1982**, *77*, 407–415.

42. Koenker, R. *Quantile Regression*; Cambridge University Press: Cambridge, UK, 2005.

43. Takeuchi, I.; Le, Q.V.; Sears, T.D.; Smola, A.J. Nonparametric Quantile Estimation. *J. Mach. Learn. Res.* **2006**, *7*, 1231–1264.

44. Lamboni, M. Uncertainty quantification: A minimum variance unbiased (joint) estimator of the non-normalized Sobol' indices. *Stat. Pap.* **2018**, *61*, 1939–1970. [CrossRef]

45. Lamboni, M. Weak derivative-based expansion of functions: ANOVA and some inequalities. *Math. Comput. Simul.* **2022**, *194*, 691–718. [CrossRef]

46. Lamboni, M. Global sensitivity analysis: An efficient numerical method for approximating the total sensitivity index. *Int. J. Uncertain. Quantif.* **2016**, *6*, 1–17. [CrossRef]

47. Owen, A.B. Sobol' Indices and Shapley Value. *Siam/Asa J. Uncertain. Quantif.* **2014**, *2*, 245–251. [CrossRef]

48. Sugiura, N. Multisample and multivariate nonparametric tests based on U-statistics and their asymptotic efficiencies. *Osaka J. Math.* **1965**, *2*, 385–426.

**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.