Bayesian estimation of probabilistic sensitivity measures

TL;DR

This paper introduces a Bayesian framework for estimating probabilistic sensitivity measures in computer experiments, enabling uncertainty quantification with minimal additional simulations, especially useful for complex, high-dimensional models.

Contribution

It presents novel Bayesian estimators, including piecewise constant priors and non-parametric density methods, for efficient sensitivity analysis at low computational cost.

Findings

Bayesian methods provide accurate uncertainty quantification.

Proposed estimators outperform traditional point estimators.

Effective in synthetic and real-world complex simulations.

Abstract

Computer experiments are becoming increasingly important in scientific investigations. In the presence of uncertainty, analysts employ probabilistic sensitivity methods to identify the key-drivers of change in the quantities of interest. Simulation complexity, large dimensionality and long running times may force analysts to make statistical inference at small sample sizes. Methods designed to estimate probabilistic sensitivity measures at relatively low computational costs are attracting increasing interest. We propose a fully Bayesian approach to the estimation of probabilistic sensitivity measures based on a one-sample design. We discuss, first, new estimators based on placing piecewise constant priors on the conditional distributions of the output given each input, by partitioning the input space. We then present two alternatives, based on Bayesian non-parametric density estimation, which bypass the need for predefined partitions. In all cases, the Bayesian paradigm guarantees the quantification of uncertainty in the estimation process through the posterior distribution over the sensitivity measures, without requiring additional simulator evaluations. The performance of the proposed methods is compared to that of traditional point estimators in a series of numerical experiments comprising synthetic but challenging simulators, as well as a realistic application. $\textit{An Updated Version of the Manuscript is Forthcoming in Statistics and Computing.}$

Tables

Table 1. Table 1: Some probabilistic sensitivity measures

Measure	$\zeta ({\mathbb{P}}_Y,{\mathbb{P}}_{Y|X^i})$	${\xi }_i$
${\eta }_i$	${\left(𝔼\left[Y|X^i\right]-𝔼\left[Y\right]\right)}^2/𝕍\left[Y\right]$	$𝕍\left[𝔼\left(Y|X^i\right)\right]/𝕍\left[Y\right]$
${\delta }_i$	$\frac{1}{2}{\int }_{𝒴}|f_{Y|X^i}\left(y\right|X^i)-f_Y\left(y\right)|\text{d}y$	$\frac{1}{2}𝔼\left[{\int }_{𝒴}\right|f_{Y|X^i}\left(y\right|X^i)-f_Y\left(y\right)\left|\text{d}y\right]$
${\beta }_i$	${sup}_{y\in 𝒴}|F_{Y|X^i}\left(y\right|X^i)-F_Y\left(y\right)|$	$𝔼\left[{sup}_{y\in 𝒴}\right|F_{Y|X^i}\left(y\right|X^i)-F_Y\left(y\right)\left|\right]$

Table 2. Table 2: Analytical values of η i subscript 𝜂 𝑖 \eta_{i} , δ i subscript 𝛿 𝑖 \delta_{i} and β i subscript 𝛽 𝑖 \beta_{i} for the two test simulators used in this section

	2-input		21-input
Measure	$X_1$	$X_2$	$X_1-X_7$	$X_8-X_{14}$	$X_{15}-X_{21}$
${\eta }_i$	0.496	0.496	0.309	0.064	0.092
${\delta }_i$	0.315	0.315	0.212	0.084	0.102
${\beta }_i$	0.289	0.289	0.205	0.083	0.101

Table 3. Table 3: Simulator inputs for the LevelE code. U ​ ( ⋅ , ⋅ ) 𝑈 ⋅ ⋅ U(\cdot,\cdot) and L ​ U ​ ( ⋅ , ⋅ ) 𝐿 𝑈 ⋅ ⋅ LU(\cdot,\cdot) stand for the uniform and log-uniform distributions respectively

Input	Meaning	Distribution
$X_1$	Containment time	U(100,1000)
$X_2$	Iodine Leach rate	LU(${10}^{-3},{10}^{-2}$)
$X_3$	Neptunium chain Leach rate	LU(${10}^{-6},{10}^{-5}$)
$X_4$	Iodine retention factor (1st layer)	LU(${10}^{-3},{10}^{-1}$)
$X_5$	Geosphere water velocity 1st layer	U(100,500)
$X_6$	Geosphere Length 1st layer	U(1,5)
$X_7$	Factor to compute Neptunium retention coefficients Layer 1	U(3,30)
$X_8$	water velocity in geosphere’s 2nd layer	LU(${10}^{-2},{10}^{-1}$)
$X_9$	Length of geosphere’s 2nd layer	U(50,200)
$X_{10}$	Retention factor for I (2nd layer)	U(1,5)
$X_{11}$	Factor to compute Neptunium retention coefficients Layer 2	U(3,30)
$X_{12}$	Stream flow rate	LU(${10}^5,{10}^7$)

Equations

$${\xi }_i:=\mathbb{E}[\zeta ({\mathbb{P}}_Y,{\mathbb{P}}_{Y|X^i})]$$

$${\hat{\xi }}_i=\sum _{m=1}^M{\mathbb{P}}_{X^i}({\mathcal{X}}_m^i){\hat{\zeta }}_m^i,$$

$${\hat{\eta }}_i^{\star }=\sum _{m=1}^M\frac{n_m^i}{n}\frac{({\bar{y}}_m^i-\bar{y}{)}^2}{s_y^2}.$$

$${\hat{\delta }}_i^{\star }=\sum _{m=1}^M\frac{n_m^i}{n}{\int }_{\mathcal{Y}}|{\hat{f}}_Y^{\star }(y)-{\hat{f}}_m^i(y)|\text{d}y,$$

$${\hat{\delta }}_i^{\diamond }=\sum _{m=1}^M\frac{n_m^i}{n}\sum _{t=1}^{T_m^i}\left({\hat{F}}_m^i({\hat{b}}_m^i(t))-{\hat{F}}_m^i({\hat{a}}_m^i(t))\right)-\left({\hat{F}}_Y({\hat{b}}_m^i(t))-{\hat{F}}_Y({\hat{a}}_m^i(t))\right).$$

$${\hat{\beta }}_i^{\diamond }=\sum _{m=1}^M\frac{n_m^i}{n}\underset{j\in \{1,\dots ,n\}}{\max }\left|{\hat{F}}_Y(y_j)-{\hat{F}}_m^i(y_j)\right|,$$

$${\hat{F}}_Y(y)=\frac{1}{n}\sum _{j=1}^n\text{mathds}{1}_{(-\infty ,y_j]}(y);{\hat{F}}_m^i(y)=\frac{1}{n_m^i}\sum _{y_j\in {\mathbf{y}}_m^i}\text{mathds}{1}_{(-\infty ,y_j]}(y),$$

$${\hat{\eta }}_i^{◊}=\sum _{m=1}^M\frac{n_m^i}{n}\frac{{\left({\int }_{\mathcal{Y}}{\hat{F}}_m^i(y)-{\hat{F}}_Y(y)\text{d}y\right)}^2}{{\hat{\sigma }}_Y^2}.$$

$$Y|({\mathbb{P}}_m^i,X^i\in {\mathcal{X}}_m^i)\sim {\mathbb{P}}_m^i;{\mathbb{P}}_m^i\stackrel{iid}{\sim }\mathcal{DP}(\alpha G),$$

$${\mathbb{P}}_Y(\cdot |{\mathbb{P}}_{1:M}^i)=\sum _{m=1}^M{\mathbb{P}}_m^i(\cdot )\mathbb{P}(X^i\in {\mathcal{X}}_m^i)=\frac{1}{M}\sum _{m=1}^M{\mathbb{P}}_m^i(\cdot )\text{.}$$

$${\mathbb{P}}_Y(\cdot |\alpha G)=\frac{1}{M}\sum _{m=1}^M\int {\mathbb{P}}_m^i(\cdot )\text{d}\mathcal{DP}({\mathbb{P}}_m^i|\alpha G)=\int \mathbb{P}(\cdot )\text{d}\mathcal{DP}(\mathbb{P}|\alpha G),$$

$$Y|(X^i\in {\mathcal{X}}_m^i,{\mathbb{P}}_m^i)\sim {\mathbb{P}}_m^i;{\mathbb{P}}_m^i|Data\stackrel{ind}{\sim }\mathcal{DP}\left((\alpha +n_m^i){\tilde{G}}_m^i\right),$$

$${\tilde{G}}_m^i=\mathbb{E}[{\mathbb{P}}_m^i|Data]=\frac{\alpha }{\alpha +n_m^i}G+\frac{n_m^i}{\alpha +n_m^i}\sum _{y\in {\mathbf{y}}_m^i}\frac{1}{n_m^i}{\delta }^{Dirac}(y).$$

$${\mathbb{P}}_Y(\cdot |\alpha G,Data)=\frac{1}{M}\sum _{m=1}^M\int {\mathbb{P}}_m^i(\cdot )\text{d}\mathcal{DP}({\mathbb{P}}_m^i|(\alpha +n_m^i){\tilde{G}}_m^i),$$

$$\tilde{G}=\frac{\alpha }{\alpha +n}G+\frac{n}{\alpha +n}{\hat{\mathbb{P}}}_n,$$

$${\tilde{\xi }}_i=\sum _{m=1}^M\frac{n_m^i}{n}\zeta (\tilde{G},{\tilde{G}}_m^i)$$

$${\hat{\xi }}_i^{Bb}=\frac{1}{S}\sum _{s=1}^S{\hat{\xi }}_i^{Bb,s}$$

$${\tilde{Y}}_{j+1}^i|\left({\tilde{y}}_{n_m^i+1:j}^i,Data\right)\sim \frac{\alpha }{\alpha +j}G+\frac{j}{\alpha +j}{\hat{\mathbb{P}}}_j\forall j\ge n_m^i.$$

$$Y=\frac{X^1}{X^1+X^2},$$

$$Y=\sum _{i=1}^{21}{a}_i{X}^i,$$

$$f_{X^i,Y}(\cdot ,\cdot )|P=\int \mathcal{K}(\cdot ,\cdot |\theta )\text{d}P(\theta ),$$

$$P(\theta )=\sum _{\ell =1}^{\infty }{w}_{\ell }{\delta }^{Dirac}({\theta }_{\ell }),$$

$$W_1=V_1,W_{\ell }=V_{\ell }\prod _{{\ell }^{\prime }=1}^{\ell -1}(1-V_{{\ell }^{\prime }})$$

$$f_{X^i,Y}(\cdot ,\cdot )|P=\sum _{\ell =1}^{\infty }{w}_{\ell }\cdot \mathcal{N}(\cdot ,\cdot |{\mu }_{\ell },{\Sigma }_{\ell }),$$

$$f_{Y|X}(y|x)=\int \mathcal{K}(y|x,\theta )\text{d}{P}_x(\theta ),$$

$${\mathbb{P}}_x(\theta )=\sum _{\ell =1}^{\infty }{w}_{\ell }(x){\delta }_{{\theta }_{\ell }(x)}.$$

$$w_{\ell }(x)=\frac{{\omega }_{\ell }\mathcal{K}(x|{\psi }_{\ell })}{{\sum }_{{\ell }^{\prime }=1}^{\infty }{\omega }_{{\ell }^{\prime }}\mathcal{K}(x|{\psi }_{{\ell }^{\prime }})}.$$

$$f_{Y|X^i}(y|x^i)=\sum _{\ell =1}^{\infty }{w}_{\ell }(x^i)\mathcal{N}(y|a_{\ell }+b_{\ell }{x}^i,{\sigma }_{\ell }),$$

$${\text{RMSE}}_i(n)\approx \sqrt{\frac{{\sum }_{s=1}^S{\left({\hat{\xi }}_s^i(n)-{\xi }_i\right)}^2}{S}}$$

$${\hat{\eta }}^C=\frac{1}{S}\sum _{s=1}^S{\eta }^{C,s},\text{with }{\eta }^{C,s}=\sum _{m=1}^M\frac{n_m}{N}\frac{{\left({\bar{y}}_m^{C,s}-\bar{y}\right)}^2}{s_y^2},$$

Full text

An Updated Version of the Manuscript is Forthcoming in Statistics and Computing.

1 Introduction

The use of computer simulations is becoming increasingly important in broad areas of science (Lin et al., 2010; Wong et al., 2017). High-fidelity mathematical models allow analysts to perform virtual (or in silico) experiments on complex natural or societal phenomena of interest (see Smith et al., 2009, among others). Predictions are often used to support policy-making. However, the level of sophistication of the models is often too high for analytical solutions to be available. In these cases, the only way to obtain a quantitative solution may be to encode complex mathematical equations in a computer software; so that the input-output mapping remains a black-box to the analyst. It then becomes important to carefully design and execute the computer experiment. The design and analysis of computer experiments (DACE) has entered the statistical literature with the seminal work of Sacks et al. (1989) (see also the monographs of Santner et al. (2003); Kleijnen (2008)). Since then, researchers have studied the creation of space-filling designs (Pronzato and Müller, 2012; He, 2017), the calibration of computer codes with real data (Tuo and Wu, 2015), their emulation (Marrel et al., 2012), the quantification of uncertainty in their output (Oakley and O’Hagan, 2002; Ghanem et al., 2016) and their sensitivity analysis (Oakley and O’Hagan, 2004; Borgonovo et al., 2014). These areas are intertwined. A given design may allow, for instance, not only an uncertainty quantification, but also the creation of an emulator and the analysis of sensitivity.

Probabilistic (or global) sensitivity measures are an indispensable complement of uncertainty quantification, as they highlight which areas should be given priority when planning data collection or further modelling efforts. International agencies such as the US Environmental Protection Agency (U.S. Environmental Protection Agency, 2009) or the British National Institute for Health Care Excellence (NICE, 2013) and the European Commission (2009) have issued guidelines recommending the use of probabilistic sensitivity analysis methods as the gold standard for ensuring reliability and transparency when using the output of a computer code for decision-making under uncertainty. Over the years, several probabilistic sensitivity measures have been proposed. Different measures enjoy alternative properties making them preferable in different contexts and for different purposes. We recall regression-based (Helton and Sallaberry, 2009), variance-based (Saltelli and Tarantola, 2002; Jiménez Rugama and Gilquin, 2018) and moment-independent measures (Borgonovo et al., 2014), all of which offer alternative ways to quantify the degree of statistical dependence between the simulator inputs and the output. A transversal issue in realistic applications is that analytical expressions of these measures are unavailable and analysts must resort to estimation. This, however, is a challenging task, especially for simulators with a high number of inputs (the curse of dimensionality) or with long running times (high computational cost).

Recent results (Strong et al., 2012; Strong and Oakley, 2013) evidence the one-sample (or given-data) approach as an attractive design, which allows analysts to estimate global sensitivity measures from a single probabilistic sensitivity analysis sample, i.e., a sample generated for propagating uncertainty in the simulator. Thus, a one-sample approach has a nominal cost equal to the sample size and independent of the number of inputs, a feature that potentially reduces the impact of the curse of dimensionality. Related works such as Strong et al. (2012); Strong and Oakley (2013); Plischke et al. (2013); Strong et al. (2014, 2015) and Borgonovo et al. (2016) provide advances on theoretical and numerical aspects of the methodology, while at the same time, evidencing some open research questions. One-sample estimation procedures are closely related to scatter-plot smoothing, where partitioning of the covariate space plays a central role (Hastie and Tibshirani, 1990). Strong and Oakley (2013) show that the choice of partition size affects estimation, especially when the sample size is small (see Fig. 1 of (Strong and Oakley, 2013, p. 759)) . In the literature, some heuristics for determining a partition selection strategy which is optimal in some sense have been proposed, but finding a universally valid heuristic seems out of reach (see Appendix A.1 for numerical experiments illustrating this issue). Moreover, the literature is concerned with point estimators and uncertainty regarding the estimated value of a sensitivity measure is not an intrinsic part of the analysis. Because most one-sample estimators are consistent (in the frequentist sense), an accurate estimation of the error is often overlooked. However, especially at small sample sizes, it is essential for transparency that interval estimates become part of result communication (see Janon et al. (2014b) among others).

We propose to enrich the one-sample design through the use of Bayesian non-parametric (BNP) methods, aiming to reduce and even eliminate the partition selection problem, while making uncertainty in the estimates a natural ingredient. First, we extend the partition approach, using Bayesian non-parametric models to augment the output sample within each partition set, by adequately generating additional synthetic data according to two alternative schemes. The Bayesian intuition supporting these designs can be interpreted as setting a prior on the conditional distribution of the output, given that the input realization falls within a given set of the partition. We build estimators based on this intuition for variance-based, density (pdf)-based and cumulative distribution function (cdf)-based global sensitivity measures. We compare the results with given-data estimators currently in use at low sample sizes, through numerical experiments. The results show that our estimators recover the correct ranking of the inputs, while providing an appropriate quantification of the estimation uncertainty. However, the results may be strongly influenced by the partition choice. Therefore, we investigate two additional classes of estimators based on Bayesian non-parametric joint and conditional density estimation methods. These estimators eliminate the partition selection problem and, at the same time, enable error quantification. Finally, we discuss the application of all the new estimators to the global sensitivity analysis of the benchmark computer code known as LevelE (Saltelli and Tarantola, 2002). Results show that the estimators correctly identify the key drivers of uncertainty at sample sizes lower than the ones used in previous literature. Additionally, the analyst obtains a quantification of the uncertainty in the estimates in the form of a posterior distribution which can be used to determine whether the available sample is large enough to make robust conclusions about the simulator input ranking.

The reminder of the paper is organized as follows. We begin in Section 2 by introducing the framework of global sensitivity analysis and the one-sample estimation approach for probabilistic sensitivity measures. Section 3 combines Bayesian non-parametric methods and the one-sample approach to create two new partition-dependent estimators. Section 4 derives two classes of Bayesian partition-free estimators by adopting Bayesian a non-parametric density estimation approach. Section 5 presents numerical results for the LevelE code. Section 6 offers discussion and conclusions. The appendices illustrate the algorithms of the proposed estimators and discuss additional implementation details.

2 Probabilistic sensitivity analysis of computer experiments

Formally, the sensitivity analysis framework considers a multivariate mapping $g:\mathcal{X}\mapsto\mathcal{Y}$ with input space $\mathcal{X}\subseteq\mathbb{R}^{k}$ and output space $\mathcal{Y}\subseteq\mathbb{R}^{d}$, denoted as $\mathbf{y}=g(\mathbf{x})+\epsilon(\mathbf{x})$ in its more general form. In the DACE set-up, $g$ represents a set of operations performed by a computer code which processes a set $\mathbf{x}$ of inputs, resulting in a set $\mathbf{y}$ of outputs of interest. The term $\epsilon(\mathbf{x})$ represents a zero-mean error term, which is present when the simulator response is stochastic. For simplicity, we focus on deterministic univariate responses, with $\epsilon(\mathbf{x})\equiv 0$ and $d=1$. When information is not sufficient to fix the values of the inputs, it is common to “assume to have information about the factors’ probability distribution, either joint or marginal, with or without correlation, and that this knowledge comes from measurements, estimates, expert opinion…”(Saltelli and Tarantola, 2002, p. 704). We denote the input probability space by $(\mathcal{X},\mathcal{B(X)},\mathbb{P}_{\mathbf{X}})$, where $\mathbb{P}_{\mathbf{X}}$ represents the joint probability measure of $\mathbf{X}=(X^{1},\ldots,X^{k})$, assumed known. Similarly, $(\mathcal{Y},\mathcal{B(\mathcal{Y})},\mathbb{P}_{Y})$ denotes the output probability space, where $\mathbb{P}_{Y}$ represents the distribution of $Y$ induced by $\mathbb{P}_{\mathbf{X}}$ through $g$.

It has been recently shown that several probabilistic sensitivity measures frequently used in practice can be expressed as expectations of measures of discrepancy between $\mathbb{P}_{Y}$ and $\mathbb{P}_{Y|X^{i}}$. In particular, we focus on probabilistic sensitivity measures of the form:

[TABLE]

where the expectation is calculated with respect to the marginal distribution of $X^{i}$ and $\zeta$ is a pre-metric on the space of probability measures over $\mathcal{Y}$, and $\xi_{i}$ is called the probabilistic sensitivity measure of $X^{i}$ with inner operator $\zeta$ (Borgonovo et al., 2014).

Table 1 reports three probabilistic sensitivity measures encompassed by this construction, namely, the variance-based sensitivity measure ($\eta_{i}$), the density-based $\delta$-importance measure ($\delta_{i}$) and the cdf-based $\beta$-importance measure ($\beta_{i}$) (Pearson, 1905; Saltelli and Tarantola, 2002; Oakley and O’Hagan, 2004). The quantity $\eta_{i}$ represents the expected fractional reduction in the simulator output variance attained by fixing $X^{i}$ and it is equal to the popular first-order variance-based sensitivity index (Iman and Hora, 1990; Sobol’, 1993) when the simulator inputs are mutually independent Saltelli and Tarantola (2002); Liu and Owen (2006).

Most global sensitivity measures are used by analysts to quantify the strength of the statistical dependence of $Y$ on $X^{i}$. It has been noted that first-order variance-based sensitivity measures do not possess the nullity-implies-independence property, known also as Rényi’s’s postulate D for measures of statistical dependence (Rényi, 1959). In simple terms, a null value of $\eta_{i}$ does not imply that $Y$ is independent of $X^{i}$. This has led to the introduction and study of sensitivity measures that satisfy such postulate and the class of distribution-based (or moment-independent) sensitivity measures is attracting increasing attention in the literature (see e.g. Gamboa et al., 2018; Da Veiga, 2015; Rahman, 2014). As representatives, we focus on the sensitivity measures $\delta_{i}$ and $\beta_{i}$ (Table 1), which quantify the expected separation between $\mathbb{P}_{Y}$ and $\mathbb{P}_{Y|X^{i}}$ through the $L^{1}$-norm between densities and the Kolmogorov-Smirnov distance between cumulative distribution functions, respectively. Besides complying with Renyi’s postulate D, these probabilistic sensitivity measures are also invariant to monotonic transformations – a property that accelerates convergence in numerical estimation (for further details see e.g. Borgonovo et al., 2014).

As mentioned in the introduction, analytical expressions for these and other popular sensitivity measures are not available in most realistic applications, and their estimation is a prolific subject of research. We now discuss some relevant aspects of numerical estimation in the next section.

2.1 One-sample approach to the estimation of probabilistic sensitivity measures

The estimation of global sensitivity measures is a challenging task and the availability of efficient designs is crucial in realistic applications. The number of simulator evaluations necessary to estimate sensitivity measures encompassed by Eq. (1) for a simulator with $k$ simulator inputs, using a brute-force approach, would be of the order of $C=kn^{2}$ simulator runs, where $n$ denotes the sample size required for Monte Carlo uncertainty quantification. The design becomes rapidly infeasible. For instance, if $k=20$ and $n=1,000$, the $C=20,000,000$ simulator runs would require a prohibitive computational effort for most complex computer codes used in practice. However, notable advances in the literature have contributed in abating this computational burden (see Tissot and Prieur, 2015; Janon et al., 2014a, for reviews). Saltelli (2002), for instance, achieved the estimation of variance-based sensitivity measures at a cost of $C=n(k+2)$ simulator runs, while the $FAST$ method of Saltelli et al. (1999) achieves a cost of order $C=nk$.

Recently, efforts have been made towards an estimation cost independent of the number of simulator inputs, $k$. Strong and Oakley (2013), Strong et al. (2014) and Strong et al. (2015) show that value-of-information measures can be estimated from a single probabilistic sensitivity analysis sample, $\{(\mathbf{x}_{j},y_{j}):j=1,\ldots n\}$, i.e., from the Monte Carlo sample generated for uncertainty quantification, thus lowering the computational cost to $C=n$ simulator runs. Röhlig et al. (2009) and Strong et al. (2012) obtain similar results for first-order variance-based sensitivity measures and Plischke et al. (2013) extend the intuition to density-based measures. The approaches proposed in these works receive the common name of one-sample or given-data estimation methods.

One-sample methods can be seen as generalizations of the intuition developed for estimating the correlation ratio (Pearson, 1905). If $X^{i}$ is a discrete random variable then, an input-output sample of (sufficiently large) size $n$ contains repeated observations of $Y=g(X^{i},X^{-i})$, for each fixed value $X^{i}=x^{i}$, while the other factors, $X^{-i}=(X^{1},\dots,X^{i-1},X^{i+1},\dots,X^{k})$, remain random. This allows the estimation of $\zeta(\mathbb{P}_{Y},\mathbb{P}_{Y|X^{i}=x^{i}})$ directly from a sample of size $n$. For a continuous $X^{i}$, a similar result may be achieved by partitioning the support $\mathcal{X}^{i}$ of $X^{i}$ into $M$ bins $\{\mathcal{X}^{i}_{m}\}_{m=1}^{M}$. The point condition $(X^{i}=x^{i})$ is then replaced by the bin condition $(X^{i}\in\mathcal{X}^{i}_{m})$. Then, for any sensitivity measure encompassed by Eq. (1), a one-sample estimator is given by (Borgonovo et al., 2016):

[TABLE]

where $\hat{\zeta}^{i}_{m}$ may be any estimator of $\zeta(\mathbb{P}_{Y},\mathbb{P}_{Y|X^{i}\in\mathcal{X}^{i}_{m}})$. Note that by using equiprobable partition sets, $\mathbb{P}_{X^{i}}(\mathcal{X}^{i}_{m})$ should reduce to $1/M$. In practice, this partition probability is estimated by the sample proportion, $n^{i}_{m}/n$, where $n^{i}_{m}$ denotes the number of realizations for which the $i$-th input falls within the $m$-th partition set of its support. Borgonovo et al. (2016) (Theorem 2) show that, under mild conditions on the inner operator $\zeta$, a consistent version of the estimator in Eq. (2) can be obtained, if the size $M$ of the partition is chosen as a monotonically increasing function of the sample size $n$, such that $\lim\limits_{n\rightarrow\infty}\frac{n}{M(n)}=\infty$.

Indeed, the most popular one-sample estimator of $\eta_{i}$ relies on a plug-in estimator of the inner statistic, based on the output sample mean and variance, $\bar{y}$ and $s^{2}_{y}$ respectively, to estimate the marginal mean and variance of $Y$. The within cluster sample mean $\bar{y}^{i}_{m}=\dfrac{1}{n^{i}_{m}}\sum_{y\in\mathbf{y}^{i}_{m}}y$ with $\mathbf{y}^{i}_{m}=\{y_{j}:x_{j}^{i}\in\mathcal{X}^{i}_{m},j=1,2,...,n\}$ is used to estimate the conditional mean of $Y|X^{i}\in\mathcal{X}^{i}_{m}$. The final expression (see e.g. Strong et al., 2012) takes the form of Eq. (2) with:

[TABLE]

The one-sample estimator for the $\delta-$importance introduced by Plischke et al. (2013) can be written as:

[TABLE]

where $\hat{f}^{\star}_{Y}$ and $\hat{f}^{i}_{m}$ denote kernel-smoothed histograms of the full output vector $\mathbf{y}=(y_{1},\ldots y_{n})$ and the within cluster output vector $\mathbf{y}^{i}_{m}$, respectively. The authors propose a quadrature method for the numerical integration required by the $L^{1}$-norm in the inner operator, but other solutions could be used, producing similar estimators. Because estimates of this type rely on the approximation or estimation of probability density functions, we refer to them as pdf-based estimators.

Plischke and Borgonovo (2017) observe that the kernel-smoothing methods commonly involved in the calculation of pdf-based estimators may induce bias, even at high sample sizes, for simulators with a sparse output. Therefore, they introduce alternative cdf-based estimators which rely on the properties of empirical cumulative distribution functions.

Scheffé’s theorem allows one to write the $L^{1}$-distance between two probability density functions in terms of the associated probability functions, as $\int_{\mathcal{Y}}|f_{1}(y)-f_{2}(y)|\text{d}y=2(\mathbb{P}_{1}(Y\in B)-\mathbb{P}_{2}(Y\in B))$, where $B$ is the set of values for which $f_{1}(y)>f_{2}(y)$. Since $B$ can be written as a union of intervals $(a(t),b(t))_{t=1}^{T}$, these probabilities can be calculated from the corresponding cumulative distribution functions. Thus, a cdf-based estimator of $\delta_{i}$ can be obtained as:

[TABLE]

For further details on the estimation of the intervals $(\hat{a}^{i}_{m}(t),\hat{b}^{i}_{m}(t))$, we refer to Plischke and Borgonovo (2017).

Since $\beta_{i}$ is itself a cdf-based sensitivity measure, the definition of a one-sample cdf-based estimator is straightforward:

[TABLE]

where $\hat{F}_{Y}$, and $\hat{F}^{i}_{m}$ are the empirical cdf’s of $\mathbf{y}$ and $\mathbf{y}^{i}_{m}$, respectively, i.e.:

[TABLE]

and $\mathds{1}_{A}(y)$ denotes the indicator function, taking the value $1$ if $y\in A$ and [math] otherwise.

Recalling that the expected value of a random variable $Y$ can be calculated as the integral of its survival function, $\mathbb{E}[Y]=\int_{\mathcal{Y}}(1-F_{Y}(y))\text{d}y$, a cdf-based one-sample estimator of the variance-based sensitivity measure, $\eta_{i}$ is given by:

[TABLE]

Notice that, since the empirical distribution functions are piece-wise constant, the integral in the above expression reduces to a sum. Plischke and Borgonovo (2017) propose an efficient way to calculate this integral.

Most of the estimators found in the literature, including those mentioned above, are constructed either as deterministic approximations or as (frequentist) point estimators. Therefore, quantification of the estimation error (or interval estimation) requires additional manipulation. Finding asymptotic distributions of the estimators in order to provide approximate confidence intervals is not straight forward, except for the variance-based estimator $\eta_{i}^{\star}$, and even in this case, they are accurate only for high sample sizes. For instance, Gamboa et al. (2016); Janon et al. (2014a); Tissot and Prieur (2015) show that variance-based estimators, calculated with a pick-and-freeze design or a replicated Latin hypercube, are asymptotically normal, but similar results are not available for other sensitivity measures. An alternative for non-deterministic sampling methods, is to replicate the estimation procedure in order to obtain a sample of estimates and corresponding sample-based confidence intervals. This, however requires a number $C>n$ of simulator evaluations and, as mentioned in the introduction, the computational cost of such an effort could be prohibitive for time-demanding realistic applications. The idea of replicates also excludes the use of quasi-random generators to create the probabilistic sensitivity analysis samples, as they are deterministic in nature. As a further alternative, bootstrap confidence intervals have been proposed in the literature in order to avoid the need for additional simulator runs (Plischke et al., 2013; Janon et al., 2014b), but these are not an integral part of the estimation process.

A second issue to consider when using partition-based one-sample methods is the sample size bias induced by the partition. Quantities related to the marginal distribution of $Y$ are estimated using the full sample size $n$, but those related to the within bin distribution of $Y|X^{i}\in\mathcal{X}_{m}^{i}$ are estimated using a smaller sample size $n^{i}_{m}\approx n/M$. While a sample size correction is implicit in the estimation of variances (see Eq. 3), the same is not true for the pdf and cdf estimates of equations (4) to (5). In other words, there is a different granularity when estimating the conditional and the unconditional distributions. In Section 3, we propose two partition-dependent Bayesian estimators which mitigate the sample size bias, while providing a natural way to quantify the estimation error, allowing interval estimation.

Within the Bayesian paradigm, unknown objects are treated as random, and assigned a prior probability measure which reflects the analyst’s uncertainty about their values. In this context, Oakley and O’Hagan (2004) treat the input-output mapping $g$ as unknown (at least before evaluation). Thus, they define a semi-parametric regression model with a Gaussian process prior, which allows posterior inference on variance-based sensitivity measures. In fact, it is possible to calculate posterior means for the conditional and unconditional variance of $Y$ and $Y|X^{i}$ respectively, either analytically or via numerical integration. The approach eliminates the need for a partition of the covariate space, thus solving the second issue mentioned above. However, the posterior distributions of the variance-based measures (e.g. $\eta_{i}$) are not available analytically, and finding a posteriori credibility intervals for estimation error quantification would be cumbersome and this aspect is not treated in the paper. Furthermore, it is not clear how to extend the results to the estimation of other (pdf or cdf-based) sensitivity measures. In Section 4, we present two alternative partition-free Bayesian models which allow interval estimation for these types of sensitivity measures as well. For illustrative purposes, we focus on estimation of the three measures in Table 1.

3 Bayesian non-parametric partition-dependent estimation

We propose to quantify the uncertainty about fixed but unknown sensitivity measures, $\xi_{i}$, before (a priori) and after (a posteriori) the observation of a sample, $\{(\mathbf{x}_{j},y_{j}):j=1,\ldots n\}$, within the Bayesian paradigm. The $\xi_{i}$ play the role of parameters of interest and they are linked to the data through functionals of the marginal and conditional distributions of $Y$ and $Y|X^{i}$. In view of this, it seems sensible to induce a prior on $\xi_{i}$ by assigning a prior to the family $\mathcal{P}_{i}=\{\mathbb{P}_{Y|X^{i}=x^{i}}:x^{i}\in\mathcal{X}^{i}\}$ of conditional probability measures. Notice that, since $\mathbb{P}_{X^{i}}$ is assumed known, the marginal distribution of $Y$, $\mathbb{P}_{Y}(y)=\int_{\mathcal{X}^{i}}\mathbb{P}_{Y|X^{i}=x^{i}}(y|x^{i})\text{d}\mathbb{P}_{X^{i}}(x^{i})$, is fully determined by $\mathbb{P}_{Y|X^{i}}$, so no additional prior specification is required. For each $i$, $\mathcal{P}_{i}$ is a family of probability measures on $\mathcal{Y}$, indexed by $x^{i}\in\mathcal{X}^{i}$, so defining a prior probability on this space is, in principle, not a simple task. Furthermore, the relation between $\xi_{i}$ and $\mathcal{P}_{i}$ is complex, making it difficult to conceive an adequate parametric prior. In other words, choosing a family of distributions characterized by a finite-dimensional parameter $\theta$, to express an expert’s uncertainty about $\xi_{i}$ through some prior on $\theta$ would seem overly restrictive, if not unreasonable. It is known that an inadequate prior may lead to troublesome posterior (Freedman, 1965) and hinder the properties of the proposed estimators. A natural alternative is to use a Bayesian non-parametric prior in order to ensure enough flexibility to capture complex data structures. Bayesian non-parametric methods are not restricted to a finite number of parameters to represent a distribution. Generally speaking, they rely on measure-valued stochastic processes to define priors on the space of probability measures of interest. The supports of such priors are wide, ideally covering the full range of all possible distributions, in our case, over $\mathcal{Y}$ (see e.g. Hjort et al., 2010, for an extensive discussion on bayesian non-parametric priors, their properties and their use).

Our first proposal can be interpreted as a Bayesian refinement of the cdf-based estimators introduced in the previous section and, as such, relies on a partition of the input space. We assume that the distribution of $Y|X^{i}=x^{i}$ is identical for every $x^{i}\in\mathcal{X}^{i}_{m}$, and denote it by $\mathbb{P}^{i}_{m}$. In practice, it is enough to assume that $\mathbb{P}_{Y|X^{i}=x_{i}}$ can be well approximated in this way. Prior uncertainty is expressed through a prior on the collection $\{\mathbb{P}^{i}_{m}\}_{m=1}^{M}$. For simplicity, we assume that such distributions are independent and identically distributed (i.i.d.), so the problem becomes that of finding a prior which assigns probability $1$ to a large enough set of probability distributions supported on $\mathcal{Y}$. We focus our attention on the Dirichlet Process (DP), first introduced by Ferguson (1973) and widely studied in the BNP literature (see e.g. Hjort et al., 2010, Chapter 2, for a discussion on its properties). We therefore define, for each $i=1,\ldots k$ the following Bayesian non-parametric model:

[TABLE]

where $\mathcal{DP}(\alpha G)$ denotes a Dirichlet process with base measure $G$ and concentration parameter $\alpha$. The Dirichlet process could be replaced by a more general stick-breaking process, achieving greater flexibility at a similar computational cost (see e.g. Ishwaran and James, 2001; Pitman and Yor, 1997; Lijoi et al., 2007). In this case, the algorithms and proposed estimators would maintain a similar structure so we focus on the Dirichlet process, without loss of generality, in order to use a notation more familiar to a wider audience. With regards to the hypothesis of independence between the $\mathbb{P}^{i}_{m}$, it could be removed through the application of recent developments in BNP methods (see Wood et al. (2011); Teh et al. (2006); Teh and Jordan (2010); Camerlenghi et al. (2017) and Camerlenghi et al. (2019)). This, however, would lead to a complication of the estimation algorithms which goes beyond the scope of this paper.

Note that this Bayesian model is coherent, in the sense that it induces a unique prior over the unconditional distribution of $Y$, whenever the partitions are equiprobable, that is when $\mathbb{P}(X^{i}\in\mathcal{X}^{i}_{m})=\frac{1}{M}$ for all $i=1,2,...,k$ and $m=1,2,...,M$. In fact,

[TABLE]

Then, by marginalizing, we obtain

[TABLE]

because $\int\mathbb{P}_{m}^{i}(\cdot)$d$\mathcal{DP}(\mathbb{P}_{m}^{i}|\alpha G)$ does not depend on $i$ or $m$. In other words, a priori, $\mathbb{P}_{Y}\overset{}{\sim}\mathcal{DP}(\alpha G)$, so that the prior for the marginal simulator distribution is also a Dirichlet process. This statement alone, however, provides no information on the probabilistic dependence of $Y$ on $X^{i}$. Thus, it is not meaningful, by itself, for a sensitivity analysis.

The posterior for this model, given the simulator input-output realizations ($Data$ for short), $\{(\mathbf{x}_{1},y_{1}),\ldots,(\mathbf{x}_{n},y_{n})\}$, can be written as follows:

[TABLE]

where

[TABLE]

Note that the posterior of the marginal for $Y$ can be obtained as:

[TABLE]

which may depend both on $i$ and $m$. However, the marginal coherence of the model still holds, at least asymptotically. Informally, for an equiprobable partition, $\mathbb{P}(X^{i}\in\mathcal{X}^{i}_{m})=1/M$, $n_{m}^{i}\simeq n/M$ when the sample size $n$ is sufficiently large, so $\alpha/(\alpha+n_{m}^{i})\simeq M\alpha/(M\alpha+n)$ and $n_{m}^{i}/(\alpha+n_{m}^{i})\simeq n/(M\alpha+n)$. Furthermore, $\sum(1/n_{m}^{i})\delta^{Dirac}(y)\simeq\sum(M/n)\delta^{Dirac}(y)$. Thus, asymptotically, $\mathbb{P}_{Y}(\cdot|\alpha G,Data)$ does not depend on $m$ or $i$ and $\mathbb{P}_{Y}(\cdot|\alpha G,Data)\overset{}{\sim}\mathcal{DP}((\alpha+n)\widetilde{G})$, where

[TABLE]

and $\widehat{\mathbb{P}}_{n}$ denotes the empirical distribution of $Y$ based on the full set of observations, $(y_{1},\ldots,y_{n})$. Note that this is the usual posterior corresponding to the DP prior on $\mathbb{P}_{Y}$.

The sensitivity measures we aim to estimate are functionals of the conditional and marginal distributions. The posterior means in eqs. (11) and (13), respectively, may be proposed as Bayesian estimators of such densities. Thus, a Bayesian point estimator of $\xi_{i}$ may be given by:

[TABLE]

Unfortunately, the direct calculation of $\widetilde{\xi}_{i}$ is impractical. Moreover, our purpose is to provide interval estimation, so as to quantify the uncertainty associated to point estimates. A way out is to sample observations (i.e., predicted realizations of the output) from $\widetilde{G}$ and $\widetilde{G}_{m}^{i}$, in order to enrich the sample. More specifically, we have a vector $\mathbf{y}$ of $n$ observations from the original simulator used to estimate $\mathbb{P}_{Y}$, but only $n_{m}^{i}$ of these belong to $\mathbf{y}^{i}_{m}$ and are therefore used to estimate $\mathbb{P}_{m}^{i}$. Because $n_{m}^{i}<n$, the precision issue discussed in Section 2.1 emerges, causing a bias in the empirical estimation of $\xi_{i}$. By re-sampling from $\widetilde{G}$ and $\widetilde{G}_{m}^{i}$ we can enlarge both vectors, making them of the same size and, potentially, arbitrarily large. Our proposal here is simply to sample $n-n_{m}^{i}$ observations from $\widetilde{G}_{m}^{i}$, thus obtaining two vectors of size $n$. The intuition underlying this corresponds to the non-parametric Bayesian bootstrap (Bb) (Hjort, 1985, 1991). In our case, for each $m$ a sample $\widetilde{\mathbf{y}}_{m}^{i}=\{\widetilde{y}_{n_{m}^{i}+1}^{i},\dots,\widetilde{y}_{n}^{i}\}$ of size $n-n_{m}^{i}$ is obtained from $\widetilde{G}_{m}^{i}$. A value of $\widehat{\xi}^{Bb,s}_{i}$ in Eq. (2) can be calculated through any of the methods discussed in Section 2.1, using $\mathbf{y}$ to estimate all quantities related to the marginal distribution of $Y$ and the extended vector $\mathbf{y}_{m}^{Bb,i,s}=(\mathbf{y}_{m}^{i},\widetilde{\mathbf{y}}_{m}^{i,s})$ to estimate all quantities related to the conditional distribution of $Y|X^{i}\in\mathcal{X}_{m}^{i}$. Informally, the weighted average over $m$ can be seen as approximately simulated from the posterior distribution of $\xi_{i}$. By repeating this procedure $S$ times, we obtain a Bb sample $\{\widehat{\xi}_{i}^{Bb,s}:s=1,2,...,S\}$. We propose the Monte Carlo average:

[TABLE]

as a point estimator of $\xi_{i}$. Approximate credibility intervals can be obtained from the empirical quantiles. Note that, because each $\widetilde{y}_{j}^{i}$ is simulated from a single distribution, $\widetilde{G}_{m}^{i}$, the sampling process can be done in parallel and the method is computationally fast. However, the uncertainty is underestimated because the additional variability captured by the posterior distribution of Eq. (10) is ignored.

A more accurate alternative is to sample $\widetilde{\mathbf{y}}_{m}^{i}$ jointly from the Dirichlet process posterior distribution (10), instead of sampling each $\widetilde{y}_{j}^{i}$ from the posterior mean. This can be done via the Pólya Urn scheme (Pu) of Blackwell and MacQueen (1973). Specifically, $\widetilde{\mathbf{y}}_{m}^{i}$ is generated as a realization of the Pólya sequence:

[TABLE]

Once again, the extended samples $\mathbf{y}_{m}^{Pu,i,s}=(\mathbf{y}_{m}^{i},\widetilde{\mathbf{y}}_{m}^{i,s})$ can be used to obtain a value $\widehat{\xi}^{Pu,s}_{i}$ by any available method to calculate the expression in Eq. (2). We use $\widehat{\xi}^{Pu}_{i}$ to denote the Monte Carlo average of a sample of size $S$ generated in this way. Note that this is a point estimator with the same expectation as $\widehat{\xi}^{Bb}_{i}$. However, a greater variability which fully accounts for the uncertainty on $\mathbb{P}_{m}^{i}$ results in wider credibility intervals. The sampling procedure is now sequential for $s=1,\ldots,S$, so the price for greater accuracy in uncertainty estimation is a slightly higher computational time.

The technical details for Bb and Pu estimators are presented in Section A.2.1 (Appendix A.2).

3.1 Simulation study

We illustrate the performance of the two classes of estimators proposed above, via two toy examples for which the sensitivity measures can be calculated analytically (see Table 2 columns $1$ and $2$ ). The first example is the 2-input simulator

[TABLE]

where $X^{1}$,$X^{2}\overset{iid}{\sim}\text{Gamma}(3,1)$, so that the output $Y$ follows a Beta distribution. The second example is the 21-input simulator

[TABLE]

where $X^{i}\sim\text{Normal}(1,1)$, with $a_{1}=\dots=a_{7}=-4$, $a_{8}=\dots=a_{14}=2$, and $a_{15}=\dots=a_{21}=1$. The 21 inputs are correlated with $\rho_{i,j}=0.5$. Therefore, Model inputs with indices in $1,2\dots 7$ are the most important, followed by inputs with indices in $15,16,\dots,21$, followed by inputs with indices in $8,9\dots 14$. Columns $3$ to $6$ in Table 2 displays the analytical values of the sensitivity measures.

We are interested in small sample sizes, which make the estimation of global sensitivity measures challenging. In particular, we consider $n=\{300,600,900\}$. The input data, $\mathbf{x}$, is generated via Quasi-Monte Carlo. For each $n$, alternative choices of the partition size, $M$, are explored. The mass parameter, $\alpha$, for the DP prior is set equal to $0.1n/M$ throughout. The base measure, $G$, is chosen in correspondence with the support of $Y$: a Beta distribution for the first example and a Normal distribution for the second; the hyper-parameters are fixed through an empirical approach, based on the available sample $\mathbf{y}$. Note that this choice centres the prior distribution for $Y|X^{i}\in\mathcal{X}_{m}^{i}$ roughly around the marginal distribution of $Y$, thus favouring, a priori, independence between the $Y$ and $X^{i}$, with a precision proportional to the number of observations in each partition set. In practical applications, prior information elicited from experts may be expressed through different choices of $\alpha$ and $G$.

We compare the Bayesian bootstrap and Pólya urn estimators to traditional point estimators for three global sensitivity measures. Results are reported in Fig.s 1 and 2: the first row corresponds to $\eta_{i}$, the second to $\delta_{i}$ and the third to $\beta_{i}$. Columns, from left to right, correspond to increasing sample sizes. Each graph is divided into three blocks displaying importance measures estimates based on alternative choices of $M$. The dotted lines display the analytical values.

We first consider the left-most panel of Fig. 1(a). At $n=300$ the estimates vary notably with the partition size: they are downward biased for $M=3$ and upward biased for $M=21$. Observe that at $M=21$, we have $n_{m}^{i}\simeq 9$, a number too small to be reasonably chosen by the analyst. However, the bias is systematic, that is, it affects identically all estimates. Estimates are less affected by the partition choice as the sample size increases. Recall that in realistic applications, where an analyst would not know the true values of the sensitivity measures, the main interest is on the ordinal ranking of the inputs. In this example, $X^{1}$ and $X^{2}$ are equally important. However, looking at the point estimators $\widehat{\eta}_{i}^{\star}$ and $\widehat{\eta}_{i}^{\diamond}$ the analyst would rank $X^{2}$ as more important than $X^{1}$ for most combinations of $n$ and $M$. The credibility intervals for $\widehat{\eta}_{i}^{\text{Pu}}$ display a large overlapping that would prevent the analyst from ranking $X^{2}$ above $X^{1}$: there is too much uncertainty in the estimates to make such conclusion. Notice the underestimation of the uncertainty surrounding $\widehat{\eta}_{i}^{\text{Bb}}$. Rows (b) and (c) of Fig. 1 show a similar behaviour for the $\delta_{i}$ and $\beta_{i}$ sensitivity measures.

Fig. 2 shows the estimates for the $21-$input simulator in (16). For a better display clarity, instead of reporting seven sensitivity measures per group, we show numerical values for a representative of each input group, namely $X^{3},X^{10}$ and $X^{18}$. For variance-based sensitivity measures (Fig. 2(a) ) for all $n$ and $M$ considered, $\widehat{\eta}_{i}^{\star}$, $\widehat{\eta}_{i}^{\diamond}$, $\widehat{\eta}_{i}^{\text{Bb}}$ and $\widehat{\eta}_{i}^{\text{Pu}}$ are able to correctly identify $X^{3}$ as the most relevant variable. Regarding $\delta_{i}$ and $\beta_{i}$ (Fig. 2(b) and Fig. 2(c), respectively), one identifies $X_{3}$ as the most important input in almost all combinations of sample sizes and partition selections. The ranking becomes unclear for $n=300$ and $M=21$. However, this choice would leave about $9$ realizations per partition, a number too small to be reasonably chosen by the analyst. For the remaining group of inputs, the overlapping error bands for the $Bb$ and $Pu$ estimates would not allow us to deem $X^{10}$ more relevant than $X^{18}$, with either $\eta_{i}$,$\delta_{i}$ or $\beta_{i}$. Thus higher sample sizes would be needed for neatly ranking the second and third most important groups of simulator inputs.

Overall, Fig.s 1 and 2 suggest that the proposed estimators allow the identification of the most important inputs, even at small sample sizes, and, most relevantly, they provide a measure of the uncertainty in the assessment. However, the results also display a strong dependence on the partition size $M$. While i) as observed in Strong and Oakley (2013) (see their Fig. 1, at p. 759), the importance of selecting an optimal partition size diminishes as the sample size increases and ii) a suboptimal partition selection has in most cases an identical impact on the sensitivity measures (i.e., the sensitivity measures of all inputs are simultaneously upward or downward biased), the analyst is still left with the question of what is the optimal partition size for a given sample. Unfortunately, there seems to be no universally optimal selection rule (see Appendix A.1). Clearly, the problem would be solved if partition-independent estimators were available. In the next section, we study two proposals of Bayesian estimators that avoid the partition choice problem.

4 Bayesian non-parametric partition-free estimation

In this section, we propose two classes of Bayesian partition-free estimators. The first is based on the use of an infinite mixture model to estimate the joint density of $Y$ and $X^{i}$. The second, uses a Bayesian non-parametric regression model to estimate the conditional density of $Y$ given $X^{i}$.

4.1 Joint density-based estimation

The intuition is that all sensitivity measures under consideration can be recovered from the joint distribution of $Y$ and $X^{i}$. Therefore, in order to do Bayesian inference on $\xi_{i}$ it suffices to place a prior on the joint density $f_{X^{i},Y}$. We propose to do so by means of a nonparametric mixture model (see, e.g. Ferguson (1983); Lo (1984)). In other words, we consider $f_{X^{i},Y}$ to be defined as a mixture:

[TABLE]

where $\mathcal{K}$ is a parametric bivariate density and the mixing measure $P$ is a probability distribution over an appropriate space of parameters. The model is completed by assigning a non-parametric prior, $\Pi$, on $P$. Most common choices of $\Pi$ assign probability one to discrete distributions of the form

[TABLE]

placing mass $w_{\ell}$ on locations $(\theta_{\ell})$. In the literature, particular attention has been paid to nonparametric priors admitting a stick-breaking construction (Pitman, 1996; Sethuraman, 1994) where the weights $\underline{w}=(w_{1},w_{2},...)$ are defined as realization of random variables satisfying

[TABLE]

and independent of $\underline{\theta}=(\theta_{1},\theta_{2},...)\stackrel{{\scriptstyle iid}}{{\sim}}G$. Rich families of stick-breaking priors can be defined via different distributional assignments for the sequence $(V_{1},V_{2},\ldots)$ (Favaro et al., 2012; Ishwaran and James, 2001, see, e.g.,). The main advantage over other types of construction is that the stick-breaking representation of the random weights allows for efficient simulation algorithms, specially in the context of nonparametric mixture models (Ishwaran and James, 2001; Papaspiliopoulos and Roberts, 2008; Kalli et al., 2011; Yau et al., 2011). However, the most popular stickbreaking prior remains the Dirichlet process, well known even outside the specialized community of Bayesian nonparametrics. For this reason, we will focus our analysis on DP mixtures, thus letting $P\sim\Pi=\mathcal{DP}(\alpha G)$. Additionally, for simplicity, we choose $\mathcal{K}$ to be a bivariate normal density, following the density estimation scheme of Escobar and West (1995). In this case, $\theta_{\ell}=(\mu_{\ell},\Sigma_{\ell})$ and, to simplify calculations, we select $G$ as a conjugate Normal inverse-Wishart distribution. Thus, the the integral in (17) reduces to a sum and the joint density can be written as:

[TABLE]

where the weights follow (19), with $V_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\text{Beta}(1,\alpha)$.

Inference on this model is usually achieved via an MCMC scheme resulting in a sample from the posterior distribution of $f_{X^{i},Y}$ given the $Data$. In the case of the DP-mixture, the function DPdensity from the R package DPpackage provides an off-the-rack solution. In practice, the MCMC scheme generates, at each iteration $s=1,\ldots,S$, values $(\underline{w}^{s},\underline{\mu}^{s},\underline{\Sigma}^{s})$ which, substituted in expression (20), produce a density function, $f^{BNJ,s}_{X^{i},Y}$. Analytical expressions for the marginal and conditional densities, $f^{BNJ,s}_{Y}$ and $f^{BNJ,s}_{Y|X^{i}}$ as mixtures of normal distributions are made easily available by the choice of the Gaussian kernel. Clearly, it is also possible to evaluate the corresponding cumulative distribution functions. Thus, it is possible compute the global sensitivity measures of interest, $\eta_{i}^{BNJ,s}$, $\delta_{i}^{BNJ,s}$, $\beta_{i}^{BNJ,s}$ from their definitions (Table 1), obtaining a posterior sample of each. We denote the sample means by $\widehat{\eta}^{BNJ}_{i}$, $\widehat{\delta}^{BNJ}_{i}$ and $\widehat{\beta}^{BNJ}_{i}$, respectively, proposing them as Bayesian point estimators. Approximate credibility intervals can be obtained from the empirical quantiles of the samples. The procedure is summarized in Section A.2.2 of Appendix A.2, to which we refer for further details.

It is important to observe that the known marginal distribution for $X$ does not, in general, coincide with the marginal distribution for $X$ derived from each $f^{BNJ,s}_{X^{i},Y}$. Thus, by using only the joint density $f_{X^{i},Y}$ to estimate the sensitivity measures, important information, standard in global sensitivity analysis is wasted. In fact, inference for conditional densities based in the joint model is known to be approximate (see e.g. Müller and Quintana, 2004). In the next section, we present an alternative estimation method which avoids this problem through a recent Bayesian approach to conditional density estimation.

4.2 Conditional density-based estimation

We now propose to use a Bayesian non-parametric regression model to do inference directly on the conditional density of $Y|X^{i}$, thus using all of the information contained in the $Data$ to estimate the relationship between the variables and exploiting the knowledge of the marginal distribution of $X$ to obtain the marginal distribution of $Y$. The idea is to transform the non-parametric mixture of equation (20) into a mixture of conditional densities:

[TABLE]

This time a non-parametric prior, $\Pi$, is placed on the family, $\{P_{x}\}_{x\in\mathcal{X}}$ of mixing distributions indexed by $x$. Analogous to the DP mixture model of the previous section, a dependent DP mixture model or DDP mixture (MacEachern, 1999, 2000) is obtained when $P_{x}$ follows a DP prior, marginally for every $x$, so that:

[TABLE]

The random covariate-dependent weights $W_{\ell}(x)$ follow the stick-breaking construction of Eq. (19), for i.i.d. random processes $\{V_{\ell}(x):x\in\mathcal{X}\}$. In other words, $\mathbf{V}(x)\sim\mathcal{DP}$ for every $x$. It has been proved sufficient flexibility is achieved through models in which only the particles $\theta_{\ell}$ or the weights $w_{\ell}$ depend on the covariate $x$ (Barrientos et al., 2012), the second option being favoured due to better predictive capabilities. Several proposals have been studied in the literature, focusing on alternative definitions of the random functional weights $w_{\ell}(x)$ (e.g. Dunson and Park, 2008; Griffin and Steel, 2006; Dunson and Rodriguez, 2011).

The stick-breaking structure of the weights, which imposes a geometric decay, may be bypassed through an alternative construction allowing further flexibility:

[TABLE]

The denominator of this expression is, again, an infinite mixture of parametric kernels, $\mathcal{K}$, this time with support $\mathcal{X}$. Each $\omega_{\ell}$ can be interpreted as the probability that a realization of $Y$ comes from the $\ell$-th regression component regardless of the value of $X$, just as $\omega_{\ell}$ is the conditional probability given $X=x$. Such density regression model, where the weights $w_{\ell}$ in (23) follow the stick-breaking representation of (19) and the extended parameters $(\theta_{\ell},\psi_{\ell})$ are i.i.d. from some adequate base measure, $G$, was proposed by Antoniano-Villalobos et al. (2014), to which we refer the reader for additional details on the role and choice of hyper parameters, as well as the algorithm used for inference.

We adopt this construction to estimate the conditional density $f_{Y|X^{i}}(y|x_{i})$ as a mixture of linear regression models:

[TABLE]

where $w_{\ell}(x_{i})$ is given by Eq. (23), with a DP prior. Once again, a MCMC approach is used to generate a sample, this time from the posterior distribution of $f_{Y|X^{i}}$. Each $f^{BNC,s}_{Y|X^{i}}(y|x_{i}),s=1,\ldots,S$, together with the known marginal for $X^{i}$ can be used to calculate (e.g. by numerical integration) a corresponding marginal for $Y$. As discussed in Section 4.1, this is all that is needed to compute the global sensitivity measures of interest, $\eta_{i}^{BNC,s}$, $\delta_{i}^{BNC,s}$ and $\beta_{i}^{BNC,s}$. These, again allow point estimation, e.g. via the Monte Carlo averages, which we denote by $\widehat{\eta}_{i}^{BNC}$, $\widehat{\delta}_{i}^{BNC}$ and $\widehat{\beta}_{i}^{BNC}$, and interval estimation, via empirical quantiles. Section A.2.3 (Appendix A.2) summarizes the procedure and offers additional technical details.

4.3 Simulation study

We examine the performance of the classes of partition-independent estimators proposed in Sections 4.1 and 4.2, via the $2-$input and $21-$input simulators introduced in Section 3.1. For both joint and conditional density-based estimation, we set a burn-in period as $10n$ and the stored MCMC samples size $S=1000$. Results are illustrated in Fig. 3.

Figure 3(ii) reports results for the $21-$input simulator. The Bayesian non-parametric joint estimators $\widehat{\eta}^{BNJ}_{i}$, $\widehat{\delta}^{BNJ}_{i}$ and $\widehat{\beta}^{BNJ}_{i}$ correctly recover the true values of the parameters and, as the sample size increases from $n=300$ to $n=600$, the credibility intervals become narrower. At $n=900$, there is no more overlap among the three groups of sensitivity measures, allowing the analyst to rank the inputs neatly. Regarding the Bayesian non-parametric conditional estimates, we observe that $X_{3}$ is correctly identified as the most relevant input at all sample sizes. The values $\eta_{10}^{BNC},\delta_{10}^{BNC}$, $\beta_{10}^{BNC}$ as well $\eta_{18}^{BNC},\delta_{18}^{BNC}$, $\beta_{18}^{BNC}$ are overestimated by the $BNC$ estimators at $n=300$. However, the bias is reduced as $n$ increases and at $n=900$ the credibility intervals appear centred around the analytical values of the sensitivity measures. We also observe that for both classes of estimates the analytical value of the sensitivity measures falls within the $95\%$ credible intervals.

For this example joint Bayesian estimators seem to outperform their conditional counterpart. This is to be expected, because the joint Gaussian structure of the data is more easily recovered by the joint model in this case, so the loss due to ignoring the true distribution of $X^{i}$ has a lesser effect on the results. However, we can appreciate a reassuring improvement of the BNC estimates as the sample size increases. One may argue that, in a situation in which the true conditional distribution of $Y$ given $X^{i}$ is unknown and may be complex, estimation based on the conditional density model may be preferred, as more robust; the price to pay is that a larger sample size may be required, specially in high-dimensional situations. We then challenge these results for $2-$input simulator in Eq. (15), in which the distributions are not normal.

Assume for the moment that the analyst does not know the true value of the sensitivity measures. In terms of ordinal ranking, Fig. 3(i) suggests that the two simulator inputs are equally important. The credibility intervals of $X_{1}$ and $X_{2}$ obtained with both the BNJ and BNC estimators are overlapping at the all sample sizes and for all sensitivity measures, so that the analyst cannot deem one of them more important than the other. The performance of the two estimators is similar. However, note that the credible intervals of the joint model (BNJ) are wider than those for the conditional model especially for variance-based sensitivity measures. As expected, for a non-normal distribution, the joint model resents from the wrongful introduction of the marginal distribution for $X_{i}$. We analyze this behavior further in addressing results for the LevelE case study.

5 Case study: LevelE simulator

In this section, we evaluate the performance of the proposed estimators through the benchmark simulator of sensitivity analysis, LevelE. The LevelE code simulates the release of radiological dose from a nuclear waste disposal site to humans over geological eras. The code has been developed in an international exercise launched by the Nuclear Energy Agency (NEA) in the mid 1980’s Nuclear Energy Agency (1989). Goal of the exercise was the realization of a reference simulator for the prediction of flow and transport of radionuclides in actual geologic formations against which to compare other simulators developed internationally to support the selection of radioactive wast management policies. Since then, LevelE has become the benchmark simulator of sensitivity analysis (Saltelli et al., 2000; Saltelli and Tarantola, 2002). During the international exercise, distributions for the uncertain simulator inputs were assessed (Table 3), and have become the reference for analysis on this code. From a technical viewpoint, the LevelE code solves of a set of nested partial differential equations that compute the released radiological dose in Sievert/year over a time range of $t=10,000$ to $2\times 10^{9}$ years. The detailed equations of the code are reported in Saltelli and Tarantola (2002).

Previous works have discussed the sensitivity analysis of this simulator using alternative sampling methods and sizes. For instance, Saltelli et al. (2000) employ $3,084$ simulator evaluations to obtain point estimates of the first and total order variance-based sensitivity indices. Saltelli and Tarantola (2002) employ $10,000$ simulator runs for the point estimation of first-order variance-based sensitivity indices, a second experiment with $16,384$ runs for the point of the first and total order sensitivity indices according to the design in Saltelli (2002) (no uncertainty in the estimates is provided). In Ratto et al. (2007), stable patterns for the estimation of variance-based sensitivity measures are obtained at a cost of about $1,024$, after the input-output dataset has been used to train an emulator. In Castaings et al. (2012), design based on substituted columns sampling and permuted columns sampling are used, with convergence at about $10^{4}$ runs. Wei et al. (2014) propose a copula-based estimation methods that reduces the cost to about $1,000$ runs for point estimates, with $20$ replicates for obtaining confidence intervals. Plischke and Borgonovo (2017) apply a given-data design for the point estimators $\widehat{\eta}^{\diamond}_{i}$, $\widehat{\delta}^{\diamond}_{i}$ and $\widehat{\beta}^{\diamond}_{i}$ using a sample up to size $n=5,000$, with estimates becoming stable for $n>1,000$ runs. Thus, a sample of size $n=1,000$ can be considered reflective of state of art for the identification of the key-uncertainty drivers of LevelE.

We report results for the calculation of global sensitivity measures using all classes of estimators discussed in the present work for samples of sizes $n=600$ and $n=900$. Figures 4 and 5 display the results.

The graphs in Fig. 4 report the Bayesian bootstrap and Pólya urn estimators, vis-á-vis the point estimators for variance-based (graphs in row a), density-based (graphs in row b) and cdf-based (graphs in row c) sensitivity measures. The results show that already at $n=600$ the two most important simulator inputs are correctly identified. However, the estimates are sensitive to the partition size. Consider the right graph in row a). The credibility intervals of the variance-based Pólya urn estimators with $M=26$ are completely overlapping. This signals that, had the analyst chosen such partition size, the estimates would not be meaningful. The separation becomes, instead, clearer at smaller partition sizes with $M=9$ being possibly the optimal choice. Note that the estimates tend to be upward biased as the partition size increases, in agreement with our previous experiments and also with previous literature findings.

We then come to the joint and conditional partition-independent Bayesian density estimators (Fig. 5).

The two graphs in row (a) display the estimates and credibility intervals for variance-based sensitivity measures ($\widehat{\eta}^{BNJ}_{i}$,$\widehat{\eta}^{BNC}_{i}$), the two graphs in row (b) for density-based sensitivity measures ($\widehat{\delta}^{BNJ}_{i}$, $\widehat{\delta}^{BNC}_{i}$) and the two graphs in row (c) for cdf-based ($\widehat{\beta}^{BNJ}_{i}$, $\widehat{\beta}^{BNC}_{i}$) sensitivity measures. Figure 5 shows that the two key-uncertainty drivers are correctly identified already at $n=600$, by $\widehat{\eta}^{BNC}_{i}$, $\widehat{\delta}^{BNC}_{i}$ and $\widehat{\beta}^{BNC}_{i}$, as the credibility intervals of the associated sensitivity measures separate from the credibility intervals of the remaining simulator inputs. The $\widehat{\delta}^{BNJ}_{i}$ and $\widehat{\beta}^{BNJ}_{i}$ correctly identify the two most influential simulator inputs. However, $\widehat{\eta}^{BNJ}_{i}$ fails to produce meaningful results for variance-based sensitivity measures at either sample sizes. This confirms the results of Section 4.3. The deviation from normality strongly affects the ability of BNJ to capture the conditional density of $Y$ given $X_{i}$, since much of the information contained in the data goes into the estimation of unnecessary components of the density mixture of the marginal distribution of $X_{i}$. This reduces the estimation precision, leading to the wider confidence intervals.

Let us consider the perspective of an analyst interpreting the results overall. From the available Data, the analyst is able to obtain alternative estimators for representatives of three categories of sensitivity measures, with display of credibility intervals. With the exception of $\widehat{\eta}^{BNJ}_{i}$, the estimators communicate that uncertainty in the simulator response is mostly driven by two simulator inputs, with the remaining ones being of lower significance. Thus, the analyst is allowed to confidently report the key-uncertainty drivers to the decision-maker even if the sample size is limited. At the same time, Fig.s 4 and 5 communicate that the sample is not sufficient to rank the medium and low-important simulator inputs with confidence. If the decision-maker (modeler) wished sharper estimates of the sensitivity measures of these inputs, the analyst would need a larger sample size. This could be obtained either through additional runs of the original simulator or by fitting an emulator and, in case the fit is accurate, running the emulator instead of the original code.

6 Conclusions

This work has presented a fully Bayesian approach to the estimation of probabilistic sensitivity measures from a given sample. The proposed algorithms yield credibility intervals for the estimates without increasing computational burden. We have studied four classes of estimators. The first two find their theoretical ground in non-parametric Bayesian estimation based on the Dirichlet process. These estimators run in parallel with one-sample frequentist estimators currently in use, produce uncertainty in the estimates and are computationally simple to implement. However, they leave the analyst with the problem of choosing the optimal partition. The introduced conditional and unconditional non-parametric Bayesian estimators eliminate the partition selection problem, while producing uncertainty in the estimates. However, their numerical implementation needs to be carefully executed, as it requires a combination of numerical integration and MCMC. Algorithms are available, but their convergence might take a longer time than the Bayesian bootstrap and Pólya urn estimators. Then, how should one proceed in a practical situation? The several numerical experiments performed by the authors (of which a subset was reported in the paper) evidence that the estimators succeed in identifying key-uncertainty drivers at small sample sizes in most situations. Then, a suggested approach would be to apply first the Bayesian bootstrap and/or Pólya urn estimators on the available sample for computing an ensemble of sensitivity measures (e.g., $\eta,\delta,\beta$). If the sensitivity measure estimates and credibility intervals yield a clear picture of the simulator inputs influence, then the analysis could be considered satisfactory. However, the analyst ought to test this assertion repeating the estimates at alternative partition sizes. In case results are strongly dependent on the partition size, the analyst can invest in the Bayesian non-parametric estimation. If these estimators yield a clear picture about the simulator input influence, the analysis is conclusive. Conversely, a larger sample is needed and the analyst ought to plan for additional simulator runs.

While we have discussed three well-known global sensitivity measures, the paradigm presented here can be applied to the estimation of any global sensitivity measure, including, among others, value of information, sensitivity measures based on any discrepancy between densities or cumulative distribution functions.

From a more general perspective, the work shows that combining recent advances in Bayesian non-parametric density estimation with probabilistic sensitivity analysis in DACE may lead to improvements in the estimation of global sensitivity measures. Research in Bayesian non-parametric density estimation is active in Statistics and Machine Learning, but the advances in this discipline are not directly known to the DACE community. This work represents a first systematic bridge between these two closely related areas of Statistics, and we hope it could favour further research for transferring findings in Bayesian-non parametric estimation to the field of computer experiments. At the same time, exposing Bayesian estimation to the demands coming from probabilistic sensitivity analysis of realistic simulators may challenge state of the art and stimulate further research in Bayesian-non parametric estimation.

Appendix A Appendix

A.1 Numerical experiments for the partition selection problem

The authors performed several thought experiments on test cases. The results show the difficulty, maybe impossibility, of finding a universally valid rule for linking the partition size $M$ to the sample size $n$. We report some experiments results.

Assume the analyst wants to find an “optimal ”(in some sense) partition refining strategy, i.e., a relationship that produces the partition size $M$ that minimizes the estimation error at sample size $n$ for the pdf-based point estimators $\widehat{\eta}_{i}^{\star}$, $\widehat{\delta}_{i}^{\star}$ and cdf-based point estimator $\widehat{\beta}_{i}^{\diamond}$ [eqs. (3), (4) and (6)]. We focus on one estimator type for simplicity and also because Borgonovo et al. (2016) propose an heuristic inspired by the rule of histogram partitioning of Freedman-Diaconis (Freedman and Diaconis, 1981), in which $M\sim\sqrt[3]{n}$.

To evaluate the estimators’ performance at fixed values of $M$ and $n$, we use the Root Mean Square Error (RMSE):

[TABLE]

where $S$ is the number of bootstrap replicates. $\hat{\xi}_{i,l}$ is the $l$-th bootstrap replicates of $\xi_{i}$.

We estimate the sensitivity measures with sample sizes varying from $300$ to $900$, and partition sizes covering the natural numbers between $5$ and $35$. Then we calculate the RMSEs with $S=100$ bootstrap replicates. Figures 7 and 7 present the heatplot of RMSEs in percentage (RMSE${}_{i}/\xi_{i}\cdot 100\%$). The horizontal axis indicates the sample size, and the vertical axis the partition size. The darker the color of a region in the plot, the lower the estimation error. For example, in Fig. 7(a), dark (blue) refers to low RMSE (less than 10 percent), and light (red) to relative high RMSE (higher than 14 percent). The magenta line maps n into M using the previously mentioned heuristic function. Figure 7 shows that the proposed heuristic works well on the 2-input simulator (Eq. (15)), with the magenta line falling mainly into dark coloured regions. However, for the 21-input simulator (Eq. (16)) we would incur in high errors at small sample sizes. For instance consider graph a) in Fig. 7. The graph reports the error in the estimates of $\delta_{3}$ for the second model. The heuristic would propose values of $M$ at about $20$ for all values of $n$ as optimal partition sizes. However, the partition size that minimizes the error is at about $M=10$ or lower. The different behavior here could also be related to the differences in structure and dimensionality of the models. However, even for the same model, the heatplots differ significantly across the sensitivity measures. For the first simulator (Eq. (15)), the ideal partition size for the variance-based estimator is between 10 to 15 (Fig. 7 (b)), while for mutual information, if falls between 20 to 25 (Fig. 7 (c)).

These results show that aiming at postulating a universally valid heuristic might be a cumbersome task.

A.2 Implementation details for the Bayesian non-parametric estimators

We present further details regarding the implementation of the Bayesian non-parametric estimation methods in Sections 3 and 4. Inference on the three selected sensitivity measures $\eta_{i}$, $\beta_{i}$ and $\delta_{i}$ is performed independently for each $i=1,\ldots,k$. Therefore, in order to simplify the notation, we will leave out the index $i$ throughout this appendix, considering its value fixed. Throughout this section, all the integrals are approximated numerically using trapezoidal rule, and all the supremes are approximated by the maximum on a predetermined grid over $\mathcal{Y}$.

A.2.1 Partition-dependent bootstrap and Pólya urn estimation

Recall that in Sections 3, given $M$, we have the partition $\{\mathcal{X}_{m}\}^{M}_{m=1}$ of $\mathcal{X}$ according to the sample proportion and corresponding $\{\mathbf{y}_{m}\}$.

Within each partition set, we generate $n-n_{m}$ new points $\widetilde{\mathbf{y}}_{m}^{s}$ and obtain the extended vector $\mathbf{y}_{m}^{C,s}=(\mathbf{y}_{m},\widetilde{\mathbf{y}}_{m}^{s})$ with $C\in\{Bb,Pu\}$, where $\widetilde{\mathbf{y}}_{m}^{s}$ is sampled from the posterior mean $\widetilde{G}_{m}$ for $C=Bb$, and is generated through Pólya urn scheme when $C=Pu$. The superscript $s$ is used to indicate the $s-$th replicate.

After repeating the sampling procedure for $S$ times, we obtain the partition-depended Bayesian estimator of $\eta$ by calculating the Monte Carlo average:

[TABLE]

where $\bar{y}_{m}^{C,s}$ is the sample mean of $\mathbf{y}_{m}^{C,s}$; $\bar{y}$ and $s_{y}^{2}$ are the sample mean and variance of $\mathbf{y}$. Approximate credibility intervals of $\eta$ can be obtained from the empirical quantiles of $\{\eta^{C,s},s=1\dots,S\}$. The same intuition is used for $\delta$ and $\beta$. Specifically, we use

[TABLE]

where $\hat{f}^{\star}_{Y}$ and $\hat{f}_{m}^{C,s}$ are kernel smoothing functions of $\mathbf{y}$ and $\mathbf{y}_{m}^{C,s}$, respectively; $\hat{F}_{Y}$, and $\hat{F}_{m}^{C,s}$ are the empirical cdf’s of $\mathbf{y}$ and $\mathbf{y}_{m}^{C,s}$, respectively.

Note that the calculations of $\eta^{C,s}$, $\delta^{C,s}$ and $\beta^{C,s}$ are equivalent to the pdf-based estimators in eqs. (3), (4), (6) but with the enriched samples. Alternatively, the cdf-based estimators in eqs. (8) and (5) could be used for $\eta^{C,s}$ and $\delta^{C,s}$.

A.2.2 Partition-free joint density-based estimation

Following the proposal in Jara et al. (2011), we fix $\alpha=1$, and choose $G$ to be a Normal-Inverse Wishart distribution

[TABLE]

where $\mathcal{N}(\cdot|m,A)$ denotes a bivariate normal distribution with mean $m$ and covariance matrix $A$, and $IW(\cdot|4,\psi)$ denotes an Inverse-Wishart distribution with mean $\psi^{-1}$. A hyper-prior is assigned to the parameters of the base measure, with hyperparameters determined empirically:

[TABLE]

where $\text{Gamma}(\cdot|a_{1},a_{2})$ denotes the Gamma distribution with mean $a_{1}/a_{2}$; $m_{2}=(\mu_{X},\bar{y})$ and $s_{2}=\text{diag}(\sigma_{X}^{2},s_{y}^{2})$.

Inference is achieved through the function DPdensity from the DPpackage in R. The output is a MCMC posterior sample $\underline{\theta}^{s}=(\underline{w}^{s},\underline{\mu}^{s},\underline{\Sigma}^{s})$, $s=1,\ldots,S$. In practice, the number $J_{s}$ of components with non-zero weights is finite, thus we have

[TABLE]

Given the posterior realizations, the corresponding joint density can be obtained:

[TABLE]

By the properties of the bivariate Normal distribution, the marginal and conditional distributions, $f^{BNJ,s}_{Y}$ and $f^{BNJ,s}_{Y|X^{i}}$ recpectively, are also mixtures of Normal distributions:

[TABLE]

where $\nu_{\ell}^{s}=\mu_{2,\ell}^{s}+\sigma_{3,\ell}^{s}(x-\mu_{1,\ell}^{s})/\sigma_{1,\ell}^{s}$ and $\tau_{\ell}^{s}=\sigma_{2,\ell}^{s}-(\sigma_{3,\ell}^{s})^{2}/\sigma_{1,\ell}^{s}$. Clearly, the corresponding cdfs, $F^{BNJ,s}_{Y}$ and $F^{BNJ,s}_{Y|X}$, as well as the marginal mean and variance can be calculated trivially. In particular,

[TABLE]

Thus, MCMC samples of the sensitivity measures of interest can be obtained as follows:

[TABLE]

where

[TABLE]

[TABLE]

Point estimators of interest are obtained as Monte Carlo averages:

[TABLE]

A.2.3 Partition-free conditional density-based estimation

Following the proposal of Antoniano-Villalobos et al. (2014), we fix $\alpha=1$ and choose $\mathcal{K}(x|\psi_{\ell})$ to be a Normal kernel, with $\psi_{\ell}=(\mu_{\ell},\tau)$. The base measure $G$ is given by:

[TABLE]

where $\mathbf{b}_{\ell}=(a_{\ell},b_{\ell})$.The hyperparameters are chosen empirically. As an illustration, consider the 21-input simulator. Figure 8 shows the scatter-plot of $(\mathbf{x}^{3},\mathbf{y})$ and the convex hull, i.e. the smallest convex set containing all points. In this case, we fix $\mathbf{b}_{0}=(-1.5,-5.5)$ and $C^{-1}=\text{diag}(43^{2},11^{2})$, in order to allow each local linear component to lie between the blue and red lines in the figure, which represent the main behaviour of the data.

We use the Matlab subroutine provided by Antoniano-Villalobos et al. (2014) to generate an MCMC posterior sample $(\underline{\theta}^{s},\underline{\psi}^{s})=(\underline{a}^{s},\underline{b}^{s},\underline{\sigma}^{s},\underline{\omega}^{s},\underline{\mu}^{s},\tau^{s})$, $s=1\ldots S$, where

[TABLE]

Given the a posterior realization $(\underline{\theta}^{s},\underline{\psi}^{s}))$, a conditional density can be obtained from eqs. (23) and (24):

[TABLE]

The corresponding marginal pdf $f^{BNC,s}_{Y}$ of $Y$ is obtained by integrating with respect to the true $f_{X}$:

[TABLE]

Clearly, the corresponding marginal and conditional cdfs, $F^{BNC,s}_{Y|X}$ and $F^{BNC,s}_{Y}$, respectively can be obtained trivially. In particular, posterior realizations of the marginal mean and variance of $Y$ are given by

[TABLE]

Thus, MCMC samples of the sensitivity measures of interest can be obtained as follows:

[TABLE]

where

[TABLE]

[TABLE]

Finally, point estimators of interest are obtained as Monte Carlo averages:

[TABLE]

**Note ** The code can be downloaded from:

https://github.com/LuXuefei/Nonparametric-estimation-of-probabilistic-sensitivity-measures, along with the simulated data to reproduce results.