## Appendix A Polynomial Chaos Expansion
The author leverage this fact to propose an arbitrary polynomial chaos expansion (aPCE), that can operate with probability measures that may be implicitly and incompletely defined via their statistical moments. Using aPCE, one can build the orthonormal polynomials even in the absence of the exact probability density function \(f_{\mathbf{X}}(x)\). Calculation of the coefficientsTo approximate the \(c_{\mathbf{\alpha}}\) coefficients in the expansion in Equation (10), one can use a least square minimization approach. But many predictors (\(\Psi_{\mathbf{\alpha}}(\mathbf{X})\) in (10)) make little to no contribution to fitting PCE to the computational model responses and this can result in over-fitting. To circumvent this issue, one can exploit the concept of _sparsity_. It is an attractive concept, which offers elegant complexity control, over-fitting control, feature extraction, the potential for characterization of meaningful input variables along with the practical computational speed and compactness. As stated earlier, PCEs also belong to linear regression models and employing the concept of sparsity can lead to zero values for many predictors in Equation (10). Here, we use a Bayesian sparse learning approach, in that we impose a sparsity-inducing prior on the coefficients (weights) of the predictors (\(\Psi_{x_{i}}\) in (10)), whose parameters are considered to be random variables with a hyperprior. Then, the posterior of the weights are inferred e.g., using a fast marginal likelihood maximization algorithm [61]. This learning process leads to extremely sparse inferred predictors, since they yield relatively few non-zero \(c_{i}\) coefficients. That means a significant number of the predictors give posterior distributions centered at zero. Let the target variable be \(\mathrm{Y}\), which is given by a deterministic function \(y(\mathbf{X},\mathbf{c})\) with an additive Gaussian noise which reads:  
\[\mathrm{Y}=y(\mathbf{X},\mathbf{c})+\epsilon, \tag{13}\]where \(\epsilon\) is a zero mean Gaussian random variable with precision (inverse variance) \(\beta\). Hence, the equation can be cast, in the probabilistic sense, as:  
\[p(\mathrm{Y}|\mathbf{X},\mathbf{c},\beta)=\mathcal{N}(\mathrm{Y}|y(\mathbf{X}, \mathbf{c}),\beta^{-1}). \tag{10}\]  
Additionally, let \(\mathbf{X}=\{x_{1},...,x_{N}\}\) be a data set of inputs with the corresponding model responses \(\mathrm{Y}=\{\mathcal{Y}_{1},...,\mathcal{Y}_{N}\}\). One can group the model responses into a column vector denoted by \(\mathbf{Y}\) to be distinguished from a single observation of a multivariate response, which would be denoted by \(\mathrm{Y}\). Assuming these data points are drawn independently from the distribution in Equation (10), and using Equation (11), a multivariate Gaussian likelihood function can be derived as:  
\[ \begin{split} p(\mathbf{Y}|\mathbf{X},\mathbf{c},\beta)& =\prod_{n=1}^{N}\mathcal{N}(y_{n}|\mathbf{c}^{\top}\Psi(\mathbf{x} _{n}),\beta^{-1})\\ &=\Big{(}2\pi\beta^{-1}\Big{)}^{-N/2}\exp\biggl{\{}-\frac{\beta}{2 }||y_{n}-\mathbf{c}^{\top}\Psi(\mathbf{x}_{n})||^{2}\biggr{\}},\end{split}  \tag{11}\]  
where this is a function of the parameter \(\mathbf{c}\) and \(\beta\). We introduce a Gaussian prior distribution over the parameter vector \(\mathbf{c}\) by giving each of the weight parameters \(\mathbf{c}_{i}\) a separate hyper-parameter \(\alpha_{i}\). Thus, the prior of the weights reads as:  
\[ \begin{split} p(\mathbf{c}|\mathbf{\alpha})&=\prod_{i= 1}^{P}\mathcal{N}(c_{i}|0,\alpha_{i}^{-1})\\ &=\prod_{i=1}^{P}\biggl{[}2\pi^{-1/2}\alpha_{i}^{1/2}\exp\biggl{\{} -\frac{1}{2}\alpha_{i}c_{i}^{2}\biggr{\}}\biggr{]},\end{split}  \tag{12}\]  
where \(\mathbf{\alpha}=\{\alpha_{i},...,\alpha_{P}\}^{\top}\) denotes the precision of the prior over its associated weight parameter \(\mathbf{c}\). The form of prior is ultimately responsible for the sparsity properties of the model (for more details, see Reference [62]).