Nonparametric estimation of probabilistic sensitivity measures

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 first, propose new estimators based on a one-sample design and building on the idea of placing piecewise constant Bayesian priors on the conditional distributions of the output given each input, after partitioning the input space. We then present two alternatives, based on Bayesian non-parametric density estimation, which bypass the need for predefined partitions. Quantification of uncertainty in the estimation process through is possible without requiring additional simulator evaluations via Bootstrap in the simplest proposal, or from the posterior distribution over the sensitivity measures, when the entire inferential procedure is Bayesian. 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.

Appendix

1.1 Given data estimators for the sensitivity measures in Table 1

The one-sample estimator of $\eta _i$ used here 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,\ldots ,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:

$$\begin{aligned} {\widehat{\eta }}_i^\star = \sum _{m=1}^{M}\frac{n^i_m}{n} \frac{({\bar{y}}^i_m - {\bar{y}})^2}{s^2_y}. \end{aligned}$$

(17)

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

$$\begin{aligned} {\widehat{\delta }}_i^{\star } = \sum ^{M}_{m=1} \frac{n^i_m}{n} \int _{{\mathcal {Y}}}|{\hat{f}}^{\star }_{Y}(y)-{\hat{f}}^i_m (y)|\text {d}y, \end{aligned}$$

(18)

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:

$$\begin{aligned} {{\widehat{\delta }}}_i^\diamond= & {} \sum _{m=1}^M \frac{n^i_m}{n}\,\sum _{t=1}^{T^i_m} \left( \hat{F}^i_m(\hat{b}^i_m(t)) -\hat{F}^i_m(\hat{a}^i_m(t)) \right) \nonumber \\&-\left( \hat{F}_Y(\hat{b}^i_m(t)) -\hat{F}_Y(\hat{a}^i_m(t)) \right) . \end{aligned}$$

(19)

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:

$$\begin{aligned} {\widehat{\beta }}_{i}^{\diamond }=\sum _{m=1}^{M}\frac{n^i_m}{n}\max _{j\in \{1,\ldots ,n\}} \left| {\hat{F}}_Y(y_j)-{\hat{F}}^i_m(y_j) \right| , \end{aligned}$$

(20)

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

$$\begin{aligned} \hat{F}_{Y}(y)= & {} \frac{1}{n}\sum _{j=1}^{n}{\mathbb {1}}_{(-\infty ,y_j]}(y);\nonumber \\ {\hat{F}}^i_m(y)= & {} \frac{1}{n^i_m}\sum _{y_j \in {\mathbf {y}}^i_m }\mathbb {1}_{(-\infty ,y_j]}(y), \end{aligned}$$

(21)

and $\mathbb {1}_A(y)$ denotes the indicator function, taking the value 1 if $y\in A$ and 0 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:

$$\begin{aligned} {\widehat{\eta }}_{i}^{\Diamond } = \sum _{m=1}^{M}\frac{n^i_{m}}{n}\, \frac{\left( \int _{\mathcal {Y}} {\hat{F}}^i_m(y)- \hat{F}_{Y}(y) \text {d}y \right) ^2 }{{{\widehat{\sigma }}}^2_Y}. \end{aligned}$$

(22)

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.

1.2 Numerical experiments for the partition selection problem

Fig. 5

RMSE of sensitivity measures estimates for $X^1$ of the 2-input simulator in Eq. (16). Magenta lines correspond to $M=2.5\root 3 \of {n}$; $n \in [300,900]$, $M\in [5,34]$

Fig. 6

RMSE of sensitivity measures estimates for $X^3$ of the 21-input simulator in Eq. (7). Magenta lines correspond to $M=2.5\root 3 \of {n}$; $n \in [300,900]$, $M\in [5,34]$

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. (17), (18) and (20)]. 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 \root 3 \of {n}$.

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

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

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 5 and 6 present the heatplot of RMSEs in percentage ($\hbox {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. 5a, dark (blue) refers to low RMSE (less than 10%), and light (red) to relative high RMSE (higher than 14%). The magenta line maps n into M using the previously mentioned heuristic function. Figure 5 shows that the proposed heuristic works well on the 2-input simulator [Eq. (16)], with the magenta line falling mainly into dark coloured regions. However, for the 21-input simulator [Eq. (7)] we would incur in high errors at small sample sizes. For instance consider graph a) in Fig. 6. 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. (16)], the ideal partition size for the variance-based estimator is between 10 to 15 (Fig. 5b), while for mutual information, if falls between 20 to 25 (Fig. 5c).

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

1.3 Details on the non-parametric estimators

We present further details regarding the implementation of the Bayesian non-parametric estimation methods in Sects. 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}}$.

Partition-dependent Bu and Pu estimation. Note that model 3 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,\ldots ,k$ and $m=1,2,\ldots ,M$. In fact,

$$\begin{aligned} {\mathbb {P}}_{Y}(\cdot |{\mathbb {P}}_{1:M}^{i})=\sum \limits _{m=1}^{M}{\mathbb {P}} _{m}^{i}(\cdot ){\mathbb {P}}(X^i\in {\mathcal {X}}^i_m)=\frac{1}{M} \sum \limits _{m=1}^{M}{\mathbb {P}}_{m}^{i}(\cdot )\text {.} \end{aligned}$$

Then, by marginalizing, we obtain

$$\begin{aligned} {\mathbb {P}}_{Y}(\cdot |\alpha G)= & {} \frac{1}{M}\sum \limits _{m=1}^{M}\int {\mathbb {P}}_{m}^{i}(\cdot )\text {d}{{\mathcal {D}}}{{\mathcal {P}}}({\mathbb {P}}_{m}^{i}|\alpha G)\\= & {} \int {\mathbb {P}}(\cdot )\text {d}{{\mathcal {D}}}{{\mathcal {P}}}({\mathbb {P}}|\alpha G), \end{aligned}$$

because $\int {\mathbb {P}}_{m}^{i}(\cdot )$d${{\mathcal {D}}}{{\mathcal {P}}}({\mathbb {P}}_{m}^{i}|\alpha G)$ does not depend on i or m. In other words, a priori, ${\mathbb {P}}_{Y} \overset{}{\sim }{{\mathcal {D}}}{{\mathcal {P}}}(\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 of the marginal for Y can be obtained as:

$$\begin{aligned} {\mathbb {P}}_{Y}(\cdot |\alpha G,D^i)=\frac{1}{M}\sum \limits _{m=1}^{M}\int {\mathbb {P}}_{m}^{i}(\cdot )\text {d}{{\mathcal {D}}}{{\mathcal {P}}}({\mathbb {P}}_{m}^{i}|(\alpha +n_{m}^{i}) {\widetilde{G}}_{m}^{i}),\nonumber \\ \end{aligned}$$

(23)

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 _m (1/n_{m}^{i})\delta ^{Dirac}_y\simeq \sum _m (M/n)\delta ^{Dirac}_y$. Thus, asymptotically, ${\mathbb {P}}_{Y}(\cdot |\alpha G,D^i)$ does not depend on m or i and ${\mathbb {P}}_{Y}(\cdot |\alpha G,D^i)\overset{}{\sim } {{\mathcal {D}}}{{\mathcal {P}}}((\alpha +n){\widetilde{G}})$, where

$$\begin{aligned} {\widetilde{G}}=\dfrac{\alpha }{\alpha +n}G+\dfrac{n}{\alpha +n} {\widehat{\mathbb {P}}}_n, \end{aligned}$$

(24)

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

For the numerical experiments in Sects. 3 and 4.3, the input data, ${\mathbf {x}}$, are obtained by transforming a 21-dimensional standard Gaussian sample generated via a Halton Quasi random sequence via the Cholesky decomposition of the covariance matrix. Specifically, MATLAB^© functions haltonset and chol are used. The mass parameter, $\alpha$, for the DP prior is set equal to $0.1\,n/M$ throughout. The base measure, G, is a Normal distribution with hyper-parameters fixed via an empirical approach, based on the available sample ${\mathbf {y}}$. ${\mathbf {y}}_m^i$. Overall, these choices centre 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. Algorithm 1 details the inferential procedure. Note that the calculations of $\eta ^{C,s}$, $\delta ^{C,s}$ and $\beta ^{C,s}$ are equivalent to the pdf-based estimators in Eqs. (17)–(20) but with the enriched samples. Alternatively, the cdf-based estimators in Eqs. (22) and (19) could be used for $\eta ^{C,s}$ and $\delta ^{C,s}$.

Partition-free BNJ estimation. Following the proposal in Jara et al. (2011), we choose G to be a Normal-Inverse Wishart distribution

$$\begin{aligned}&(\mu _\ell , \varSigma _\ell )| (m_1, \gamma , \psi _1){\mathop {\sim }\limits ^{iid}}{\mathcal {N}}\left( \mu _{\ell } \vert m_1, \frac{1}{\gamma } \varSigma \right) IW(\varSigma _{\ell } \vert 4, \psi _1),\\&\ell =1,2,\ldots \end{aligned}$$

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:

$$\begin{aligned}&\gamma \sim \text {Gamma} \left( \cdot |0.5, 0.5\right) ,\; m_1 \vert (m_2, s_2) \sim {\mathcal {N}}(\cdot |m_2, s_2),\\&\psi _1 \vert (s_2) \sim IW (\cdot |4, s_2^{-1}), \end{aligned}$$

where $\text {Gamma} (\cdot |a_1,a_2)$ denotes the Gamma distribution with mean $a_1/a_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{\varSigma }}^s)$, $s=1,\ldots ,S$. In practice, the number $J_s$ of components with non-zero weights is finite, thus we have

$$\begin{aligned} {\underline{w}}^s= & {} (w_1^s, \dots , w_{J_s}^s), \quad \nonumber \\ {\underline{\mu }}^s= & {} (\mu _1^s, \dots , \mu _{J_s}^s), \quad {\underline{\varSigma }}^s = (\varSigma _1^s, \dots , \varSigma _{J_s}^s), \nonumber \\ \text {with} \qquad \mu _{\ell }^s= & {} \left[ {\begin{array}{c} \mu _{1,\ell }^s \\ \mu _{2,\ell }^s \\ \end{array} } \right] , \quad \varSigma _{\ell }^s = \left[ {\begin{array}{cc} \sigma _{1,\ell }^s &{} \sigma _{3,\ell }^s \\ \sigma _{3,\ell }^s &{} \sigma _{2,\ell }^s\\ \end{array} } \right] . \end{aligned}$$

(25)

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

$$\begin{aligned} f^{BNJ,s}_{X,Y} (x,y|{\underline{\theta }}^s) = \sum _{\ell =1}^{J_s} w_\ell ^s \cdot {\mathcal {N}} (x,y|\mu _\ell ^s,\varSigma _\ell ^s). \end{aligned}$$

(26)

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

$$\begin{aligned} f^{BNJ,s}_{Y} (y \vert {\underline{\theta }}^s)= & {} \sum _{\ell =1}^{J_s} w_\ell ^s \cdot {\mathcal {N}} (y|\mu _{2,\ell }^s,\sigma _{2,\ell }^s), \nonumber \\ f^{BNJ,s}_{Y\vert x} (y \vert x,{\underline{\theta }}^s )= & {} \sum _{\ell =1}^{J_s} w_\ell ^s \cdot {\mathcal {N}} \left( \cdot |\nu _{2,\ell }^s,\tau _{2,\ell }^s \right) \end{aligned}$$

(27)

where $\nu _{2,\ell }^s = \mu _{2,\ell }^s + \sigma _{3,\ell }^s (x - \mu _{1,\ell }^s )/\sigma _{1,\ell }^s$ and $\tau _{2,\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 \vert X}$, as well as the marginal mean and variance can be calculated trivially. In particular,

$$\begin{aligned} \mu _{Y}^s :={\mathbb {E}}[Y \vert {\underline{\theta }}^s]= & {} \sum _{\ell =1}^{J_s} w_{\ell }^s \mu _{2,\ell }^s, \nonumber \\ V_Y^s :={\mathbb {V}}[Y \vert {\underline{\theta }}^s]= & {} \sum _{\ell =1}^{J_s} w_{\ell }^s \left( \sigma _{2,\ell }^s + \left( \mu _{Y}^s - \mu _{2,\ell }^s \right) ^2 \right) . \end{aligned}$$

(28)

For the numerical experiments in Sects. 3 and 4.3, the mass parameter is fixed to $\alpha = 1$, and the remaining hyperparameters are set to $m_2 = (\mu _{X}, {{\bar{y}}})$ and $s_2 = \text {diag}(\sigma _{X}^2, s_y^2)$. The Pseudo-code of BNJ estimation for $\eta , \delta$ and $\beta$ is illustrated in Algorithm 2.

Partition-free BNC 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:

$$\begin{aligned}&({\mathbf {b}}_\ell ,\sigma _{\ell },\mu _\ell ){\mathop {\sim }\limits ^{iid}}\\&{\mathcal {N}}\left( {\mathbf {b}}_\ell \mid {\mathbf {b}}_0, \sigma _{\ell } C^{-1} \right) \text {Gamma} (\sigma _{\ell }^{-1}\mid 1,1) {\mathcal {N}}\\&\quad \left( \mu _\ell \mid \mu _0, \left( \tau /10 \right) ^{-1} \right) , \end{aligned}$$

where ${\mathbf {b}}_\ell = (a_\ell ,b_\ell )$ and $\tau \sim \text {Gamma} (\cdot \mid 1,1)$. The hyperparameters are chosen empirically. As an illustration, consider the 21-input simulator. Figure 7 shows the scatter-plot of $({\mathbf {x}}^3,{\mathbf {y}})$ and the corresponding 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

Fig. 7

Scatter plot of Y and $X^3$ for the 21-input simulator. Red lines constitute the convex hull of $\{(x^3_1, y_1),\ldots ,(x^3_n, y_n)\}$. The bold blue and red lines are used for prior specification

$$\begin{aligned} {\underline{a}}^s= & {} (a_1^s, \dots , a_{J_s}^s), \nonumber \\ {\underline{b}}^s= & {} (b_1^s, \dots , b_{J_s}^s), \quad {\underline{\sigma }}^s = (\sigma _1^s, \dots , \sigma _{J_s}^s),\nonumber \\ {\underline{\omega }}^s= & {} (\omega _1^s, \dots , \omega _{J_s}^s), \quad {\underline{\mu }}^s = (\mu _1^s, \dots , \mu _{J_s}^s). \end{aligned}$$

(29)

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

$$\begin{aligned} f^{BNC,s}_{Y|X} (y|x,{\underline{\theta }}^s,{\underline{\psi }}^s) = \sum _{\ell =1}^{J_s} w_\ell ^s(x) {\mathcal {N}}\left( y|a_\ell ^s+b_\ell ^s x, \sigma _\ell ^{s} \right) . \end{aligned}$$

(30)

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

$$\begin{aligned} f^{BNC,s}_{Y}(y \vert {\underline{\theta }}^s) \approx \int _{{\mathcal {X}}} f^{BNC,s}_{Y|X} f_{X} \text {d}x. \end{aligned}$$

(31)

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

$$\begin{aligned} \mu _Y^s:= & {} {\mathbb {E}}[Y\vert {\underline{\theta }}^s,{\underline{\psi }}^s] \approx \int _{{\mathcal {Y}}} y f^{BNC,s}_{Y} \text {d}y, \\ V_{Y}^s:= & {} {\mathbb {V}}[Y \vert {\underline{\theta }}^s,{\underline{\psi }}^s] \approx \int _{{\mathcal {Y}}} \left( y - \mu _Y^s \right) ^2 f^{BNC,s}_{Y} \text {d}y \end{aligned}$$

The Pseudo-code of BNC estimation for $\eta , \delta$ and $\beta$ is illustrated in Algorithm 3.