Multi Fidelity Moving Particles Algorithm Considering Aleatory and Epistemic Uncertainty

Multi Fidelity Moving Particles Algorithm Considering Aleatory and Epistemic Uncertainty⋆ Jonas Kaupp1, Carsten Proppe1 Institute of Engineering Mechanics, Karlsruhe Institute of Technology, Kaiserstr. 10, Bdg. 10.23, 76131 Karlsruhe, Germany, {jonas.kaupp, proppe}@kit.edu Abstract. Calculating the probability of failure in structural reliability analysis often faces high computational costs. The costs increase further if imprecise probabilities are present, since an optimization problem must be solved for each evaluation of the model. In this paper, a multi fidelity algorithm based on the moving particles method is proposed to reduce the number of cost intensive model evaluations needed. It uses either additive or multiplicative information fusion to combine the different model fidelity outputs. Furthermore, probability boxes are used to represent the imprecise probabilities. The resulting multi fidelity estimator is tested on two examples. It is shown that a reduction of relative error and variance is achieved, but depends on the posed problem. This results in fewer evaluations of the high fidelity model if the same relative error is used as a termination threshold. Keywords: reliability estimation, information fusion, multi fidelity, moving particles algorithm, imprecise probabilities, probability boxes 1 Introduction In structural reliability analysis the probability of failure pfis evaluated with respect to a quantity of interest (e.g. displacement or stress). For given random input parameters xof dimension dwith joint probability density function p(x) and input domain D ∈ Rd, the failure probability is calculated as pf=ZFF p(x)dx(1) with failure domain FF={x∈ D|g(x)<0}. The limit state function (LSF) g(x) becomes negative if the model response for the input parameters is greater than a given threshold for the quantity of interest. Usually a Monte Carlo Simulation (MCS) is carried out, since the integral over FFis not analytically tractable or the LSF defining the bounds of the failure domain is highly nonlinear. MCS is independent of the dimension of the input domain, but converges rather slowly, ⋆DFG SPP 1886 subproject 24 2 Kaupp, Proppe especially if pfis small. [17] Furthermore, the evaluation of the LSF can be computationally intensive. Therefore, advanced sampling techniques, also known as variance reduction techniques, which evaluate less samples for the same estimator accuracy, have been introduced. For importance sampling, samples are drawn from a better suited distribution to cover the failure domain and those samples are re-weighted afterwards such that the bias from sampling a different distribution vanishes. [18] Additionally, there are so called importance splitting techniques which calculate the probability of failure by using intermediate failure domains which create conditional distributions. From these, samples are drawn via Markov Chain Monte Carlo simulation (MCMC). To compute pf, all conditional probilities are multiplied. This leads to a computational advantage since the conditional probabilities are large compared to pfand thus it is possible to estimate them with a low coefficient of variation from rather few samples. Importance splitting methods comprise subset simulation with a fixed intermediate failure threshold (e.g. p0= 0.1) [3] and the moving particle algorithm with a maximum number of intermediate failure domains [11]. Another way to improve the estimator is by combining multiple models of either varying levels with discretization parameter h, also known as multi level estimators, introduced by [9] and [12], or of different fidelities, where no restrictions on the model types are made [20]. In the latter approach surrogate models can be implemented easily, however the correlation of the models has to be high in order to achieve an efficient estimator. [22] introduced a multi level version of the moving particles method leading to a considerable reduction of computational cost. The majority of the scientific literature considers aleatory uncertainties, but there exist uncertainties caused by not knowing the underlying system and its parameters completely, referred to as epistemic uncertainties. The latter can be reduced with more effort, e.g. more measurements, and therefore higher costs. [6] Another term for the combination of aleatory and epistemic uncertainties is imprecise probability. These imprecise probabilities can be captured using intervals [29], fuzzy variables [27], probability boxes (p-boxes) [8] or combinations thereof. Importance sampling has also been extended to capture intervals [28]. [2] applied importance splitting for discretized p-boxes. [24] presented a combination of surrogate models and p-boxes to compute the probability of failure. Despite these advances, a multi fidelity approach including imprecise probabilites does not exist. The moving particles algorithm is suitable in particular since it is easily adaptable for a multi fidelity setting and has a sound mathematical foundation rooted in Poisson process theory. The paper is structured as follows: at first, the applied methods are explained in detail in Section 2. This includes the moving particles algorithm, multi fidelity estimators and p-boxes. Second, in Section 3 the multi fidelity moving particles estimator (MFMP) is introduced and its adaption to p-boxes is described. Third, the estimator is applied to a simple academic example and a single de- MFMP Considering A. & E. Uncertainty 3 gree of freedom (SDOF) oscillator in Section 4. Lastly, conclusions are drawn and future research is pointed out in Section 5. 2 Methods Starting with the moving particles algorithm the fundamental concepts and methods are introduced. Herein, additive and multiplicative multi fidelity estimators are presented, afterwards the expectation maximization (EM) algorithm for a multivariate Poisson distribution is explained and lastly, imprecise probabilities in the form of p-boxes are described. 2.1 Moving Particles Algorithm The moving particles algorithm gives an estimate of the probability of failure pf. The algorithm starts with an initial MCS in the input domain D ∈ Rdwith a prior joint distribution. dis the number of input dimensions. Each sample xi∈ D, i = 1, ..., N – called particle in this section – and its LSF evaluation g(xi) are sorted such that the one with the highest value gmax (farthest away from the failure domain FF) is taken as the first one to move. Now the aim is to reduce the value of the LSF, such that g(xi,new)< g(xi). All other particles currently fulfill this criterion and are part of a conditional distribution resulting from the prior distribution with condition g(xi)< gmax, which needs to be sampled from. This is done by performing MCMC with one of the particles in the remaining set being the starting point of the Markov Chain, which is generated e.g. by the Metropolis-Hastings Kernel or improvements thereof. If the new state xi,new is accepted, the particle is moved to its new location and xi=xi,new. Counting the number of moves Miof each particle until it reaches FFfollows a Poisson distribution with parameter λ=PN i=1 Mi N(2) for a total number of Nparticles. Lastly, the probability of failure is linked to λ via pf= exp(−λ). [11,26] The independence of the generated Poisson process trajectories must be ensured. This is achieved by introducing seed avoidance [26] and a burn-in period for the Markov Chain, albeit a rather small one [22]. On average, the number of times the LSF is evaluated is Ntot =N(1 −Nblog(pf)) (3) with Nbbeing the burn-in period and the estimator’s coefficient of variation results in δ=r−log(pf) N(4) 4 Kaupp, Proppe 2.2 Multi Fidelity Estimators Multi fidelity estimators are obtained by additive or multiplicative information fusion such that the output of different models is combined in an efficient manner. In the first approach, the output of one model is used and refined by adding evaluations of different fidelity models. [23] [20] used the control variate framework to additively combine the models and distributed the number of model evaluations optimally with the following formula: EMMF=EMhi+ L−1 X l=1 αlEMlo l,Nl−EMlo l,Nl+1  (5) Here, Mhi denotes the output of the high fidelity model, Mlo l,Nlthe output of the l-th low fidelity model for a total of Lmodels with ML=Mhi. For the determination of the number of evaluations of each model Nland the weighting factor αlthe interested reader is referred to [20] or the review of multi fidelity estimators published in [21]. For multiplicative fusion either a Bayesian regression can be performed and a surrogate model is used to link the output of the different models in a (non) linear way [4] or expectation maximization is applied [23]. As the distribution of the moves is Poisson, model dependencies in the moving particles method are taken into account by a multivariate Poisson distribution with parameters λl, l = 1, ..., L with Lbeing the number of models. The algorithm is described in the next section. 2.3 Expectation Maximization for multivariate Poisson distributions First, the structure of the multivariate Poisson distribution and its decomposition into independent univariate Poisson distributions is introduced. Afterwards, the expectation maximization algorithm for multivariate Poisson distributions as a combination of the ideas of [15] and [1] is presented. Multivariate Poisson Distribution Generally, a random variable X∈Rk following a multivariate Poisson distribution can be described by independent, univariate Poisson variables Ym∼Po(θm), m = 1, ..., M with a link matrix Φ∈Rk×Mand M= 2k−1. X=ΦY (6) The expected value and the variance result in E(X) = Φθ(7) V ar(X) = ΦΣΦT(8) with Σ= diag(θ1, ..., θm), because of the independence of Ym. [16,15] MFMP Considering A. & E. Uncertainty 5 Expectation Maximization In general, the expectation maximization algorithm finds maximum likelihood estimates of the parameters of probability distributions for a given number of realizations iteratively. It is divided into two main steps starting with the expectation step (E-step), where the expectation of the log-likelihood function is computed depending on the current distribution parameters and realizations. Afterwards, in the maximization step (M-step) the log-likelihood is maximized with respect to the distribution parameters. The algorithm alternates between those steps and terminates after the relative change of the parameters passes a given threshold. [5] The EM algorithm for multivariate Poisson distributions analogously finds the most likely parameters θmof the univariate Poisson distributions Ymfollowing the above mentioned decomposition for a drawn number of realizations. Afterwards, the parameters of the multivariate distribution are obtained using the link matrix Φ. Given Zsamples of X, being the multivariate Poisson distribution, a pseudo value sis derived in the E-step which represents the expected value of the loglikelihood of each independent Poisson distribution in the presence of the sample xi[15]: sim =E(Ym|xi,θt−1) = θmP(X=xi−ϕm) P(X=xi).(9) Here, the index i= 1, ..., Z denotes the sample xiused and mrefers to the independent univariate Poisson distributed random variable Ymwith parameter θmat the (t−1)-th step of the EM algorithm. P(·) are the probabilities of the multivariate Poisson distribution and ϕmthe corresponding column vector of Φ. The probabilities can be computed from a recurrence scheme [14]. Note that the probability of the numerator is readily available when computing the denominator. For the M-step, the ideas of [15] and [1] are combined to compute the expected values of the pseudo values svia E(sm|X,θt−1) = Z X i=1 sim = Z X i=1 E(Ym|xi,θt−1).(10) Finally, the new parameters of the independent Poisson distributions for step t result in θt m=1 ZE(sm|X,θt−1).(11) The Eand M-step are repeated until ε<εth with εbeing M X m=1 |θt m−θt−1 m| θt−1 m .(12) A suitable value for εth might be 1 ·10−3. 2.4 Probability Boxes P-boxes can be used to describe imprecise probabilities. Instead of using one cumulative distribution function F(CDF) which is known or estimated, an upper 6 Kaupp, Proppe and lower bound on the CDF is introduced, Fand Frespectively (see Fig 1). The true CDF is not known. [8,7] P-boxes are classified as distributional or distribution-free. A distributional pbox has a known distribution family with interval valued parameters. The distribution-free p-box has no restrictions and the bounds are inferred from data. Both can be discretized and treated as random sets afterwards, but the information on the family is lost in the case of distributional p-boxes [2]. Figure 1 shows a distributional or parametric p-box on the left and a distribution-free p-box on the right. For the former, a normal family with interval mean and standard deviation is used (X∼ N([−1,1],[1,2])) and for the latter only minimum (min), maximum (max), mean and standard deviation (std) are known, with the tuple (min, max, mean, std) being (-4, 4, 0, 2). Here, Fis in red and Fin black. The toolbox by [10] is used to calculate the p-boxes and create figures thereof. If d > 1 for the input domain, the dimensions can be treated as independent or coupled. To couple the dimensions, a copula, also called a dependence structure, is implemented. There are a plethora of different types of copulae. For strictly aleatory uncertainties the copula and the marginal distributions form the joint probability density function (Sklar’s Theorem). The copula with standard uniform marginals Uconnects the physical domains xivia (x1, x2, ..., xd) = (F−1 1(u1), F−1 2(u2), ..., F −1 d(ud)) with ui∈ U. [19] To include p-boxes the principle remains the same, but the quantiles uilead to random sets [F(ui), F (ui)] [2]. −7.5 −5.0 −2.5 0.0 2.5 5.0 7.5 0.0 0.2 0.4 0.6 0.8 1.0 −4 −3 −2 −1 0 1 2 3 4 0.0 0.2 0.4 0.6 0.8 1.0 Fig. 1. a distributional p-box of the normal family is seen on the left, a distribution-free p-box with only known minimum, maximum, mean and standard deviation is shown on the right 3 Multi Fidelity Moving Particles Algorithm As mentioned in Section 2.1 the moving particles algorithm leads to a Poisson process with parameter λfor a single fidelity (SF) problem. Applying the same algorithm to a multi fidelity setting needs adaptions. MFMP Considering A. & E. Uncertainty 7 3.1 Data structure At the lowest fidelity the SF algorithm is performed. With regards to the computational costs, the model is usually cheap. The MCMC resulting from the lower fidelity model must be fed to the next higher fidelity model at each move. This is achieved using the scheme by [13], where every generated random number of the Markov Chain is stored. The new particle location is accepted if the LSF value of the current model level is reduced. Since the different models are interdependent, the number of moves of the particles are realisations of a multivariate Poisson distribution. Another goal is the reduction of particles for the higher fidelity models. As described in Section 2.2, the lowest fidelity model has N0=Nparticles, the higher fidelity models have Nlparticles respectively and Nl≥Nl+1. This leads to an upper triangular structure of particles if Nl> Nl+1 (see Fig. 2). The structure is then split into tuple of particles existent for all models. Each tuple forms a category cwith category 1 containing NLparticles of every model. The subsequent categories consist of Nl−Nl+1 particles with less and less columns until the last category Cwith only lowest fidelity particles is reached. Overall, it must be ensured that enough particles are in each category to maintain independence. The different fidelities can be combined using additive or multiplicative information fusion which is shown in the next sections. Lastly, the algorithm is extended to account for imprecise probabilities. Fig. 2. The particles of given fidelities lead to an upper triangular particle structure. A category is formed if particles are available for all models. 8 Kaupp, Proppe 3.2 Additive Fusion For additive fusion the multi level version of the moving particles algorithm developed by [22] E(ML) = E(M0) + L X l=1 E(Ml−Ml−1) =1 N0 N0 X i=1 M(i) 0+ L X l=1 1 Nl Nl X i=1 M(i) l−M(i) l−1(13) can be combined with the multi fidelity estimator given by [20]. The resulting estimator reads ˆ ML=1 NL NL X i=1 M(i) L+ L−1 X l=1 αl  1 Nl Nl X i=1 M(i) l−1 Nl+1 Nl+1 X i=1 M(i) l (14) =E(ML) + L−1 X l=1 αl(E(Ml)−E(Ml+1)) with αl=ρL,lσL σl (15) being the weighting factor which is calculated from a few pilot samples. ρL,l is the Pearson correlation coefficient of model lto high fidelity model Land σis the standard deviation of the corresponding model. Furthermore, M(i) ldenotes the number of moves of the i-th particle for the l-th model and Nlthe total number of particles of the l-th model. 3.3 Multiplicative Fusion To properly apply the EM algorithm presented in section 2.3 a complete set of samples needs to exist in order to estimate the parameters for the next iteration. Since this constraint cannot be satisfied the parameters are computed in a different way, which was suggested by [1] for the bivariate Poisson distribution and is extended here to a multivariate Poisson distribution. To make the transition more comprehensible, the bivariate version of [1] is stated first and then extended to the general multivariate case. Here, the expected values of the pseudo values sare calculated as follows: E(sm|X,θt−1) = Z1 X i=1 E(Ym|xcomplete,i,θt−1) + Z2 X i=1 E(Ym|x2,i,θt−1) + Z3 X i=1 E(Ym|x1,i,θt−1) (16) MFMP Considering A. & E. Uncertainty 9 For a total of Zsamples in X, the samples are divided into three categories: complete sample, missing X1value or missing X2value. The number of samples in each category is denoted with Z1,Z2and Z3respectively. If applied to a bivariate moving particles algorithm, the last term vanishes, since the adaptions lead to a structure, where the higher fidelity model has less evaluations than the model with lower fidelity. Using this information, the generalized calculation of the expected values for a multivariate problem simplifies to: E(sm|X,θt−1) = Z1 X i=1 E(Ym|xc=1,l=1,i, ..., xc=1,l=L,i,θt−1) + Z2 X i=1 E(Ym|xc=2,l=1,i, ..., xc=2,l=L−1,i,θt−1) +... + ZC X i=1 E(Ym|xc=C,l=1,i,θt−1) (17) with cfor the category and lfor the model. The number of samples in each category is Zc=Nl−Nl+1 as explained in section 3.1. The parameters θt mcan then be determined from Eq. 11. Since this is the general multivariate case, the interested reader can find the computations necessary for the expected values of the trivariate case in appendix A. 3.4 Extension of the algorithm to p-boxes P-boxes can be implemented easily using the approach proposed by [2]. Here, the p-boxes are discretized and treated as random sets. In the physical input space these sets form hypercubes, also called focal elements. These can be transformed into the copula space, where they equal points again. For these points, the MCMC for the moving particles algorithm can be performed as usual. It only differs in the acceptance/rejection of the proposal state. Here, a transformation back to the physical space and the solution of an optimization problem is needed to get the upper and lower bounds of the LSF, gand grespectively. As g > 0 implies g > 0, gmust only be evaluated if g < 0. At the end, λand λare obtained which finally result in the upper and lower bounds for the probability of failure pfand pf. If a multi fidelity problem is faced, the implementation of an additive or multiplicative version is needed and the bounds of pfcan be estimated separately. 3.5 Construction of the estimator The estimator is divided into two parts. It starts with the MFMP which provides the moves of each particle for each fidelity. In Appendix B the resulting MFMP estimator is sketched. At first, a pilot run is done to get the optimal samples and correlations per fidelity according to [20]. Second, the categories are looped, 16 Kaupp, Proppe of further research to decrease the computation time of the recursive algorithm and to be a competition to the additive information fusion. Acknowledgements This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project no. 428469142 – as part of the priority programme ”Polymorphic uncertainty modelling for the numerical design of structures” (SPP 1886) which is gratefully acknowledged by the authors. CRediT author statement Jonas Kaupp: Conceptualisation, Methodology, Software, Validation, Formal Analysis, Investigation, Visualisation, Writing - Original Draft - Review & Editing Carsten Proppe: Funding acquisition, Project administration, Supervision, Writing - Review & Editing A Expected Values for the determination of independent univariate Poisson distributions Recall that a trivariate Poisson distribution can be written as X=ΦY ⇐⇒   X1 X2 X3 =  1000111 0101011 0011101            Y1 Y2 Y3 Y4 Y5 Y6 Y7           (26) If now an incomplete sample is present, the expected values for Yiresult in E(Y1|X2,i, X3,i) =E(Y1|X3,i) = E(Y1) = λ1 E(Y2|X2,i, X3,i) = X2,i X j=0 j·P(Y2=j|Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X2,i X j=0 j·P(Y2=j)· P(Y4+Y6+Y7=X2,i −j∩Y3+Y4+Y5+Y7=X3,i) P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X2,i X j=0 j·Po(Y2=j)BP o(Z1+Z3=X2,i −j, Z2+Z3=X3,i) BPo(Z4+Z6=X2,i, Z5+Z6=X3,i) MFMP Considering A. & E. Uncertainty 17 with Po(·) and BP o(·) being a univariate and bivariate Poisson distribution, respectively. λiis the parameter to describe the independent univariate Poisson distribution Yi. The remaining parameters are given by: Z1=Y6, Z2=Z5=Y3+Y5, Z3=Z6=Y4+Y7, Z4=Y2+Y6. Note, Zis not the number of samples as in Eq. 17, but a random variable. The other expected values are calculated analogously: E(Y2|X3,i) =E(Y2) = λ2 E(Y3|X2,i, X3,i) = X3,i X j=0 j·P(Y3=j|Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X3,i X j=0 j·P(Y3=j)· P(Y2+Y4+Y6+Y7=X2,i ∩Y4+Y5+Y7=X3,i −j) P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X3,i X j=0 j·Po(Y3=j)BP o(Z1+Z3=X2,i, Z2+Z3=X3,i −j) BPo(Z4+Z6=X2,i, Z5+Z6=X3,i) with: Z1=Z4=Y2+Y6, Z2=Y5, Z3=Z6=Y4+Y7, Z5=Y3+Y5 E(Y3|X3,i) = EBX3,i,λ3 λ3+λ4+λ5+λ7=X3,i λ3 λ3+λ4+λ5+λ7 Note, here B(·) is equal to a Bernoulli distribution. Furthermore: E(Y4|X2,i, X3,i) = min(X2,i,X3,i) X j=0 j·P(Y4=j|Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = min(X2,i,X3,i) X j=0 j·P(Y4=j)· P(Y2+Y6+Y7=X2,i −j∩Y3+Y5+Y7=X3,i −j) P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = min(X2,i,X3,i) X j=0 j·Po(Y4=j)BP o(Z1+Z3=X2,i −j, Z2+Z3=X3,i −j) BPo(Z4+Z6=X2,i, Z5+Z6=X3,i) with: Z1=Z4=Y2+Y6, Z2=Z5=Y3+Y5, Z3=Y7, Z6=Y3+Y4 E(Y4|X3,i) =EBX3,i,λ4 λ3+λ4+λ5+λ7=X3,i λ4 λ3+λ4+λ5+λ7 18 Kaupp, Proppe E(Y5|X2,i, X3,i) = X3,i X j=0 j·P(Y5=j|Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X3,i X j=0 j·P(Y5=j)· P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y7=X3,i −j) P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X3,i X j=0 j·Po(Y5=j)BP o(Z1+Z3=X2,i −j, Z2+Z3=X3,i −j) BPo(Z4+Z6=X2,i, Z5+Z6=X3,i) with: Z1=Z4=Y2+Y6, Z2=Y3, Z3=Z6=Y4+Y7, Z5=Y3+Y5 E(Y5|X3,i) =EBX3,i,λ5 λ3+λ4+λ5+λ7=X3,i λ5 λ3+λ4+λ5+λ7 E(Y6|X2,i, X3,i) = X2,i X j=0 j·P(Y6=j|Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X2,i X j=0 j·P(Y6=j)· P(Y2+Y4+Y7=X2,i −j∩Y3+Y4+Y5+Y7=X3,i) P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = X2,i X j=0 j·Po(Y6=j)BP o(Z1+Z3=X2,i −j, Z2+Z3=X3,i) BPo(Z4+Z6=X2,i, Z5+Z6=X3,i) with: Z1=Y2, Z2=Z5=Y3+Y5, Z3=Z6=Y4+Y7, Z4=Y2+Y6 E(Y6|X3,i) =E(Y6) = λ6 E(Y7|X2,i, X3,i) = min(X2,i,X3,i) X j=0 j·P(Y7=j|Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = min(X2,i,X3,i) X j=0 j·P(Y7=j)· P(Y2+Y4+Y6=X2,i −j∩Y3+Y4+Y5=X3,i −j) P(Y2+Y4+Y6+Y7=X2,i ∩Y3+Y4+Y5+Y7=X3,i) = min(X2,i,X3,i) X j=0 j·Po(Y7=j)BP o(Z1+Z3=X2,i −j, Z2+Z3=X3,i −j) BPo(Z4+Z6=X2,i, Z5+Z6=X3,i) with: Z1=Z4=Y2+Y6, Z2=Z5=Y3+Y5, Z3=Y4, Z6=Y4+Y7 E(Y7|X3,i) =EBX3,i,λ7 λ3+λ4+λ5+λ7=X3,i λ7 λ3+λ4+λ5+λ7 MFMP Considering A. & E. Uncertainty 19 B Flow chart of the estimator Start Models Distributions/P-boxes Moving Particles Options First run? Initial interval MCS for every model c= 1 Pilot run for ρL,l, Nl, αland used models with the same number of particles (see [20]) c > C? Calculate gmax c=c+ 1 Calculate gmax evaluate moves per model concatenate categories to form Fig. 2 calculate λ, λ via Eq. 14 or EM separately and pf= exp(−λ), pf= exp(−λ) Stop false true false true Fig. 10. Main scheme of the estimator 20 Kaupp, Proppe References 1. Adamids, K., Loukas, S.: Ml estimation in the bivariate poisson distribution in the presence of missing values via the em algorithm. Journal of Statistical Computation and Simulation 50(3-4), 163–172 (1994). https://doi.org/10.1080/00949659408811608 2. Alvarez, D.A., Uribe, F., Hurtado, J.E.: Estimation of the lower and upper bounds on the probability of failure using subset simulation and random set theory. Mechanical Systems and Signal Processing 100, 782–801 (2018). https://doi.org/10.1016/j.ymssp.2017.07.040 3. Au, S.K., Beck, J.L.: Estimation of small failure probabilities in high dimensions by subset simulation. Probabilistic Engineering Mechanics 16(4), 263–277 (2001). https://doi.org/10.1016/S0266-8920(01)00019-4 4. Biehler, J., Gee, M.W., Wall, W.A.: Towards efficient uncertainty quantification in complex and large-scale biomechanical problems based on a bayesian multi-fidelity scheme. Biomechanics and modeling in mechanobiology 14(3), 489–513 (2015). https://doi.org/10.1007/s10237-014-0618-0 5. Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society. Series B (Methodological) 39(1), 1–38 (1977), http://www.jstor.org/stable/2984875 6. Der Kiureghian, A., Ditlevsen, O.: Aleatory or epistemic? does it matter? Structural Safety 31(2), 105–112 (2009). https://doi.org/10.1016/j.strusafe.2008.06.020 7. Faes, M.G., Daub, M., Marelli, S., Patelli, E., Beer, M.: Engineering analysis with probability boxes: A review on computational methods. Structural Safety 93, 102092 (2021). https://doi.org/10.1016/j.strusafe.2021.102092 8. Ferson, S., Kreinovich, V., Ginzburg, L., Myers, D.S., Sentz, K.: Constructing probability boxes and dempster-shafer structures (2003), https://citeseerx.ist.psu. edu/viewdoc/download?doi=10.1.1.3.4623&rep=rep1&type=pdf 9. Giles, M.B.: Multilevel monte carlo path simulation. Operations Research 56(3), 607–617 (2008). https://doi.org/10.1287/opre.1070.0496 10. Gray, N., Ferson, S., de Angelis, M., Gray, A., Baumont de Oliveira, F.: Probability bounds analysis for python. Software Impacts 12, 100246 (2022). https://doi.org/10.1016/j.simpa.2022.100246 11. Guyader, A., Hengartner, N., Matzner-Løber, E.: Simulation and estimation of extreme quantiles and extreme probabilities. Applied Mathematics & Optimization 64(2), 171–196 (2011). https://doi.org/10.1007/s00245-011-9135-z 12. Heinrich, S.: Multilevel monte carlo methods. In: Margenov, S., Wa´sniewski, J., Yalamov, P. (eds.) Large-Scale Scientific Computing, Lecture Notes in Computer Science, vol. 2179, pp. 58–67. Springer Berlin Heidelberg, Berlin, Heidelberg (2001). https://doi.org/10.1007/3-540-45346-6 5 13. Jasra, A., Law, K.J.H., Xu, Y.: Markov chain simulation for multilevel monte carlo. Foundations of Data Science 3(1), 27 (2021). https://doi.org/10.3934/fods.2021004 14. Kano, K., Kawamura, K.: On recurrence relations for the probability function of multivariate generalized poisson distribution. Communications in Statistics - Theory and Methods 20(1), 165–178 (1991). https://doi.org/10.1080/03610929108830490 15. Karlis, D., Meligkotsidou, L.: Multivariate poisson regression with covariance structure. Statistics and Computing 15(4), 255–265 (2005). https://doi.org/10.1007/s11222-005-4069-4 MFMP Considering A. & E. Uncertainty 21 16. Kawamura, K.: The structure of multivariate poisson distribution. Kodai Mathematical Journal 2(3) (1979). https://doi.org/10.2996/kmj/1138036064 17. Lemaire, M.: Structural reliability. ISTE and John Wiley & Sons, London and Hoboken, New Jersey (2010). https://doi.org/10.1002/9780470611708, http://site. ebrary.com/lib/alltitles/docDetail.action?docID=10361215 18. Melchers, R.E., Beck, A.T.: Structural reliability analysis and prediction. Wiley, Hoboken, NJ, third edition edn. (2018), http://search.ebscohost.com/login.aspx? direct=true&scope=site&db=nlebk&AN=1619149 19. Nelsen, R.B.: An Introduction to Copulas. Springer Series in Statistics, Springer Science+Business Media Inc, New York, NY, second edition edn. (2006). https://doi.org/10.1007/0-387-28678-0, http://site.ebrary.com/lib/ alltitles/docDetail.action?docID=10187506 20. Peherstorfer, B., Willcox, K., Gunzburger, M.: Optimal model management for multifidelity monte carlo estimation. SIAM Journal on Scientific Computing 38(5), A3163–A3194 (2016). https://doi.org/10.1137/15M1046472 21. Peherstorfer, B., Willcox, K., Gunzburger, M.: Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Review 60(3), 550– 591 (2018). https://doi.org/10.1137/16M1082469 22. Proppe, C.: A multilevel moving particles method for reliability estimation. Probabilistic Engineering Mechanics 59, 103018 (2020). https://doi.org/10.1016/j.probengmech.2020.103018 23. Proppe, C., Kaupp, J.: On information fusion for reliability estimation with multifidelity models. Probabilistic Engineering Mechanics 69, 103291 (2022). https://doi.org/10.1016/j.probengmech.2022.103291 24. Sch¨obi, R., Sudret, B.: Structural reliability analysis for p-boxes using multilevel meta-models. Probabilistic Engineering Mechanics 48, 27–38 (2017). https://doi.org/10.1016/j.probengmech.2017.04.001 25. Schueremans, L., van Gemert, D.: Benefit of splines and neural networks in simulation based structural reliability analysis. Structural Safety 27(3), 246–261 (2005). https://doi.org/10.1016/j.strusafe.2004.11.001 26. Walter, C.: Moving particles: A parallel optimal multilevel splitting method with application in quantiles estimation and meta-model based algorithms. Structural Safety 55, 10–25 (2015). https://doi.org/10.1016/j.strusafe.2015.02.002 27. Zadeh, L.A.: The concept of a linguistic variable and its application to approximate reasoning—i. Information Sciences 8(3), 199–249 (1975). https://doi.org/10.1016/0020-0255(75)90036-5 28. Zhang, H.: Interval importance sampling method for finite element-based structural reliability assessment under parameter uncertainties. Structural Safety 38, 1–10 (2012). https://doi.org/10.1016/j.strusafe.2012.01.003 29. Zhang, H., Mullen, R.L., Muhanna, R.L.: Interval monte carlo methods for structural reliability. Structural Safety 32(3), 183–190 (2010). https://doi.org/10.1016/j.strusafe.2010.01.001