{'fraction_non_alphanumeric': 0.03674233559560028, 'fraction_numerical': 0.010680169350892495, 'mean_word_length': 4.475131042516017, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic polynomial chaos expansions to emulate stochastic simulators

In the context of uncertainty quantification, computational models are required to be repeatedly evaluated. This task is intractable for costly numerical models. Such a problem turns out to be even more severe for stochastic simulators, the output of which is a random variable for a given set of input parameters. To alleviate the computational burden, surrogate models are usually constructed and evaluated instead. However, due to the random nature of the model response, classical surrogate models cannot be applied directly to the emulation of stochastic simulators. To efficiently represent the probability distribution of the model output for any given input values, we develop a new stochastic surrogate model called stochastic polynomial chaos expansions. To this aim, we introduce a latent variable and an additional noise variable, on top of the well-defined input variables, to reproduce the stochasticity. As a result, for a given set of input parameters, the model output is given by a function of the latent variable with an additive noise, thus a random variable. In this paper, we propose an adaptive algorithm which does not require repeated runs of the simulator for the same input parameters. The performance of the proposed method is compared with the generalized lambda model and a state-of-the-art kernel estimator on two case studies in mathematical finance and epidemiology and on an analytical example whose response distribution is bimodal. The results show that the proposed method is able to accurately represent general response distributions, i.e., not only normal or unimodal ones. In terms of accuracy, it generally outperforms both the generalized lambda model and the kernel density estimator.


Introduction
In modern engineering, computational models, a.k.a. simulators, are commonly used to simulate different operational scenarios of complex systems in silico. These models help engineers assess the reliability, control the risk, and optimize the system components in the design phase. Conventional simulators are usually deterministic: a given set of input parameters has a unique corresponding model response. In other words, repeated model evaluations with the same input values will always give identical results. In contrast, stochastic simulators return different outcomes of the model response when run twice with the same input parameters.
Stochastic simulators are widely used in engineering and applied science. The intrinsic stochasticity typically represents some uncontrollable effect in the system [1,2]. For example, in mathematical finance, Brownian motions are commonly introduced to represent stochastic effects and volatility of the stock market [1]. In epidemic simulations, additional random variables on top of the well-defined characteristic values of the population are used to simulate the stochastic spread of a disease [2].
Mathematically, a stochastic simulator can be viewed as a function where D X is the domain of the input parameters, and Ω denotes the probability space that represents the internal stochasticity. The latter is due to some latent random variables Ξ(ω) which are not explicitly considered as a part of the input variables. The stochastic simulator can then be considered as a deterministic function of the input vector x and the latent variables Ξ.
However, it is assumed that one can only control x but not Ξ when evaluating the model. Hence, when the value of x is fixed but Ξ is generated randomly following the underlying probability distribution, the output remains random.
In practice, each model evaluation for a fixed vector of input parameters x 0 uses a particular realization of the latent variables, i.e., a particular ω 0 ∈ Ω that is usually controlled by the random seed. Thus, it provides only one realization of the output random variable. In order to fully characterize the associated distribution of M s (x 0 , ·), it is necessary to repeatedly run the stochastic simulator with the same input parameters x 0 . The various output values obtained by this procedure are called replications in the sequel.
In the context of uncertainty quantification or optimization, various input values should be investigated. To this aim, multiple runs of the simulator are needed for many different inputs and for many replications. This becomes impracticable for high-fidelity costly numerical models.
In this context, surrogate models have received tremendous attention in the past two decades. A surrogate model is a proxy of the original model constructed from a limited number of model runs. However, standard surrogate models such as polynomial chaos expansions [3] and Gaussian processes [4] that have been successfully developed for deterministic simulators are not directly applicable to emulating stochastic simulators due to the random nature of the latter.
In the past decade, large efforts have been dedicated to estimating some summary quantities of the response distribution which are deterministic functions of the input.
For the mean and variance of the response distribution, Ankenman et al. [5] proposed using replications to estimate the mean and variance for various input values. The mean function is represented by a Gaussian process, for which the variance estimated from the replications is cast as a heteroskedastic effect. Marrel et al. [6] modeled both the mean and variance by Gaussian processes. The estimation procedure is similar to the feasible generalized least-squares [7] that consists in alternating between fitting the mean from the data and the variance from the residuals. This approach does not require replications. Binois et al. [8] proposed jointly optimizing the likelihood to represent the mean and variance by Gaussian processes, which is mainly designed for data with replications.
To estimate the quantiles of the response distribution, Koenker and Bassett [9] proposed optimizing the check function, which established the quantile regression method. Plumlee and Tuo [10] suggested estimating the quantiles by performing replications and building a Gaussian process from the estimated quantiles. The reader is referred to Torossian et al. [11] for a detailed review.
The methods listed above produce only targeted summary quantities. However, far less literature has been devoted to the emulation of the entire probability distribution function of the response random variable for a given input. Three types of methods can be found in the literature.
Moutoussamy et al. [12] proposed using replications to characterize the response distribution for different input values. Then, the fitted distributions (based on replications) for the discrete input values can be extended to the entire input space by parametric or nonparametric techniques.
Since this approach capitalizes on replications for local inference, it is necessary to generate many replications to obtain an accurate surrogate [13], i.e., in the order of 10 3 − 10 4 [14].
In the second approach, a stochastic simulator is considered as a random field indexed by the input variables [15,16]. When fixing the internal stochasticity ω in Eq. (1), the stochastic simulator is a mere deterministic function of x, called a trajectory. This function can be emulated by standard surrogate methods. Collecting different trajectories, one can approximate the underlying random field using Karhunen-Loève expansions. Therefore, it is necessary to fix the internal randomness to apply this approach, which is practically achieved by controlling the random seed.
The third type of methods is referred to as the statistical approach and does not require replications or manipulating the random seed. If the response distribution belongs to the exponential family, generalized linear models [17] and generalized additive models [18] can be efficiently applied. For arbitrary types of response distributions, nonparametric estimators developed in statistics can be applied, namely kernel density estimators [19,20] and projection estimators [21]. However, nonparametric estimators are known to suffer from the curse of dimensionality, which indicates that the necessary amount of data increases drastically with increasing input dimensionality. To balance between very restrictive parametric assumptions and nonparametric approaches, Zhu and Sudret [22,23] proposed using generalized lambda distributions to approximate the response distributions. The four distribution parameters are seen as functions of the input and further represented by polynomial chaos expansions. The main limitation of this approach is that it cannot produce multimodal distributions, however.
In this paper, we develop an original approach that directly emulates the functional representation in Eq. (1). More precisely, we extend the classical polynomial chaos expansions to emulating stochastic simulators. We introduce a latent variable and a noise variable to reproduce the random behavior of the model output. We develop an adaptive method to construct such a surrogate model. This novel stochastic surrogate is parametric and shown to be not limited to unimodal distributions.
The remainder of the paper is organized as follows. In Section 2, we first review the standard polynomial chaos representations. In Section 3, we present a novel formulation named stochastic polynomial chaos expansions which is meant for stochastic simulators. In Section 4, we present the algorithms to adaptively build such a surrogate from data without the need for replications.
We illustrate the performance of the proposed method on a complex analytical example and on case studies from mathematical finance and epidemiology in Section 5. Finally, we conclude the main findings of the paper and provide outlooks for future research in Section 6.

Reminder on polynomial chaos expansions
Polynomial chaos expansions (PCEs) have been widely used in the last two decades to emulate the response of deterministic simulators in many fields of applied science and engineering. Consider a deterministic model M d which is a function that maps the input parameters x = (x 1 , x 2 , . . . , x M ) T ∈ D X ⊂ R M to the scalar output y = M d (x) ∈ R. In the context of uncertainty quantification, the input vector x is affected by uncertainties and thus modeled by a random vector X with prescribed joint probability density function (PDF) denoted by f X . In the sequel, we focus on the case where the input parameters are independent for simplicity. Therefore, the joint PDF is expressed by where f X j is the marginal PDF of the input random variable X j . Note that in the case where the input vector X has dependent components, it is always possible to transform them into independent ones using the Nataf or Rosenblatt transform [24][25][26].
Because of the randomness in the input, the model response Under certain conditions on the joint PDF f X [27], the Hilbert space H possesses a polynomial basis. As a result, M d can be represented by an orthogonal series expansion where c α is the coefficient associated with the basis function ψ α that is defined by the multi-index α. More precisely, the multivariate basis function ψ α is given by a tensor product of univariante where α j indicates the degree of ψ α (x) in its j-th component x j , and φ (j) k : k ∈ N is the orthogonal polynomial basis with respect to the marginal distribution f X j of X j , which satisfies In the equation above, the Kronecker symbol δ kl is such that δ kl = 1 if k = l and δ kl = 0 otherwise.
Following Eq. (5), the multivariate polynomial basis is defined from univariate orthogonal polynomials that depend on the corresponding marginal distribution. For uniform, normal, gamma and beta distributions, the associated orthogonal polynomial families are known analytically [28].
For arbitrary marginal distributions, such a basis can be iteratively computed by the Stieltjes procedure [29].
The spectral representation in Eq. (4) involves an infinite sum of terms. In practice, the series needs to be truncated to a finite sum. The standard truncation scheme is defined by selecting all the polynomials whose total degree is small than a given value p, i.e., A p,M = α ∈ N M , M j=1 α j ≤ p . However, this will provide a large number of terms for big values of p and M . A more flexible scheme is the hyperbolic (q-norm) truncation scheme [30]: where p is the maximum polynomial degree, and q ∈ (0, 1] defines the quasi-norm α q = M j=1 |α j | q 1/q . This truncation scheme allows excluding high-order interactions among the input variables but keeps univariate effects up to degree p. Note that with q = 1, we recover the full basis of total degree less than p. To estimate the coefficients in Eq. (4), one popular approach relies on minimizing the mean-squared error between the model response and the surrogate model. The basic method applies ordinary least-squares (OLS) with a given set of basis (e.g., defined by a truncation scheme) [31].
In this approach, the model is evaluated on a number of points called the experimental design (1) , . . . , x (N ) . The associated model responses are gathered into y = y (1) , . . . , y (N ) with y (i) = M x (i) . The basis functions (and thus the coefficients) can be arranged by ordering the multi-indices {α j } P j=1 . The regression matrix Ψ is defined by Ψ ij = ψ α j x (i) . By minimizing the mean-squared error between the original model and the surrogate on the experimental design, the OLS estimator is given bŷ With increasing polynomial degree or input dimension, the number of coefficients increases drastically. As a consequence, a large number of models runs are necessary to guarantee a good accuracy, which becomes intractable for costly simulators. To solve this problem, Blatmann and Sudret [30], Doostan and Owhadi [32], Babacan et al. [33] developed methods to build sparse PCEs by only selecting the most influential polynomials. The reader is referred to the review papers by Lüthen et al. [34,35] for more details.

Introduction
Let us now come back to stochastic simulators. It would be desirable to have a spectral expansion such as Eq. (4) for stochastic simulators. Indeed, the standard PCE has numerous features such as close-to-zero-cost model evaluations, and clear interpretation of the coefficients in terms of sensitivity analysis [36]. However, because the spectral expansion in Eq. (4) is a deterministic function of the input parameters, it cannot be directly used to emulate stochastic simulators.
Considering the randomness in the input variables, the output of a stochastic simulator is a random variable. The randomness of the latter comes from both the intrinsic stochasticity and the uncertain inputs. When fixing the input parameters, the model response remains random.
For the purpose of clarity, we denote by Y x the random model response for the input parameters x and by Y the model output containing all the uncertainties: following Eq. (1), we have From a probabilistic perspective, Y x is equivalent to the conditional random variable Y | X = x.
Let F Y |X (y | x) denote the associated cumulative distribution function (CDF). By using the probability integral transform, we can transform any continuous random variable Z to the desired distribution, that is where F Z is the CDF of Z. The equality in Eq. (10) is to be understood in distribution, meaning that two random variables on the left-and right-hand side follow the same distribution. In Eq. (10), the right-hand side is a deterministic function of both x and z. As a result, assuming that Y has a finite variance, we can represent this function using a PCE in the (X, Z) space, that is, For a given vector of input parameters x, the expansion is a function of the artificial latent variable Z, thus a random variable Then, we apply a truncation scheme A (e.g., Eq. (7)) to reduce Eq. (12) to a finite sum Even though Eq. (13) is derived from Eq. (11), it is more general. Equation (10) offers one way to represent the response distribution by a transform of a latent variable. But many other transforms can achieve the same goal. For example, using Z ∼ N (0, 1), both µ(x) + σ(x)Z and ). Because we are interested in the response distribution, Eq. (13) only requires that the polynomial transform of the latent variable produces a distribution that is close to the response distribution, but the transform does not need to follow Eq. (11) exactly. Note that the latent variable Z is only introduced to reproduce the stochasticity, but it does not allow us to represent the detailed data generating process of the simulator though. In other words, the PCE in Eq. (13) cannot emulate the response for a particular replication, yet it provides a representation of the distribution of Y x .

Potential issues with the formulation in Eq. (13)
Building a PCE by least-squares as presented in Section 2 requires evaluating the deterministic function to surrogate, which, in the case of stochastic simulators, is the left-hand side of Eq. (11).
However, it is practically impossible to evaluate such a function, as the response distribution F −1 Y |X is unknown. One common way to fit the latent variable model defined in Eq. (13) is maximum likelihood estimation [37,38]. In this section, we show some potential problems associated with a standard use of this method for building Eq. (13), which calls for a novel fitting algorithm.
According to the definition in Eq. (13),Ỹ x is a function of Z. Denote f Z (z) the PDF of Z and D Z the support of Z. Based on a change of variable [39], we can obtain the PDF ofỸ x , which is denoted by fỸ x (y; x, c). As a result, the (conditional) likelihood function of the coefficients c for a data point (x, y) is given by Now, let us consider an experimental design X = x (1) , . . . , x (N ) . The stochastic simulator is assumed to be evaluated once for each point Note that here we do not control the random seed, so the model outcomes for different values of x are independent. Thus, the likelihood function can be computed by the product of l c; x (i) , y (i) over the N data points. As a result, the maximum likelihood estimator is given byĉ Equation (15) commonly serves as a basic approach for fitting parametric statistical models (including stochastic surrogates) [17,22,40]. However, the likelihood function of the latent PCE defined in Eq. (13) is unbounded and can reach +∞, making the maximization problem Eq. (15) ill-posed.
To illustrate the issue, let us consider a simple stochastic simulator without input variables, which gives a realization of Y upon each model evaluation. Hence, the surrogate in Eq. (13) contains only the latent variable Z, that is, be a second-degree polynomial expressed by monomials g(z) = a 1 z 2 + a 2 z + a 3 . Note that there is a one-to-one mapping between monomials and full polynomial chaos basis, so one can map a = (a 1 , a 2 , a 3 ) T to c through a change of basis. Using a change of variable [39], the PDF ofỸ is where 1 is the indicator function, and g denotes the derivative of g. For a given y 0 , certain choices of a can make any given z 0 with f Z (z 0 ) = 0 satisfy g(z 0 ) = y 0 and g (z 0 ) = 0: The system of equations in Eq. (17) is underdetermined for a. Therefore, there are infinite combinations of the coefficients a, and therefore of c, such that the denominator of Eq. (16) is zero and the numerator is non-zero, which gives fỸ (y 0 ) = +∞. Consequently, the maximum likelihood estimation will always produce a certain vector c that makes the likelihood reach +∞.
As a conclusion, the surrogate ansatz of Eq. (13) can produce non-smooth conditional PDFs with singularity points where fỸ x tends to infinity. Consequently, the standard maximum likelihood estimation would fail.

Formulation of stochastic polynomial chaos expansions
In the previous section, we discussed some potential problems of the model defined in Eq. (13).
To regularize the optimization problem in Eq. (15) and smooth out the produced PDFs, we introduce an additive noise variable , and define the stochastic surrogate as follows: where is a centered Gaussian random variable with standard deviation σ, i.e., ∼ N (0, σ 2 ).
With this new formulation, the response PDF of the stochastic surrogate is a convolution of that of the PCE and the Gaussian PDF of .
Using Hölder's inequality, the above integral is bounded from above by meaning that the PDF ofỸ x and the associated likelihood function are bounded.
To illustrate the role of the additive noise variable in Eq. (18), let us consider a random variable Y with bimodal distribution to be represented by where the latent variable Z follows a standard normal distribution and ∼ N (0, σ). In the case σ = 0 (the noise term vanishes), we build the model by applying a standard algorithm such as least-angle regression (LAR) [30] to the probability integral transform F −1 Y (F Z (Z)). When the regularization term is added, maximum likelihood estimation can be used (see Section 4.1 for details) to construct the surrogate. Figure 1 shows the original (reference) PDF, and the ones obtained by LAR (σ = 0) and by the stochastic PCE for two different values of σ. It is observed that the PDF obtained by LAR has singularity points, which confirms the analysis in Section 3.2, whereas the proposed noise term regularizes the PDFs. Moreover, LAR is applied directly to the probability integral transform which in practice is unknown. In contrast, the maximum likelihood estimation does not require knowing the values of Z (in this example, only the realizations of Y are used). Finally, the value of σ affects the accuracy of the model. Hence, σ is an additional parameter of the model that must also be fitted to the data to get the optimal approximation. The fitting procedure is detailed in the next section.

Fitting the stochastic polynomial chaos expansion
To construct a stochastic PCE defined in Eq. (18), one needs to estimate both the coefficients c and the standard deviation σ of the noise variable. In this section, we present a method to calibrate these parameters from data without replications. Moreover, we propose an algorithm that adaptively selects an appropriate distribution for the latent variable Z and truncation scheme A.

Maximum likelihood estimation
Let us assume for a moment that the standard deviation σ of the noise variable is given (the estimation of σ will be investigated separately in Section 4.4). From Eq. (18), we see that our surrogate responseỸ x is the sum of a polynomial function of (x, z) and the noise variable .
Therefore, its PDF can be computed by sinceỸ x | Z = z is a Gaussian random variable with mean value α∈A c α ψ α (x, z) and variance σ 2 according to Eq. (18). In this equation, ϕ stands for the standard normal PDF. Therefore, for a given data point (x, y), the likelihood of the parameters c conditioned on σ reads In practice, we can use numerical integration schemes, namely Gaussian quadrature [41], to efficiently evaluate this one-dimensional integral, that is where N Q is the number of integration points, z j is the j-th integration point, and w j is the corresponding weight, both associated to the weight function f Z . Based on Eq. (24) and the available data (X , y), the PCE coefficients c can be fitted using the maximum likelihood The gradient of Eq. (24), and therefore of Eq. (25), can be derived analytically. Hence, we opt for the derivative-based BFGS quasi-Newton method [42] to solve this optimization problem.

Starting point for the optimization
The objective function to optimize in Eq. (25) is highly nonlinear. As a result, a good starting point is necessary to ensure convergence. According to the properties of the polynomial chaos basis functions, the mean function of a stochastic PCE can be expressed as where α z is the degree of the univariate polynomial in Z. Equation (26) contains all the terms without Z, as indicated by α z = 0. We define this set of multi-indices as Another surrogatem(x) of the mean function can be obtained by using standard (or sparse) regression to directly fit the following expansion: The obtained coefficients c m are used as initial values for the coefficients {c α : α ∈ A m } of the stochastic surrogate in the optimization procedure, i.e., c α for α ∈ A m .
For the other coefficients {c α : α ∈ A \ A m }, we randomly initialize their value.

Warm-start strategy
Because of the form of the likelihood Eq. (23), the gradient at the starting point can take extremely large values when σ is small. In this case, the optimization algorithm may become unstable and converge to an undesired local optimum. To guide the optimization, we propose a warm-start strategy summarized in Algorithm 1. We generate a decreasing sequence σ = {σ 1 , . . . , σ Ns } with σ Ns = σ (the target value). In this paper, we choose the maximum value σ 1 of the sequence as the square root of the leave-one-out error ε LOO in the mean fitting procedure (see Appendix A.1 for the explanation of this choice). Then, σ is generated equally-spaced in the log-space between √ ε LOO and σ. Starting with σ 1 which is the largest element of σ, we build a stochastic PCE

Cross-validation
As explained in Section 3.2, the hyperparameter σ cannot be jointly estimated together with the PCE coefficients c because the likelihood function can reach +∞ for certain choices of c and σ = 0. As a result, σ should be tuned separately from the estimation of c.
In this paper, we propose applying cross-validation (CV) [40] to selecting the optimal value of σ. More precisely, the data (X , y) are randomly partitioned into N cv equal-sized groups Note that the coefficients depend on the value of σ, and thus we explicitly write them as functions of σ. The validation set V k is then used to evaluate the out-of-sample performance: We repeat this procedure for each group of the partition {V k : k = 1, . . . , N cv } and sum up the respective score to estimate the generalized performance, referred to as CV score in the sequel.
Then, the optimal value of σ is selected as the one that maximizes this CV score: Because of the nested optimization in Eq. (29), the gradient of Eq. (31) is difficult to derive.
In this paper, we apply the derivative-free Bayesian optimizer [43] to solving Eq. According to our investigations, the optimal valueσ is always within the proposed interval.
After solving Eq. (31), the selectedσ is used in Eq. (25) with all the available data to build the final surrogate.
Large value of N cv can lead to high computational cost, especially when N is big. In this paper, we choose N cv = 10 for N < 200 (small data set), N cv = 5 for 200 ≤ N < 1,000 (moderate data set) and N cv = 3 for N ≥ 1,000 (big data set).

Adaptivity
The method developed in Sections 4.1 and 4.4 allows us to build a stochastic PCE for a given distribution of the latent variable Z and truncated set A of polynomial chaos basis. In principle, one can choose any continuous probability distribution for the latent variable and a large truncated set. However, in practice, certain types of latent variables may require a lot of basis functions to approximate well the shape of the response distribution. This leads to many model parameters to estimate, which would cause overfitting when only a few data are available. In this section, we propose a procedure to iteratively find a suitable distribution for the latent variable Z and truncation scheme A.
We consider N z candidate distributions D = {D 1 , . . . , D Nz } for the latent variable, N p degrees p = p 1 , . . . , p Np and N q q-norms q = p 1 , . . . , p Np that are used to define the hyperbolic truncation scheme in Eq. (7). Both p and q are sorted in increasing order.
The adaptive procedure is shown in Algorithm 2 and described here. For each type of latent variable and truncation set A = A p,q,M , we first apply the hybrid LAR algorithm developed by Blatman and Sudret [30] to fitting the mean functionm(x) as shown in Eq. (28). This algorithm only selects the most important basis among the candidate set A m defined in Eq. (27).
To reduce the total number of unknowns in the optimization Eq. (25), we exclude from A the basis functions in A m that are not selected by hybrid LAR. In other words, we only estimate the coefficients associated with the basis functions that either have α z = 0 or are selected by the hybrid LAR when fitting the mean function m(x). Then, we use the methods presented in Sections 4.1 and 4.4 to build a stochastic PCE for A and record the CV score. The latter is used for model comparisons, and the one with the best CV score is selected as the final surrogate.

Algorithm 2 Adaptive algorithm for building a stochastic PCE
Input: (X , y)

8:
A n ← Hybrid-LAR (X , y, A m ) % Selection of the basis form(x)

9:
A ← A n ∪ A c 10: Apply the algorithm presented in Sections 4.1 and 4.4 to build a stochastic PCE with A, which gives c, σ, and the CV score l ip,iq associated with σ. 11: end for 12: end for 13: end for 14: Return the model with the maximum CV score In order to avoid going through all the possible combinations, we propose a heuristic early stopping criterion for both degree and q-norm adaptivity. If two consecutive increases of q-norm cannot improve the CV score, the inner loop for q-norm adaptivity stops. Besides, if the best model (among all the q-norms) of a larger degree decreases the CV score, the algorithm stops exploring higher degrees. Note that the early stopping is only applied to p-and q-adaptivity, but all the candidate distributions are investigated.
In summary, we sketch the overall procedure (presented in Sections 4.1 to 4.5) to adaptively build a stochastic PCE from data in Figure 2.

Post-processing of stochastic polynomial chaos expansions
In this section, we show how to post-process a stochastic PCE for various analyses. The very feature of this surrogate is that it provides a functional mapping between the input parameters X, the latent variable Z, and the noise term : To generate realizations ofỸ , we simply sample X, Z and following their distributions and then evaluate Eq. (32). To obtain samples ofỸ x for a fixed x (e.g., to plot the conditional distribution), we follow the same procedure with fixed X = x. Moreover, Eq. (32) can be easily vectorized for efficient sampling.
By generating a large number of samples, one can display the distribution ofỸ andỸ x using histograms or kernel density estimation. We can also use the quadrature version in Eq. (24) to get an explicit form of the conditional response distribution ofỸ x .
In addition, because the proposed surrogate model is derived based on PCE, it inherits all the good properties of PCE. In particular, some important quantities can be directly computed by post-processing the PCE coefficients c and the parameter σ, without the need for sampling.
Indeed, the mean and variance ofỸ are given by where c 0 is the coefficient of the constant function.
As already shown in Eq. (26), for a given value of x, the mean of the model responseỸ x can be computed as Similarly, we can compute the variance as follows:

Global sensitivity analysis
In the context of global sensitivity analysis of stochastic simulators [23], various types of Sobol' indices can also be computed analytically for the proposed surrogate model. The classical Sobol' indices are defined from the Sobol'-Hoeffding decomposition of the deterministic model given by the stochastic simulator with both the well-defined input variables X and its intrinsic stochasticity as explicit inputs ω, see Eq. (1). Since the surrogate model in Eq. (32) is also a deterministic function of X and the additional variables Z and , the Sobol' indices can be efficiently computed from the PCE coefficients, similarly to the classical PCE-based Sobol' indices [36]. For example, the first-order classical Sobol' index of the i-th input X i is given by Similarly, one can also calculate higher-order and total Sobol' indices of the model Eq. (32). Let us split the input vector into two subsets The higher-order and total Sobol' indices, denoted by S u and S T i , respectively, are given by where A u def = {α ∈ A : α i = 0, α j = 0 , α z = 0 , ∀i ∈ u, ∀j ∈ ∼u}. However, as mentioned in Section 3, the surrogate model aims only at emulating the response distribution of the simulator instead of representing the detailed data generation process. Therefore, the indices involving the artificial variables introduced in the surrogate (i.e., Z and ), e.g., the first-order Sobol' index for Z and the total Sobol' index for each component of X, do not reveal the nature of the original model [23].
The QoI-based Sobol' indices quantify the influence of the input variables on some quantity of interest of the random model response, e.g., mean, variance, and quantiles [23]. As the mean function in Eq. (26) is a PCE, the associated Sobol' indices can be computed in a straightforward way [36]. Similar to Eq. (36), the first-order index is given by while higher-order and total Sobol' indices of the mean function read In addition, the variance function in Eq. (35) is a polynomial. The associated Sobol' indices can be computed by building another PCE to represent Eq. (35) the without error.

Numerical examples
In this section, we validate the proposed method on several examples, namely case studies from mathematical finance and epidemiology and a complex analytical example with bimodal response distributions. To illustrate its performance, we compare the results obtained from the stochastic polynomial chaos expansion (SPCE) with two state-of-the-art models that are developed for emulating the response distribution of stochastic simulators. The first one is the generalized lambda model (GLaM). This surrogate uses the four-parameter generalized lambda distribution to approximate the response distribution of Y x for any x ∈ D X . The distribution parameters, as functions of the inputs, are represented by PCEs (see details in [13,22]). The second model is based on kernel conditional density estimation (KCDE) [44]. This method uses kernel density estimation to fit the joint distributionf X,Y (x, y) and the marginal distributionf X (x). The response distribution is then estimated by where K y and K j 's are the kernels for Y and X j 's, and h y and h j 's are the associated bandwidths which are hyperparameters selected by a thorough leave-one-out cross-validation [20].
Finally, we also consider a model where we represent the response with a normal distribution.
The associated mean and variance as functions of the input x are set to the true values obtained from the simulator. Therefore, the accuracy of such an approximation measures how close the response distribution is to the normal distribution. Moreover, this model represents the "oracle" of Gaussian-type mean-variance models, such as the ones presented in Marrel et al. [6] and Binois et al. [8].
To quantitatively compare the various surrogates, we define an error metric between the simulator and the emulator by where Y x is the model response,Ỹ x denotes that of the surrogate (with the same input parameters as Y x ), and Y is the model output aggregating all the uncertainties from both the input and the intrinsic stochasticity. d WS is the Wasserstein distance of order two [45] between the two probability distributions defined by where Q 1 and Q 2 are the quantile functions of random variables Y 1 and Y 2 , respectively. The error metric ε in Eq. (41) is unitless and invariant to shift and scale, i.e., To evaluate the numerator in Eq. (41), we generate a test set X test of size N test = 1,000 from the input distribution of X. The Wasserstein distance is calculated for each point x ∈ X test and then averaged over N test .
We use Latin hypercube sampling (LHS) [46] to generate the experimental design and the test set. The stochastic simulator is evaluated only once for each set of input parameters, i.e., we do not use replications. To study the convergence property of the surrogates, experimental designs of various sizes are investigated. Each scenario is run 20 times with independent experimental designs to account for the statistical uncertainty in the LHS design and also in the internal stochasticity of the simulator. As a result, error estimates for each size of experimental design are represented by box plots constructed from the 20 repetitions of the full analysis.

Geometric Brownian motion
In the first example, we consider the Black-Scholes model that is popular in mathematical Equation (44) is a stochastic differential equation used to model the evolution of a stock price S t . Here, x = (x 1 , x 2 ) T are the input variables that describe the expected return rate and the volatility of the stock, respectively. W t is a Wiener process that represents the stochastic behavior of the market. Without loss of generality, we set the initial condition to S 0 = 1.
The simulator is stochastic: for a given x, the stock price S t is a stochastic process, where the stochasticity comes from W t . In this example, we are interested in Y x = S 1 , which corresponds to the stock value at t = 1 year. We set X 1 ∼ U(0, 0.1) and X 2 ∼ U(0.1, 0.4) to represent the uncertainty in the return rate and the volatility, where the ranges are selected based on real data [47].
The solution to Eq. (44) can be derived using Itô calculus [48]: Y x follows a lognormal distribution defined by As the distribution of Y x is known analytically in this simple example, we can sample directly from the response distribution to get the model output instead of simulating the whole path of GLaM accurately represent the variation of the response PDF. Moreover, SPCE better represents the left tail in Fig. 3d. In contrast, KCDE can well approximate the response PDF for low volatility (in Fig. 3a) but exhibits unrealistic oscillations in the case of high volatility. For convergence studies, we vary the size of the experimental design N ∈ {100; 200; 400; 800; 1,600} and plot the error ε defined in Eq. (41) with respect to N in Fig. 4. In order to show more details, each subfigure in Fig. 4 compares SPCE with one competitor. We observe that the average error of KCDE built on N = 400 model runs is similar to the best normal approximation, whereas both SPCE and GLaM provide smaller errors. Compared with KCDE and GLaM, the average performance of SPCE is always the best for all sizes of experimental design. For large N , namely N = 1,600, the average error of SPCE is less than half of that of KCDE, and the spread of the error is narrower than that obtained by GLaM.

Stochastic SIR model
The Because no newborn or death is considered, the three quantities satisfy E t + I t + R t = P . As a result, any two out of the three counts, e.g., E t and I t , can characterize the configuration of the population of size P at time t.
where β is the contact rate of an infected individual, and γ is the recovery rate. The next configuration of the population is the one that comes first, i.e., for T R < T I , the system evolves to C R at t + T R with S t+T R = E t − 1 and I t+T I = I t + 1, and vice versa. We iterates this updating procedure until the time T where I T = 0 corresponding to no remaining infected individual: no infection or recovery can happen, and the outbreak stops. Since the population size is constant and recovered individuals will not be infected again, the outbreak will stop at finite time, i.e., T < +∞. The simulation process described here corresponds to the Gillespie algorithm [49].
The input variables of the simulator are the initial conditions S 0 and I 0 and the transitive rates β and γ. We are interested in the total number of newly infected individuals during the outbreak without counting the initial infections, which is an important quantity in epidemics management [8]. This can be calculated by the difference between the number of susceptibles at time 0 and T , i.e., Y = S 0 − S T . Because each updating step in Eq. (46) depends on two latent variables T I and T R , the simulator is stochastic. Moreover, the total number of latent variables is also random.
In this case study, we set P = 2,000. To account for different scenarios, the input variables The uncertainty in the first two variables is due to the lack of knowledge of the initial condition.
The two transitive rates β, γ are affected by possible interventions such as quarantine and increase of medical resources. Figure 6 illustrates the response PDF for four different sets of input parameters. Because of the transition process in Eq. (46), no analytical closed-form distribution of Y x can be derived.
Therefore, we use 10 4 replications for each input values to obtain the reference histograms.
The surrogate models are trained on an experimental design of size N = 1,600 (without any replications). We observe that the four PDFs are unimodal. The reference histogram in Fig. 6a is slightly right-skewed, while the others in Fig. 6 are symmetric. SPCE and GLaM produce similar predictions of the PDF which are very close to the reference histograms. In comparison, KCDE overestimates the spread of the distributions in. Moreover, the KCDE prediction has non-negligible probability for unrealistic negative values in Fig. 6a. Besides, it exhibits relatively poor shape representations with spurious wiggles in Fig. 6c and Fig. 6d.  normal approximation gives an error of 6 × 10 −4 which is smaller than any of the surrogates in consideration. Note that this model is not built on the training data but using the mean and variance from the 10 4 replications for each test point. This implies that the response distribution is close to normal. We do not include this error in Fig. 7 to not loose detailed comparisons of the surrogate models. Figure 7 reveals a poor performance of KCDE in this case study. This is because the example is four-dimensional, and KCDE is a kernel-based method which is known to suffer from the curse of dimensionality. In contrast, SPCE and GLaM are flexible parametric models, and both provide a much smaller error than KCDE for all values of N . Compared with GLaM, SPCE yields a similar spread of the error but demonstrates better average performance for N ≥ 400.

Bimodal analytical example
The response distributions of the previous two examples are unimodal. In the last example, we consider a complex analytical example to test the flexibility of the stochastic polynomial chaos expansion. For this purpose, we directly define the response distribution to approximate as where ϕ stands for the standard normal PDF. This response PDF is a mixture of two Gaussian  We observe that small values of x yield a bimodal distribution with the higher mode on the right. With x increasing, the two modes merge and form a unimodal distribution at x = 0.5. Then, the two modes separate again, which leads to bimodal distributions with the higher mode on the left. This shape variation can also be observed from Fig. 8.
As opposed to the previous two examples, GLaM cannot represent this evolution, since generalized lambda distributions cannot produce multimodal distributions. In contrast, SPCE and KCDE capture well the bimodality and also the shape variation. Moreover, in Fig. 9c the higher mode is moving to the left, which is a feature not exhibited by KCDE but correctly captured by SPCE.

Conclusions
In this paper, we present a novel surrogate model called stochastic polynomial chaos expansions (SPCE) to emulate the response distribution of stochastic simulators. This surrogate is an extension of the classical polynomial chaos expansions developed for deterministic simulators. In order to represent the intrinsic stochasticity of the simulator, we combine a latent variable with the well-defined inputs to form a polynomial chaos representation. In addition, we introduce an additive Gaussian noise as a regularizer. We propose using the maximum likelihood estimation for calibrating the coefficients c of the polynomial basis. The standard deviation σ of the noise variable is a hyperparameter that regularizes the optimization problem for the polynomial coefficients c and is tuned by cross-validation to avoid overfitting. The cross-validation score is also used as a model selection criterion to choose an appropriate truncation scheme for the polynomial chaos expansion in an adaptive manner, and the most suitable distribution for the latent variable. As seen from the presentation and the application examples, the proposed method does not require replications.
The performance of the developed method is illustrated on examples from mathematical finance and epidemiology and on an analytical example showcasing a bimodal response distribution.
The results show that SPCE is able to well approximate various response distributions whether unimodal or not, with a reasonable number of model runs.
Using an appropriate error measure defined in Eq. (41), SPCE is compared with the generalized lambda model (GLaM) and one state-of-the-art kernel conditional density estimator (KCDE). In the first two examples where the response distribution is unimodal, SPCE noticeably outperforms KCDE and provides slightly more accurate results than GLaM which is known for its flexibility for representing unimodal distributions. In the last example featuring bimodal distributions which cannot be well approximated by generalized lambda distributions, SPCE can still capture the complex shape variation and yields smaller errors than KCDE. All in all, SPCE generally performs as the best against the various competitors considered in this study.
Applications of the proposed method to complex engineering problems, such as wind turbine design [50] and structural dynamics [51], should be considered in future investigations. Statistical properties (e.g., consistency and asymptotics) of the maximum likelihood estimation used in SPCE remains to be studied. This will allow for assessing the uncertainty in the estimation procedure.
Finally, the proposed approach has been validated so far only for problems with small to moderate dimentionality. To improve the efficiency and performance of SPCE in high dimensions, models that have a general sparse structure (not only regarding the mean function) are currently under investigations.

Acknowledgments
This paper is a part of the project "Surrogate Modeling for Stochastic Simulators (SAMOS)" funded by the Swiss National Science Foundation (Grant #200021_175524), whose support is gratefully acknowledged.

A.1 Upper bound
In this section, we demonstrate that the leave-one-out error obtained from fitting the mean function Eq. (28) provides an upper bound for σ 2 .
Taking the expectation of Eq. (35) with respect to X, it holds The leave-one-out error ε LOO in the mean-fitting process is an estimate of E (m(X) − Y X ) 2 [52]. The latter can be decomposed as