S P R I N G E R B R I E F S I N S TAT I S T I C SLouis J. M. AslettFrank P. A. CoolenJasper De Bock EditorsUncertainty in EngineeringIntroduction to Methods and ApplicationsSpringerBriefs in StatisticsSpringerBriefs present concise summaries of cutting-edge research and practicalapplications across a wide spectrum of fields. Featuring compact volumes of 50 to125 pages, the series covers a range of content fromprofessional to academic. Typicaltopics might include:• A timely report of state-of-the art analytical techniques• A bridge between new research results, as published in journal articles, and acontextual literature review• A snapshot of a hot or emerging topic• An in-depth case study or clinical example• A presentation of core concepts that students must understand in order to makeindependent contributionsSpringerBriefs in Statistics showcase emerging theory, empirical research, andpractical application in Statistics from a global author community.SpringerBriefs are characterized by fast, global electronic dissemination, standardpublishing contracts, standardizedmanuscript preparation and formatting guidelines,and expedited production schedules.More information about this series at https://link.springer.com/bookseries/8921Louis J. M. Aslett · Frank P. A. Coolen ·Jasper De BockEditorsUncertainty in EngineeringIntroduction to Methods and ApplicationsEditorsLouis J. M. AslettDepartment of Mathematical SciencesDurham UniversityDurham, UKJasper De BockFoundations Lab for imprecise probabilitiesGhent UniversityZwijnaarde, BelgiumFrank P. A. CoolenDepartment of Mathematical SciencesDurham UniversityDurham, UKISSN 2191-544X ISSN 2191-5458 (electronic)SpringerBriefs in StatisticsISBN 978-3-030-83639-9 ISBN 978-3-030-83640-5 (eBook)https://doi.org/10.1007/978-3-030-83640-5© The Editor(s) (if applicable) and The Author(s) 2022. This book is an open access publication.OpenAccess This book is licensed under the terms of the Creative Commons Attribution 4.0 InternationalLicense (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribu-tion and reproduction in any medium or format, as long as you give appropriate credit to the originalauthor(s) and the source, provide a link to the Creative Commons license and indicate if changes weremade.The images or other third party material in this book are included in the book’s Creative Commons license,unless indicated otherwise in a credit line to the material. If material is not included in the book’s CreativeCommons license and your intended use is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyright holder.The use of general descriptive names, registered names, trademarks, service marks, etc. in this publicationdoes not imply, even in the absence of a specific statement, that such names are exempt from the relevantprotective laws and regulations and therefore free for general use.The publisher, the authors and the editors are safe to assume that the advice and information in this bookare believed to be true and accurate at the date of publication. Neither the publisher nor the authors orthe editors give a warranty, expressed or implied, with respect to the material contained herein or for anyerrors or omissions that may have been made. The publisher remains neutral with regard to jurisdictionalclaims in published maps and institutional affiliations.This Springer imprint is published by the registered company Springer Nature Switzerland AGThe registered company address is: Gewerbestrasse 11, 6330 Cham, SwitzerlandPrefaceThis book results from the Second Training School of the European Research andTrainingNetworkUTOPIAE: Uncertainty Treatment andOptimisation in AerospaceEngineering (www.utopiae.eu), at Durham University (United Kingdom) from 2 to6 July 2018. The main focus was on uncertainty quantification.This book consists of 9 chapters providing introductory overviews to topics onuncertainty quantification in engineering, written by topic experts who presentedlectures at the school. The book starts with introductions to Bayesian statisticsand Monte Carlo methods, followed by three chapters on imprecise probabilities:introductions to the general theory and to imprecise Markov chains, and a shortoverview of statistical methods. Then attention shifts to reliability theory and simula-tionmethods for complex systems. The final two chapters are on aspects of aerospaceengineering, considering stochastic model updating from an imprecise Bayesianperspective, and uncertainty quantification for aerospace flight modelling.We are grateful to the chapter authors for their contributions to this book, and fortheir enthusiastic presentations at the training school. We also thank MassimilianoVasile for leading the UTOPIAE network and bringing experts from different fieldstogether to enhance research in the fascinating field of aerospace engineering. Wethank the staff at Springer for their support in the preparation of this book for theSpringer Briefs in Statistics series.This book is sponsored by the UTOPIAE ITN, and the corresponding TrainingSchool was also organised within the UTOPIAE ITN, supported by H2020-MSCA-ITN-2016 UTOPIAE, grant agreement 722734.Durham, United KingdomGhent, BelgiumOctober 2020Louis J. M. AslettFrank P. A. CoolenJasper De BockvContents1 Introduction to Bayesian Statistical Inference . . . . . . . . . . . . . . . . . . . . . 1Georgios P. Karagiannis2 Sampling from Complex Probability Distributions: A MonteCarlo Primer for Engineers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15Louis J. M. Aslett3 Introduction to the Theory of Imprecise Probability . . . . . . . . . . . . . . . 37Erik Quaeghebeur4 Imprecise Discrete-Time Markov Chains . . . . . . . . . . . . . . . . . . . . . . . . . 51Gert de Cooman5 Statistics with Imprecise Probabilities—A Short Survey . . . . . . . . . . . . 67Thomas Augustin6 Reliability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Lisa Jackson and Frank P. A. Coolen7 Simulation Methods for the Analysis of Complex Systems . . . . . . . . . . 95Hindolo George-Williams, T. V. Santhosh, and Edoardo Patelli8 Overview of Stochastic Model Updating in AerospaceApplication Under Uncertainty Treatment . . . . . . . . . . . . . . . . . . . . . . . . 115Sifeng Bi and Michael Beer9 Aerospace Flight Modeling and Experimental Testing . . . . . . . . . . . . . 131Olivier ChazotviiChapter 1Introduction to Bayesian StatisticalInferenceGeorgios P. KaragiannisAbstract We present basic concepts of Bayesian statistical inference. We brieflyintroduce the Bayesian paradigm. We present the conjugate priors; a computationalconvenient way to quantify prior information for tractable Bayesian statistical anal-ysis. We present tools for parametric and predictive inference, and particularly thedesign of point estimators, credible sets, and hypothesis tests. These concepts arepresented in running examples. Supplementary material is available from GitHub.1.1 IntroductionStatistics mainly aim at addressing two major things. First, we wish to learn ordraw conclusions about an unknown quantity, θ ∈ called ‘the parameter’, whichcannot be directly measured or observed, by measuring or observing a sequence ofother quantities called ‘observations (or data, or samples)’ x1:n := (x1, . . . , xn) ∈ Xmwhose generating mechanism is (or can be considered as) stochastically dependenton the quantity of interest θ though a probabilistic model x1:n ∼ f (·|θ). This is aninverse problem since wewish to study the cause θ by knowing its effect x1:n .Wewillrefer to this as parametric inference. Second, we wish to learn the possible values ofa future sequence of observations y1:m ∈ Xm given x1:n . This is a forward problem,and we will call it predictive inference. Here, we present how both inferences can beaddressed in the Bayesian paradigm.1Consider a sequence of observables x1:n := (x1, . . . , xn) generated from a sam-pling distribution f (·|θ) labeled by the unknown parameter θ ∈ . The statisticalmodel m consists of the observations x1:n , and their sampling distribution f (·|θ) ;m = ( f (·|θ); θ ∈ ).1https://github.com/georgios-stats/UTOPIAE-Bayes.G. P. Karagiannis (B)Department of Mathematical Sciences, Durham University, Durham, United Kingdome-mail: georgios.karagiannis@durham.ac.uk© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_112 G. P. KaragiannisUnlike in Frequentist statistics, in Bayesian statistics unknown/uncertain parame-ters are treated as random quantities and hence follow probability distributions. Thisis justified by adopting the subjective interpretation of probability [4], as the degreeof the researcher’s believe about the uncertain parameter θ . Central to the Bayesianparadigm is the specification of the so-called prior distributions dπ(θ) on the uncer-tain parameters θ representing the degree of believe (or state of uncertainty) of theresearcher about the parameter. Different researchers may specify different priorprobabilities, as this is in accordance to the subjective nature of the probability. Thespecification of the prior is discussed in Sect. 1.2.The Bayesian model consists of the statistical model f (x1:n|θ) containing theinformation about θ available from the observed data x1:n , and the prior distributionπ(θ) reflecting the researcher’s believe about θ before the data collection. It is denotedas( f (x1:n|θ), π(θ)) or as{x1:n|θ ∼ f (·|θ)θ ∼ π(·) .Bayesian parametric inference relies on the posterior distribution π(θ |x1:n)whosedensity or mass function (PDF or PMF) is calculated by using the Bayes theoremπ(θ |x1:n) = f (x1:n|θ)π(θ)∫f (x1:n|θ)π(dθ) (1.1)as a tool to invert the conditioning from x1:n|θ to θ |x1:n . Posterior distribution (1.1)quantifies the researcher’s degree of believe after taking into account the observations.By using subjective probability arguments, we can see interpret (1.1) as amechanismthat updates the researcher’s degree of believe from the prior π(θ) to the posteriorπ(θ |x1:n) in the light of the observations collected.Bayesian predictive inference about a future observation y∗ can be addressedbased on the predictive distribution defined asp(y|x1:n) =∫f (y|θ)π(dθ |x1:n) = Eπ ( f (y|θ)|x1:n). (1.2)Essentially, it is the expected value of the sampling distribution averaging outthe uncertain parameter θ with respect to its posterior distribution reflecting theresearcher’s degree of believe.Although the posterior and predictive distributions quantify the researcher’sknowledge, they are not enough to give a solid answer about the quantity to belearned. In what follows we discuss important concepts based on decision theorywhich are used for Bayesian inference.1 Introduction to Bayesian Statistical Inference 31.2 Specification of the PriorPrior distribution π(θ) needs to reflect the researcher’s degree of believe about theuncertain parameter θ ∈ . Sophisticated prior distributions often lead to ineluctableposterior or predictive probabilities, and hence Bayesian analysis. Following, wepresent a computationally convenient class of priors applicable to several scenarios.1.2.1 Conjugate PriorsConjugate priors is a mathematically convenient way to specify the prior model incertain cases. They facilitate the tractable implementation of the Bayesian statisticalanalysis, by leading to computationally tractable posterior distributions.Formally, if F = { f (·|θ); ∀θ ∈ } is a class of parametric models (samplingdistributions), and P = {π(θ |τ); ∀τ } is a class of prior distributions for θ , then theclass P is conjugate for F ifπ(θ |x1:n) ∈ P, ∀ f (·|θ) ∈ F andπ(·) ∈ P.It is straightforward to specify a conjugate prior when the sampling distributionis member of the exponential family. Consider observation xi generated from a sam-pling distribution in the exponential familyxi |θ IID∼ Efk(u, g, h, φ, θ, c); i = 1, . . . , nwith density Efk(x |u, g, h, φ, θ, c) = u(x)g(θ) exp(∑kj=1 c jφ j (θ)(∑ni=1 h j (x)))and g(θ) = 1/ ∫ u(x) exp(∑kj=1 c jφ j (θ)(∑ni=1 h j (x)))dx . The likelihood functionis equal tof (x1:n|θ) =n∏i=1u(xi )g(θ)n exp(k∑j=1c jφ j (θ)(n∑i=1h j (xi ))). (1.3)The conjugate prior, corresponding to likelihood (1.3), admits density of the formπ(θ |τ) = 1K (τ )g(θ)τ0 exp(k∑j=1c jφ j (θ)τ j ) (1.4)where τ = (τ0, . . . , τk) is such that K (τ ) =∫g(θ)τ0 exp(∑kj=1 c jφ j (θ)τ j )dθ <∞. The resulting posterior of θ has the form4 G. P. Karagiannisπ(θ |x1:n, τ ) = 1K (τ ∗)g(θ)τ∗0 exp(k∑j=1c jφ j (θ)τ∗j ))with τ ∗ = (τ ∗0 , τ ∗1 , . . . , τ ∗k ), τ ∗0 = τ0 + n, and τ ∗j =∑ni=1 h j (xi ) + τ j for j = 1,. . . , k.It is easy to see that (1.4) is conjugate to (1.3) as the posterior can be re-writtenas π(θ |x1:n, τ ) = π(θ |τ ∗) where τ ∗ = τ + tn(x1:n), and tn(x1:n) = (n,∑ni=1 h1(xi ),. . . ,∑ni=1 hk(xi )). You can check the demo in.2Example: Bernoulli model (Cont.)Consider observations x1:n = (x1, . . . , xn) ∈ Rn generated from a Bernoulli distri-bution with success rate θ ∈ [0, 1]; i.e., xi |θ ∼ Br(θ), i = 1, . . . , n. Interest lies inspecifying a conjugate prior for θ .The sampling distribution is member of the exponential family, with u(x) = 1,g(θ) = (1 − θ), c1 = 1, φ1(θ) = log( θ1−θ ), h1(x) = x , becausef (x |θ) = Br(x |θ) = θ x (1 − θ)1−x = (1 − θ) exp(log( θ1 − θ )x).The corresponding conjugate prior has PDF such asπ(θ |τ) ∝ (1 − θ)τ0 exp(log( θ1 − θ )τ1) = θ(τ1+1)−1(1 − θ)(τ0−τ1+1)−1,where we recognize Beta distribution π(θ |τ) = Be(θ |a, b), with a = τ1 + 1, b =τ0 − τ1 + 1. Therefore, the posterior distribution isπ(θ |x1:n, τ ) = π(θ |τ0 + n, τ +n∑i=1h(xi )) ∝ θ(τ1+nx̄+1)−1(1 − θ)(τ0+n−τ1−nx̄+1)−1which is Be(θ |a∗, b∗), with a∗ = a + nx̄ , and b∗ = b + n − nx̄ .1.3 Point EstimationOften interest lies in learning the ‘true’ value of the unknown parameter θ ∈ , or thefuture values of a future sequence of observations y1:m ∈ Xm ; this is performed viathe Bayesian point estimator. Here, we demonstrate the theory of the Bayesian pointestimator in parametric inference, and leave the extension to the predictive inferenceto the reader.2 Web-applet: https://georgios-stats-1.shinyapps.io/demo_conjugatepriors/.1 Introduction to Bayesian Statistical Inference 5Bayes (parametric) point estimator of θ ∈ with respect to the loss function(θ, δ) and the posterior distribution π(θ |x1:n) is an Bayes rule δπ which minimizes∫(θ, δ)π(dθ |x1:n); i.e.,δπ (x1:n) = arg min∀δ∈Eπ ((θ, δ)|x1:n) = arg min∀δ∈∫(θ, δ)π(dθ |x1:n). (1.5)Often the accuracy of the Bayes point estimator is represented by its standard error.A commonly accepted metric for the standard error of the j-th dimension of theestimator δπ isseπ (δ j |x1:n) =√MSEπ (δ j |x1:n)where MSEπ (δ j |x1:n) = [Eπ((θ − δ)(θ − δ)|x1:n)] j, j is the mean squared error ofδ j .A number of standard Bayesian point estimates, under different loss functions,are location summary statistics of the posterior distribution (mean, median, mode,quantiles, etc.) You can check the demo in.3The Bayesian estimate of θ with respect to the linear loss (θ, δ) = c1(δ −θ)1θ≤δ(δ) + c2(θ − δ)1{θ≤δ}c(δ) is the c2c1+c2 -th posterior quantile; i.e., π(θ ∈ (−∞,δ(x1:n))|x1:n) = c2c1+c2 . The linear loss function essentially allows the adjustment ofthe penalty between over-estimating and under-estimating θ , by adjusting c1and c2.In particular, for c1 = c2, we get the absolute loss (θ, δ) = |θ − δ| and the posteriorestimator is the posterior medianδ(x1:n) = medianπ (θ |x1:n). (1.6)The absolute loss is more appropriate when over-estimation and under-estimationare of the same concern (as penalized the same).The Bayes estimate δπ (x1:n) of θ with respect to the quadratic loss function(θ, δ) = (θ − δ)2 isδπ (x1:n) = Eπ (θ |x1:n). (1.7)The posterior mean of θ as an estimator of θ essentially minimizes the estimator errorseπ (δ|x1:n), which is equal to the posterior standard error. Obviously, the standarderror of the estimator (1.7) is equal to the posterior standard error. Compared to theabsolute loss, the quadratic loss aims at over-penalizing large but unlikely errors. Infact, quadratic loss aims at minimizing the standard error seπ (δ|x1:n).Finally, the Bayesian estimate of θ with respect to the zero-one loss (θ, δ) =1 − 1Bε (δ)(θ) is the posterior modeδ(x1:n) = modeπ (θ |x1:n) (1.8)as ε → 0.3 Web-applet: https://georgios-stats-1.shinyapps.io/demo_PointEstimation/.6 G. P. KaragiannisExample: Bernoulli model (Cont.)Interest lies in calculating the Bayesian point estimator under the absolute lossfunction. This is the Maximum A posteriori Estimator (the posterior mode). It islog(π(θ |x1:n)) ∝ (nx̄ + a − 1) log(θ) + (n − nx̄ + b − 1) log(1 − θ).For a > 0, b > 0, ddθ log(π(θ |x1:n))|p=δ(x) = 0 implies δ(x) = nx̄+a−1n+a+b−2 . Note that(a.) if a → 1, b → 1 (aka π(θ |a, b) ∝ 1), then δπ (x) = x̄ similar to frequentistsstats; (b.) if a → 0, b → 0 (aka π(θ |a, b) ∝ θ−1(1 − θ)−1), then δ(x) = nx̄−1n−2 ; ifa → 1/2, b → 1/2 (akaπ(θ |a, b) ∝ θ−1/2(1 − θ)−1/2), then δ(x) = nx̄−1/2n−1 ; if n →∞, a > 0, b > 0, then δ(x) = x̄ .1.4 Credible SetsInstead of just reporting parametric (or predictive) point estimates for θ (or y1:m), itis often desirable and more useful to report a subset of values Ca ⊆ (or Ca ⊆ Xm)where the posterior (or predictive) probability that θ ∈ Ca (or y1:m ∈ Ca) is equal toa certain value a reflecting one’s degree of believe.The definition below describes the credible set [1, 5].Definition 1.1 (Posterior Credible Set) A set Ca ⊆ such thatπ(θ ∈ Ca|x1:n) =∫Caπ(dθ |x1:n) ≥ 1 − ais called ‘100(1 − a)%’ posterior credible set for θ , with respect to the posteriordistribution π(dθ |x1:n).In contrast to the frequentist stats, in Bayesian stats we can speak meaningfully ofthe probability that θ is in Ca , because probability 1 − a reflects one’s degree ofbelieve that θ ∈ Ca .Among all the credible sets Ca in Definition 1.1, we are often interested in thosethat have the minimum volume. It can be proved [2] that the highest probabilitydensity (HPD) sets have this property. HPD consider those values of θ correspondingto the highest posterior pdf/pmf (aka the most likely values of θ ).Definition 1.2 (Posterior highest probability density (HPD) set) The 100(1 − a)%highest probability density set for θ ∈ with respect to the posterior distributionπ(θ |x1:n) is the subset Ca of of the formCa = {θ ∈ : π(θ |x1:n) ≥ ka} (1.9)where ka is the largest constant such that1 Introduction to Bayesian Statistical Inference 7π(θ ∈ Ca|x1:n) ≥ 1 − a. (1.10)From the decision theory perspective, HPD set Ca is the Bayes estimate of Ca thecredible interval under the loss function (Ca, θ) = k|Ca| − 1Ca (θ), for k > 0 whichpenalizes sets with larger volumes. The proof is available in [2].Example: Multivariate Normal modelConsider observations x1, . . . , xn independently drawn fromaq-dimensional normalNq(μ,) with unknown μ ∈ Rq , q ≥ 1, and known , μ0, 0. Assume prior μ ∼Nq(μ0, 0). Interest lies in calculating the Ca parametric HPD credible interval forμ.The posterior PDF of μ isπ(μ|x1:n) ∝ f (x1:n|μ)π(μ) =n∏i=1Nq(xi |μ,)Nq(μ|μ0, 0)∝ exp(−12(μ − μ̂n)T ̂−1n (μ − μ̂n)) ∝ Nq(μ|μ̂n, ̂n)where n = (n−1 + −10 )−1, and μ̂n = ̂n(n−1 x̄ + −10 μ0). So μ|x1:n ∼ Nq(μ̂n, ̂n).From Definition 1.2, the credible set has the formCa = {μ ∈ Rq : π(μ|x1:n) ≥ ka}= {μ ∈ Rq : (μ − μ̂n)T ̂−1n (μ − μ̂n) ≤ − log(2π det(̂n)))ka = k̃a}where ka is the greatest value satisfyingπNq (μ̂n ,̂n)(μ ∈ Ca|x1:n) ≥ 1 − a ⇐⇒πχ2q ((μ − μ̂n)T ̂−1n (μ − μ̂n) ≤ k̃a) ≥ 1 − a. (1.11)Here, (μ − μ̂n)T ̂−1n (μ − μ̂n) ∼ χ2q as a sum of squares of independent standardnormal random variables, and hence k̃a is the 1 − a-th quantile of the χ2q distribution;i.e., k̃a = χ2q,1−a . Therefore, Ca parametric HPD credible set for μ isCa = {μ ∈ Rq : (μ − μ̂n)T ̂−1n (μ − μ̂n) ≤ χ2q,1−a}8 G. P. KaragiannisFig. 1.1 Schematic of 1D HPD setIn real applications, the calculation of the credible interval might be intractable,due to the inversion in (1.9) or integration in (1.10). Below, we present a Naivealgorithm [1] that can be implemented in a computer.4• Create a routine which computes all solutions θ∗ to the equation π(θ |x1:n) = ka ,for a given ka . Typically, Ca = {θ ∈ : π(θ |x1:n) ≥ ka} can be constructed fromthose solutions.• Create a routine which computes π(θ ∈ Ca|x1:n) =∫θ∈Ca π(θ |x1:n)dθ• Numerically solve the equation π(θ ∈ Ca|x1:n) = 1 − a as ka varies.Figure 1.1 demonstrates the above procedure in 1D unimodal and tri-modal cases.Specifically, the red horizontal bar denotes ka moves upwards, and intersects thedensity at locations which are the potential boundaries of Ca . The bar stops to movewhen the total density above regions of the parametric space is equal to 1 − α. TheHPD credible set results as the union of these sub-regions. You can check the demoin.5Theorem 1.1 suggests a computationally convenient way to calculate HPD cred-ible intervals in 1D, and unimodal cases. The proof is available in [3].Theorem 1.1 Let θ follows a distribution with unimodal density π(θ |x1:n). If theinterval Ca = [L ,U ] satisfies1.∫ UL π(θ |x1:n)dθ = 1 − a,2. π(U ) = π(L) > 0, and3. θmode ∈ (L ,U ), where θmode is the mode of π(θ |x1:n),then interval Ca = [L ,U ] is the HPD interval of θ with respect to π(θ |x1:n).4 Web-applet: https://georgios-stats-1.shinyapps.io/demo_crediblesets/.5 Web-applet: https://georgios-stats-1.shinyapps.io/demo_crediblesets/.1 Introduction to Bayesian Statistical Inference 9Example: Bernoulli model (Cont.)Interest lies in calculating the 2-sides 95% HPD interval for θ , given a sample withn = 30, and ∑30i=1 xi = 15, and prior hyper-parameters a = b = 2.The posterior distribution of θ is Be(a + nx̄ = 17, b + n − nx̄ = 17), which is1D and unimodal; hence we use Theorem 1.1. It is1 − a =∫ ULBe(θ |17, 17)dθ = Be(θ < U |17, 17) − Be(p < L|17, 17).Note that Beta PDF is symmetric around 0.5 when a∗ = b∗, and so is here whereBe(17, 17). Then,1 − a = Be(θ < U |17, 17) − (1 − Be(θ < U |17, 17)) = 2Be(θ < U |17, 17) − 1so Be(θ < U |17, 17) = 1 − a/2 and L = 1 −U . For a = 0.95, the 95% posteriorcredible interval for θ is [L ,U ] = [0.36, 0.64].Remark 1.1 Predictive credible sets for a future sequence of observations y1:m , aredefined and constructed as parametric ones by replacing θ with y1:m and π(x1:n|θ)with p(y1:m |x1:n) in Definitions 1.1 and 1.2, and their consequences in this section.It is left as an Exercise.1.5 Hypothesis TestOften there is interest in reducing the overall parametric space (aka the set ofpossible values of that the uncertain parameter θ can take) to a smaller subset. Forinstance; whether the proportion of Brexiters is larger than 0.5 (p > 0.5) or not(p ≤ 0.5).Such a decision can be formulated as a hypothesis test [1], namely the decisionprocedure of choosing between two non-overlapping hypothesesH0 : θ ∈ 0 vs H1 : θ ∈ 1 (1.12)where {0,1} partitions the space . Typically, hypotheses, {Hk}, are categorizedin three categories. Single hypothesis for θ is called the hypothesis where j = {θ j }contains a single element. Composite hypothesis for θ is called the hypothesis where j ⊆ contains many elements. General alternative hypothesis for θ is called thecomposite hypothesis where 1 = − {θ0} when it is compared against a singlehypothesis H0 : θ = θ0. It is denoted as H1 : θ = θ0.Based on the partitioning implied by (1.12), the overall prior π can be expressedas π(θ) = π0 × π0(θ) + π1 × π1(θ)where πk =∫kπ(dθ), and πk(θ) = π(θ)1k (θ)∫kπ(dθ) .10 G. P. KaragiannisHere, π0, and π1 describe the prior probabilities on H0 and H1, respectively, whileπ0(θ) and π1(θ) describe how the prior mass is spread out over the hypotheses H0and H1, respectively.We could see the hypothesis testing (1.12) as parametric point inference about theindicator function11(θ) ={0 , θ ∈ 01 , θ ∈ 1 . (1.13)To estimate (1.13), a reasonable loss function (θ, δ) would be the cI − cII lossfunction(θ, δ) =⎧⎪⎪⎪⎨⎪⎪⎪⎩0 , if θ ∈ 0, δ = 00 , if θ /∈ 0, δ = 1cII , if θ /∈ 0, δ = 0cI , if θ ∈ 0, δ = 1(1.14)where cI > 0 and cII > 0 are specified by the researcher. Here, cI > 0 (and cII > 0)denote the loss if we decide to accept H0 (and H1) while the correct answer would beto choose H1 (H0). According to (1.5), under (1.14), the Bayes estimator of (1.13) isδ(x1:n) ={0 , if π(θ ∈ 0|x1:n) > cIIcII+cI1 , otherwise(1.15)where π(θ ∈ 0|x1:n) =∫0π(dθ |x1:n). In other words, hypothesis H1 is acceptedif π(θ∈0|x1:n)π(θ∈1|x1:n) cIIcIπ1π0, and accept H1 if otherwise. Alternatively, Jeffreys [6] developed a scale rule(Table 1.1) to judge the strength of evidence in favor of H0 or against H0 brought bythe data. Although Jeffreys’ rule avoids the need to specify cI and cI I , it is a heuris-tic rule-of-thumb guide, not based on decision theory concepts, and hence manyresearchers argue against its use.1 Introduction to Bayesian Statistical Inference 11Table 1.1 Jeffreys’ scale rule [6]B01 log10(B01) Strength of evidence(1,+∞) (0,+∞) H0 is supported(10−1/2, 1) (−1/2, 0) Evidence against H0: notworth more than a bare(10−1, 10−1/2) (−1,−1/2) Evidence against H0:substantial(10−3/2, 10−1) (−3/2,−1) Evidence against H0: strong(10−2, 10−3/2) (−2,−3/2) Evidence against H0: verystrong(0, 10−2) (−∞,−2) Evidence against H0: decisiveExample: Bernoulli model (Cont.)We are interested in testing the hypotheses H0 : θ = 0.5 and H1 : θ = 0.5, giventhat π0 = 1/2, and using the cI − cII loss function with cI = cII. Here, 0 ={0.5} and 1 = [0, 0.5) ∪ (0.5, 1]. The overall prior is π(θ) = π01θ0(θ) + (1 −π0)Be(θ |a, b). The Bayes factor isB01(x1:n) =∏ni=1 Br(xi |θ0)∫(0,1)∏ni=1 Br(xi |θ)Be(θ |a, b)dθ= θx∗0 (1 − θ0)n−x∗B(nx̄ + a, n − nx̄ + b)/B(a, b) .Given a = b = 2, n = 30, and ∑30i=1 xi = 15, it is B01(x1:n) = 18.47 > cII/cI = 1.Hence, we accept H1.1.5.1 Model SelectionOften the researcher is uncertain which statistical model (sampling distribution) canbetter represent the real data generating process. There is a setM = {m1,m2, . . .} ofcandidate statistical models mk = { fk(·|ϕk); ϕk ∈ Φk}, where fk(·|ϕk) denotes thesampling distribution, and ϕk denotes the unknown parameters for k = 1, 2, . . . Letπk = π(mk) denote the marginal model prior and πk(ϕk) = π(ϕk |mk) denote theprior of the unknown parameters ϕk of given model mk .Selection of the ‘best’ model from a set of available candidate models can beaddressed via hypothesis testing. For simplicity, we consider there are only twomodels m0 and m1 with unknown parameters ϑ0 ∈ Φ0 and ϑ1 ∈ Φ1. Then, modelselection is performed as a hypothesis testH0 : (m, ϕ) ∈ 0 vs H1 : (m, ϕ) ∈ 1 (1.18)12 G. P. Karagianniswherek = {mk} × Φk , = ∪kk . The overall joint prior is specified asπ(m, ϕ) =π0 × π0(ϕ0) + π1 × π1(ϕ1) on (m, ϕ) ∈ where = ∪kk , where πk(ϕk) =π(m,ϕ)1mk (m)∫Φkπ(m,dϕ) on ϕk ∈ Φk , and πk =∫kπ(mk, dϕk). Now the model selection prob-lem has been translated into a hypothesis test.Example: Negative binomial vs. Poisson model [2]We are interested in testing the hypothesesH0 : xi |φ ∼ Nb(φ, 1), φ > 0, vs. H1 : xi |λ ∼ Pn(λ), λ > 0by using the cI − cII loss functionwith cI = cII. Consider two observations x1 = x2 =2 are available. Consider overall prior π(θ) with density π(θ) = π0Be(φ|a0, b0) +π1Ga(λ|a1, b1) with π0 = π1 = 0.5.This is a composite vs. composite hypothesis test. It is∫0f (x1:n|ϕ0)π0(dϕ0) = (a0 + b0)(a0)(b0)∫ 10φn+a0−1(1 − φ)nx̄+b0−1dφ= (a0 + b0)(a0)(b0)(n + a0)(nx̄ + b0)(n + nx̄ + a0 + b0)∫1f (x1:n|ϕ1)π1(dϕ1) = ba11(a1)(n + b1)nx̄+a1∫ ∞0λnx̄+a1−1 exp(−(n + b1)λ)dλ= (nx̄ + a1)(a1)(n + b1)nx̄+a11∏ni=1 xi !and hence B01(x1:n) = (a0+b0)(a0)(b0) (n+a0)(nx̄+b0)(n+nx̄+a0+b0) (a1)(n+b1)nx̄+a1(nx̄+a1)∏ni=1 xi ! . It is B01(x1:n)= 0.29 > 1, and hence I accept H1 and the Poisson model.References1. James O Berger. Statistical decision theory and Bayesian analysis. Springer Science &BusinessMedia, 2013.2. José M Bernardo and Adrian FM Smith. Bayesian theory, volume 405. John Wiley & Sons,2009.3. George Casella and Roger L Berger. Statistical inference, volume 2. Duxbury Pacific Grove,CA, 2002.4. Bruno De Finetti. Theory of probability: a critical introductory treatment, volume 6. JohnWiley& Sons, 2017.1 Introduction to Bayesian Statistical Inference 135. Morris H DeGroot. Optimal statistical decisions, volume 82. John Wiley & Sons, 2005.6. Harold Jeffreys. The theory of probability. OUP Oxford, 1998.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 2Sampling from Complex ProbabilityDistributions: A Monte Carlo Primer forEngineersLouis J. M. AslettAbstract Models which are constructed to represent the uncertainty arising in engi-neered systems can often be quite complex to ensure they provide a reasonablyfaithful reflection of the real-world system. As a result, even computation of simpleexpectations, event probabilities, variances, or integration over utilities for a decisionproblem can be analytically intractable. Indeed, such models are often sufficientlyhigh dimensional that even traditional numerical methods perform poorly. However,access to random samples drawn from the probability model under study typicallysimplifies such problems substantially. The methodologies to generate and use suchsamples fall under the stable of techniques usually referred to as ‘Monte Carlo meth-ods’. This chapter provides amotivation, simple primer introduction to the basics, andsign-posts to further reading and literature onMonte Carlo methods, in a manner thatshould be accessible to those with an engineering mathematics background. Thereis deliberately informal mathematical presentation which avoids measure-theoreticformalism. The accompanying lecture can be viewed at https://www.louisaslett.com/Courses/UTOPIAE/.2.1 MotivationThere is a natural tension when constructing a probabilistic model with the aim ofencapsulating the uncertainty in an engineered system: on the one hand, there is adesire to capture every nuance of the system to fully reflect all knowledge aboutits behaviour; on the other, there is a drive towards parsimony for reasons of inter-pretability, robustness, and computability. Interpretability and robustness are impor-tant goals and should indeed guide a reduction in model complexity, but reducingmodel complexity purely to enable computability would seem a hinderance, espe-cially if that parsimony impedes answering the research questions at hand since, putsimply, ‘reality can be complicated’ [7]. As such, the methodology of this chapterL. J. M. Aslett (B)Department of Mathematical Sciences, Durham University, Durham, United Kingdome-mail: louis.aslett@durham.ac.uk© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_21516 L. J. M. Aslettshould not be employed simply to enable an inappropriately complex model, butrather serves to facilitate the use of models which are complex enough when judgedby purely subject matter and statistical concerns.Monte Carlo methods have played a crucial role in a vast array of applicationsof statistical methodology, from the prediction of future marine species discoveries[29] through to reconstruction of the ancient climate on Earth [16]; from criminaljustice offending risk [17] to inferring networks of corporate governance throughthe financial crash [12]; and from estimating bounds on engineering system survivalfunctions [11] to the assessment of offshore oil production availability [30]. Theutility of Monte Carlo in these applications varies substantially, from estimation ofconfidence intervals and event probabilities, through optimisation methods to fullevaluation of Bayesian posterior distributions for parameter inference.With this breadth of application in mind, we may assume hereinafter that wehave a probabilistic model for some engineered system of interest which—afterconsidering all subject matter and statistical concerns—is too complex to be ableto compute relevant quantities of interest (be they event probabilities, confidenceintervals, posterior distributions, etc.). As a concrete example, if onewere to constructa Bayesian model of reliability using ideas introduced in Chap. 1, then our modelwould comprise some prior distribution over the vector of model parameters, π(θ),together with a generative model for the failure time depending on those parameters,π(t | θ). After collecting some lifetime data t = {t1, . . . , tn}, themost simple researchquestion of interest may be the posterior expected value of the parameters:Eπ [θ ] =∫θ π(θ | t) dθ = 1c∫θ π(θ)n∏i=1π(ti | θ) dθ (2.1)where is the space of all possible parameter values and c is a normalising constant.Indeed, it is traditional in Monte Carlo literature to focus attention on the compu-tation of expectations with respect to some probability density under consideration,which need not necessarily be a Bayesian posterior. That is, given a general proba-bility model π(x), x ∈ , and a functional f : → R, interest is typically in:Eπ [ f (X)] :=∫f (x)π(x) dx (2.2)and this is the perspective that will be adopted in this chapter.We complete our motivation of Monte Carlo in this Section by highlighting thegenerality of expectations of the form (2.2), followedby a short discussion of standardnumerical integration techniques. In Sect. 2.2, theMonte Carlo estimator and its erroranalysis are recapped and contrastedwith numerical integration. The coremethods ofMonte Carlo simulation are introduced in Sect. 2.3, with pointers to more advancedmaterial in Sect. 2.4. Note that we will in places abuse formal notation where webelieve it aids intuitive understanding since the goal of this chapter is to be a basic2 Sampling from Complex Probability Distributions: A Monte … 17primer, not a rigorous treatment.1 A first course in probability and statistics areassumed background.The accompanying lecture from the UTOPIAE training school can be viewed athttps://www.louisaslett.com/Courses/UTOPIAE/.2.1.1 Generality of ExpectationsThe formulation in (2.2) may appear rather restrictive to the uninitiated reader. How-ever, considering only expectations of this form does not result in any loss of gener-ality. For example, (re-)defining:π(x) := 1cπ(x)n∏i=1π(ti | x)f (x) := xmeans that (2.2) simply becomes the posterior expectation in (2.1). However, oneshould note that arbitrary statements of probability are also computable as expecta-tions. That is,P(X < a) =∫ a−∞π(x) dx =∫I(−∞,a](x)π(x) dx = Eπ [I(−∞,a](X)]where for a general set E ⊆ ,IE (x) :={1 if x ∈ E0 if x /∈ EThat is, to evaluate the probability of an arbitrary event, P(X ∈ E), simply setf (X) := IE (X) when evaluating (2.2).2.1.2 Why Consider Monte Carlo?In some special cases, the integral (2.2) may have an analytical solution and in suchsituations one should not resort to Monte Carlo or other methods. When there isno known analytical form for the integral, a reader with a general mathematical1 For example, we will write ‘P(X = x)’ even where X is continuous to emphasise the link to thedensity function and will use π(x) to reference both a target distribution or prior where the meaningis clear from context. For the more advanced reader there are already many excellent more rigoroustreatments in the literature, some of which we reference towards the end.18 L. J. M. Aslettbackground may be tempted to reach for a numerical integration method, such as asimple mid-point Riemann integral or a more sophisticated quadrature approach.Consider the mid-point Riemann integral in the simple 1-dimensional setting.Letting g(x) := f (x)π(x), then the expectation would be approximated using nevaluations by: ∫ bag(x) dx ≈ b − ann∑j=1g(x j ), (2.3)wherex j := a + b − an(j − 12).The absolute error in using (2.3) is bounded [24, Theorem 7.1]:∣∣∣∣∣∣∫ bag(x) dx − b − ann∑j=1g(x j )∣∣∣∣∣∣ ≤(b − a)324n2maxa≤z≤b |g′′(z)|.Clearly, (b−a)324 maxa≤z≤b |g′′(z)| is fixed by the problem at hand and cannot be alteredby the engineer, sowe achieve the accuracywe require by controlling n−2—that is, byusing a finer grid to compute the integral. As such, we say the error in the mid-pointRiemann integral in 1 dimension is O (n−2)—that is, if double the computationaleffort is expended by computing on a grid of twice as many points (2n), then theworst case error is reduced by a factor of 4. This fast reduction in error and an explicitbound on it are very attractive properties.However, as the dimension of x increases, the Riemann integral’s effectivenessdiminishes substantially. In general, the error of mid-point Riemann integration ind-dimensions is O (n−2/d). For example, even in a modest 10-dimensional problem,when the computational effort is doubled the worst case error is only reduced by afactor of ≈ 1.15. Put another way, to halve the worst case error in a 10-dimensionalproblem requires exp(102 log 2) = 32 times the computational effort. This problemhas been coined the ‘curse of dimensionality’.Of course, the Riemann integral is not the best numerical integration method,but even Simpson’s rule only improves this to O (n−4/d). In general Bakhvalov’sTheorem bounds all possible quadraturemethods byO(n−r/d), where r is the numberof continuous derivatives of g(·)which exist and are exploited by the quadrature rule[24].The striking result which motivates the study of Monte Carlo methods is that for ad-dimensional problem, the (mean-square)2 error is O (n−1/2). The most importantpoint to note is the absence of d in the order of the error: increasing the computa-tional effort by some fixed amount has the same relative effect on the worst caseerror regardless of dimension. Of course, the devil in the detail is that the constant2 Note that randomised simulation methods such as Monte Carlo typically report mean-square errorrather than absolute error bounds.2 Sampling from Complex Probability Distributions: A Monte … 19Fig. 2.1 The order of error reduction—that is, only the leading O (n f (d)) term—is plotted againstdifferent computational effort n. Note that all these curves would be multiplied by a different (fixed)problem dependent constantfactor which we are ignoring in that statement almost certainly has some dimensiondependence, but this is true for quadrature methods too. Figure2.1 illustrates thedifferences.Consequently, Monte Carlo methods are well suited to address the problem ofanalysing complex probabilistic models of engineered systems, since this is preciselya setting where the parameter dimension is likely to be large.2.2 Monte Carlo EstimatorsThe standard Monte Carlo estimator of the integral (2.2) isμ ∫f (x)π(x) dx ≈ 1nn∑j=1f (x j ) μ̂, (2.4)where x j ∼ π(·). In other words, the problem of integration is transformed insteadinto the problem of drawing random samples x j distributed according to the proba-bility density π(·). Importantly, this estimator is unbiased, that is, E[μ̂] = μ.If the samples x j are independently and identically distributed (iid) according toπ(·), then the root mean-square error of the estimator μ̂ isRMSE :=√√√√√Eπ⎡⎣⎛⎝∫f (x)π(x) dx − 1nn∑j=1f (x j )⎞⎠2⎤⎦ = σ√n,20 L. J. M. Aslettwhere σ 2 = Varπ ( f (X)). Again, part of this error is (mostly) inherent to the prob-lem3—σ in this case—so that we achieve desired accuracy by controlling n−1/2.There are at least three very attractive conclusions we can draw from this form:1. as mentioned already, the relative error reduction achieved by additional com-putational effort is independent of dimension;2. there is no explicit dependence on how smooth the functional, f (·), or probabilitydensity, π(·), are (though these may influence σ );3. in contrast to quadrature methods, an estimate of the error can be computedfrom the work already done to compute the integral, by computing the empiricalstandard deviation of the functional of the samples drawn from π(·).Although an absolute error is not available for a randomised method like this, asimple application of Chebyshev’s inequality does provide a probabilistic bound onthe absolute error exceeding a desired tolerance:P(|μ̂ − μ| ≥ ε) ≤ Eπ [(μ̂ − μ)2]ε2= σ2nε2.Indeed, it is also possible to invoke the iid Central Limit Theorem so that asymp-totically,P(μ̂ − μσn−1/2≤ z)n→∞−−−→ (z),where (z) denotes the standard Normal cumulative distribution function (CDF).This enables the formation of confidence intervals for μ based on large n samples.The discussion to date has tacitly assumed that simulating from arbitrary proba-bility distributions π(·) is possible and relatively efficient. In fact, most Monte Carloresearch is devoted to this effort since, as touched on above, there is rich and well-established theory when such samples are available. Therefore, for the remainder ofthis chapter, our attention turns away from discussion of the integrals which are ofprimary interest and focuses on the problem of simulating from arbitrary probabilitydistributions π(·). Once these samples are available, the results above can be used toanalyse the resulting estimators.2.3 Simple Monte Carlo Sampling MethodsIn this section we introduce some simple Monte Carlo methods which enable sam-pling from a wide array of probability distributions. Note that understanding thesesimple methods is crucial as they are extensively used as building blocks of moresophisticated sampling methodology.3 There are advanced Monte Carlo methods which can reduce this variance, but this is beyond thescope of this chapter. See for example [24, Chap. 8].2 Sampling from Complex Probability Distributions: A Monte … 21Almost all Monte Carlo procedures start from the assumption that we have avail-able an unlimited stream of iid uniformly distributed values, typically on the interval[0, 1] ⊂ R. How to generate such an iid stream is beyond the scope of this intro-ductory chapter, but the interested reader may consult [13, Chaps. 1–3] and [21].Arguably the current gold standard algorithm remains that in [22]. Typically, theaverage user of Monte Carlo need not worry about such issues and may rely on thehigh quality generators built into software such as R [26].Thus the objective hereinafter is to study how to convert a stream ui ∼ Unif(0, 1)into a stream x j ∼ π(·), where x j is generated by some algorithm depending on thestream of ui . In more advanced methods (see MCMC), x j may also depend on x j−1or even x1, . . . , x j−1.2.3.1 Inverse SamplingArguably the simplest example of generating non-uniform random variates is inversesampling, which typically applies only to 1-dimensional probability distributions(thoughhigher dimensional extensions have been studied). Let F(x) := P(X ≤ x)bethe cumulative distribution function (CDF) for the target probability density functionπ(·). Then, inverse sampling requires the inverse of the cdf, F−1(·), which is thenapplied to a uniform random draw. Precisely, see Algorithm2.1.Algorithm 2.1 Inverse sampling algorithm1: procedure Inverse sampling(F−1(·)) Generate random sample from distribution withinverse CDF F−1(·)2: u ∼ Unif(0, 1)3: x ← F−1(u)4: return x5: end procedureTo prove that the sample returned byAlgorithm2.1 is distributed according toπ(·)is straight-forward. We do so by computing the CDF, P(X ≤ x), of the X generatedby this algorithm and show that this agrees with the CDF of π(·). The first stepsubstitutes X = F−1(U ), where U ∼ Unif(0, 1), as per the algorithm:P( X ≤ x ) = P( F−1(U ) ≤ x )= P( F(F−1(U )) ≤ F(x) ) applying F(·) to both sides= P(U ≤ F(x) ) Uniform CDF P(U ≤ u) = u= F(x).Note that applying F(·) to both sides in the second line is valid, since the cumulativedistribution function is a non-decreasing function by definition.22 L. J. M. AslettFig. 2.2 Inverse sampling for both a continuous distribution function (left) and one containingjump discontinuities and regions of zero probability (right). Uniform random draws u are sampledand inverted through the distribution function in the obvious way (left), or by taking the infimumover values of x such that F(x) ≥ u (right). In the right illustration, the hypothetical (‘Hypothetical’since strictly speaking this is an event of probability zero) u1 coincides with the value at whichF(x) is constant and u2 lies within the jump discontinuityOne subtlety to be aware of is that for discrete distributions or continuous distri-butions with jump discontinuities or areas of no support, we must define:F−1(u) = inf{x : F(x) ≥ u}, ∀ u ∈ [0, 1].It may be tempting when F−1(·) is not available to use a numerical solver to solveF(x) = u in place of line 3 in Algorithm2.1. However, caution is required sincethis can result in bias [10, p. 31]. The procedure of inverse sampling is illustrated inFig. 2.2.Notice that this is univariate, yet earlierwe saw that numerical integrationwill givebetter error bounds than Monte Carlo for low dimensional problems—as such onemay choose not to use inverse sampling to actually evaluate univariate expectations.However, we often need a set of random draws from non-uniform univariate distri-butions which feed into a broader Monte Carlo algorithm, which is itself sampling inhigher dimensions: in such situations inverse sampling is very useful. Indeed, if youuse the rnorm function in R [26], it has used inverse sampling to generate randomdraws from the Normal distribution since 2003 (see /src/nmath/snorm.c lines265–270), prior to that using [18] since at least v0.62 in 1998.A final comment: inverse sampling is a special case of general transformationsampling. If one can generate samples from one distribution, there may be an appro-priate transformation to turn these into samples from another distribution that maybe more tractable or faster than inverse sampling. For further details, see for example[24, Chap. 4.6].2 Sampling from Complex Probability Distributions: A Monte … 232.3.1.1 ExampleIn order to use inversion sampling for a Weibull distribution with shape k and scaleσ , X ∼ Weibull(k, σ ), we note thatπ(x) = kσ( xσ)k−1e−(xσ )k, x ∈ [0,∞), σ > 0, k > 0F(x) = 1 − exp{−( xσ)k}.To find F−1(u), set 1 − exp{− ( xσ)k} = u and solve for x :=⇒ x = F−1(u) = σ (− log(1 − u))1/k ∼ π(·). (2.5)In order to generate samples from the Weibull we, therefore, take values, u, from aUniform random number stream and transform them using (2.5).2.3.2 Rejection SamplingOur first higher dimensional method is an elegant algorithm, which actually crops upin more advanced guises at the cutting edge of modern Monte Carlo methods (e.g.[9, 25]). Here, the goal is to find another distribution, say π̃(·), which is easier tosample from (perhaps even using inverse sampling) and where we can construct abound on the density function:π(x) ≤ cπ̃(x) ∀ x ∈ , (2.6)where c < ∞ and where π and π̃ need not be normalised probability densities. Wecall π̃(·) the ‘proposal’ density, since samples will be drawn from this and thenexactly the correct proportion of them retained in order to end up with a stream ofsamples from π(·). The full procedure is detailed in Algorithm2.2.We will proceed based on the assumption that π(·) and π̃(·) are normalised den-sities. However, note that the algorithm is also valid for un-normalised densities, solong as there still exists a c satisfying (2.6) for the un-normalised densities.The efficiency of the algorithm hinges entirely on the value of c, so that it shouldbe chosen as small as possible. This is because, letting A be the random variable foracceptance of a proposed sample X , the acceptance probability is (abusing notationto aid intuition):24 L. J. M. AslettAlgorithm 2.2 Rejection sampling algorithm1: procedure Rejection sampling(π(·), π̃(·), c) Generate random sample from distri-bution with unnormalised density π(·)2: a ← FALSE3: while a = FALSE do Repeat until acceptance4: u ∼ Unif(0, 1)5: x ∼ π̃(·) Propose a possible sample6: if u ≤ π(x)cπ̃(x) then Accept or reject proposal?7: a ← TRUE8: end if9: end while10: return x11: end procedureP(A = 1) =∫P(A = 1 | X = x)︸︷︷︸Prob line 6 of Alg 2 gives TRUE.P(X = x)︸︷︷︸Proposal density π̃ .dx=∫P(U ≤ π(x)cπ̃ (x))︸︷︷︸Uniform CDF, P(U≤u)=F(u)=u.π̃ (x) dx=∫π(x)cπ̃(x)π̃(x) dx= 1c∫π(x) dx= 1c, (2.7)where is the support of π̃(·), π̃(x) > 0, ∀ x ∈ . The final line follows becausethe integral of a density over the whole space is 1.Hence, the number of iterations of the loop on lines 3–9 in Algorithm2.2 whichmust be performed to return a single sample from π(·) is Geometrically distributedwith parameter 1c . Therefore, the expected number of random number generationsand function evaluations which must be performed is 2c per sample from π(·).To see that Algorithm2.2 does indeed give a sample from π(·), we note that thesamples returned are only those which are accepted, so we condition on this event:P(X ∈ E | A = 1) = P(A = 1 | X ∈ E) P(X ∈ E)P(A = 1) ∀ E ∈ B=∫Eπ(x)cπ̃(x) π̃ (x) dx1c=∫Eπ(x) dx .2 Sampling from Complex Probability Distributions: A Monte … 25Fig. 2.3 Geometric interpretation of rejection sampling. First y is sampled from π̃ and then auniform is sampled along the vertical dashed line at that location (i.e. between 0 and cπ̃(y)). If theuniform sample falls below π then we accept and otherwise we reject. It is therefore clear that thecloser cπ̃ ‘hugs’ π the more efficient the rejection samplerThe last line is the probability of event E under the distribution with density π(·), asrequired.There is a nice geometric interpretation of the rejection sampling algorithmwhichaids greatly with intuition. Notice that the condition in line 6 can be rewrittenucπ̃(x) ≤ π(x). This means ucπ̃ (x) is a uniform random number in the interval[0, cπ̃ (x)], so that we can view rejection sampling as first drawing a value fromπ̃(x), then moving it up to a uniformly distributed height under the curve cπ̃ (x).The consequence of this is that we are effectively sampling uniformly points underthe curve cπ̃ (x) and accepting those that fall under the curve π(x), as depicted inFig. 2.3.Care is required with rejection sampling in high dimensions because it is quiteeasy for the acceptance probability to become so small as to make the techniqueimpractical. We will see that, as well as how to implement rejection sampling, in thefollowing example.2.3.2.1 ExampleConsider the problem of sampling from a zero mean d-dimensional multivariateNormal distribution, having density:π(x) = (2π)−d/2 det()−1/2 exp(−12xT−1x), x ∈ Rd ,where is a d × d symmetric positive semi-definite covariance matrix. It is com-paratively easy to sample univariate Normal random variables (e.g. using inversesampling in R [26] as mentioned earlier, or via a transformation type approach like[3]). Thus we could consider using a multivariate Normal with diagonal covariance,26 L. J. M. Aslettσ 2 I , as a proposal, because this simply requires sampling d univariate Normal ran-dom variables.This would mean we need to determine c < ∞ such thatc det(σ 2 I )−1/2 exp(−12xTσ−2 Ix)≥ det()−1/2 exp(−12xT−1x)∀ x ∈ Rd(2.8) is symmetric, so it has eigendecomposition = Q QT =⇒ −1 = Q −1QT,where Q is an orthogonal matrix and is a diagonal matrix with entries consistingof the eigenvalues λ1, . . . , λd . The orthogonal transformation y = QTx also spansRd , so that (2.8) ⇐⇒cσ−d exp(−12yTσ−2 Iy)≥(d∏i=1λi)−1/2exp(−12yT −1y), ∀ y ∈ Rd⇐⇒ 2 log c ≥d∑i=1(σ−2 − λ−1i )y2i + 2d log σ −d∑i=1log λi .If σ−2 < λ−1i for any i , then the right-hand side cannot be bounded above (since theinequality must hold ∀ yi ∈ R), so we must have maxi λi < σ 2 and then clearly cis minimised for σ 2 = maxi λi . Since every term in the first sum of the right-handside is necessarily negative, the right-hand side is maximal for yi = 0 ∀ i , so that theoptimal c isc =(maxiλi)d/2 ( d∏i=1λi)−1/2,when σ 2 = maxi λi .In summary, there is a constraint on our proposal π̃(·) when it is an uncorrelatedmultivariateNormal density, or else it cannot boundπ(·).Moreover, we can explicitlycompute the optimal proposal variance, σ 2, to give us the highest possible acceptancerate.To make this example concrete, consider rejection sampling in this setting where =⎛⎜⎜⎜⎝1 0.9 · · · 0.90.9 1 · · · 0.9....... . ....0.9 0.9 · · · 1.⎞⎟⎟⎟⎠Note that can be written as 0.1I + B, where B is a matrix with 0.9 in everyelement. The rank of B is 1, so it has a single non-zero eigenvalue which musttherefore equal tr(B) = 0.9d, and the eigenvalues of 0.1I are all 0.1. Further-2 Sampling from Complex Probability Distributions: A Monte … 27Table 2.1 Optimal proposal variance and acceptance probability for rejection sampling a correlatedmultivariate Normal distribution using an uncorrelated multivariate Normal proposald σ 2 = maxi λi Acceptance probability 1c1 1 12 1.9 0.2293 2.8 0.0364 3.7 0.0045 4.6 4.45 × 10−5.........10 9.1 1.53 × 10−9more, 0.1I and B commute, therefore the eigenvalues of are the sum of theseeigenvalues: that is, λ1 = 0.9d + 0.1 and λi = 0.1 ∀ i = 1. As the dimension ofπ(·) increases, the spectral gap increases linearly and thus c grows very fast:c = (10λ1)(d−1)/2 = (9d + 1)(d−1)/2. Indeed, this is faster than exponential and fasterthan factorial growth! Consequently, for growing dimension, the acceptance proba-bility falls super-exponentially fast—not a desirable property. See Table2.1 for someexample values.A whimsical observation to emphasise the problem: a modern laptop can pro-duce roughly 15 million univariate Normal samples per second and the universeis estimated to be 4.32 × 1017 seconds old. Ignoring the time to evaluate the uni-form draw u or acceptance/rejection, this means the expected number of samplesthat would be generated by Algorithm2.2 for this multivariate Normal problem ind-dimensions—if run for as long as the universe has existed—would be1.5 × 107 × 4.32 × 1017d(9d + 1)(d−1)/2 .Consequently, even knowing the exactly optimal choice for σ 2 in our proposal, thiswould only be expected to render 5 samples for a 21-dimensionalmultivariateNormalwith the innocuous looking given above—rejection sampling in high dimensionscan be problematic!2.3.3 Importance SamplingThe final core standard Monte Carlo method we cover in this primer also starts fromthe perspective of having a proposal density π̃(·), thoughwe no longer require it to beable to bound π(·). Importance sampling then dispenses with the notion of directlygenerating iid samples from π(·) and focuses on their use: in computing expectationsusing those samples in (2.4).Consequently, importance samplingweights the samples28 L. J. M. Aslettfrom π̃(·) in precisely the proportion that ensures these weighted samples produceexpectations which are concordant with expectations under π(·) when used in (2.4).This is laid out precisely in Algorithm2.3.Algorithm 2.3 Importance sampling algorithm1: procedure Importance sampling(π(·), π̃(·), c) Generate random sample from dis-tribution with un-normalised densityπ(·)2: x ∼ π̃(·) Propose sample3: w ← π(x)π̃(x)4: return (x,w).5: end procedureTo see that this weighting has the desired effect, consider the expectation whichis our objective. We first consider the situation where both π and π̃ are normalised:Eπ [ f (X)] =∫f (x)π(x) dx=∫f (x)π̃(x)π̃(x)π(x) dx multiply and divide by π̃(x)=∫(f (x)π(x)π̃(x))π̃(x) dx= Eπ̃[f (X)π(X)π̃(X)].That is, we use samples directly from π̃(·), and instead adjust (2.4) to target theexpectation of the same functional under π(·).μ ∫f (x)π(x) dx ≈ 1nn∑j=1f (x j )π(x j )π̃(x j )︸︷︷︸=wj= 1nn∑j=1f (x j )wj μ̂, (2.9)where now x j ∼ π̃(·).Some care is required, because although this estimator remains unbiased, thevariance is no longer going to be the same as the usual Monte Carlo variance wherex j ∼ π(·). Indeed, nowVar(μ̂) = σ2π̃nwhere σ 2π̃ =∫( f (x)π(x) − μπ̃(x))2π̃(x)dx,which can be empirically estimated from the importance samples using,2 Sampling from Complex Probability Distributions: A Monte … 29σ̂ 2π̃ =1nn∑j=1(f (x j )wj − μ̂)2. (2.10)As such,σ 2π̃(or its empirical estimate σ̂ 2π̃) provide a guide towhenwehave a ‘good’importance sampling algorithm, since with π̃(·) fixed the only option to improve theestimate is to increase the sample size n.Indeed, it can be shown [15] that the theoretically optimal proposal distributionwhich minimises the estimator variance isπ̃(x)opt = | f (x)|π(x)∫| f (x)|π(x) dx .In particular, note that this implies that importance sampling can achieve super-efficiency whereby it results in lower variance even than sampling directly fromπ(·) when f (x) = x . Specifically, if f (x) ≥ 0 ∀x then this proposal results in azero-variance estimator! Of course, in practice we cannot usually sample from andevaluate this optimal proposal, since it is at least as difficult as the original problemwewere attempting to solve. However, even though these optimal proposals are oftenunusable, they provide guidance towards the form of a good proposal for any givenimportance sampling problem.2.3.3.1 Self-normalising WeightsThe option to use rejection samplingwith un-normalised densities is very helpful (e.g.in Bayesian settingswhere the normalising constant is often unknown).We can retainthis advantage with importance sampling by using self-normalising weights. Thealgorithm to generate the weights remains as in Algorithm2.3, but the computationof the estimator in (2.9) changes. The self-normalised version, rather than dividingby n, uses the sum of the weights,μ̂ ∑nj=1 f (x j )wj∑nj=1 wj,thereby ensuring cancellation of the unknown normalising constant from the targetand/or proposal distributions in the weights.However, it is important to note that this estimator is no longer unbiased, thoughasymptotically it is. Additionally, the variance of this estimator is more complicatedhaving only approximate form. An approximate estimate can be computed using,30 L. J. M. AslettVar(μ̂) ≈ σ̂2π̃nwhere σ̂ 2π̃ =n∑j=1w2j(f (x j ) − μ̂)2and wj =wj∑ni=1 wi.Finally, the theoretically optimal (but usually unusable) proposal in the self-normalised weight case isπ̃(x)opt ∝ | f (x) − μ|π(x).In both regular and self-normalised weight settings, one can then compute appro-priate confidence intervals in the usual manner.2.3.3.2 DiagnosticsAdditional care is required in the application of importance samplingwhen comparedto using iid samples from the distribution of interest. In particular, because importancesamplinguses aweighted collection of samples, it is not uncommon tobe in a situationwhere a small number of samples with large weight dominate the estimate, so thatsimply having many importance samples does not equate to good estimation overall.A common diagnostic for potential weight imbalance is derived by equating thevariance of aweighted importance sampling approach to the standard iidMonteCarlovariance for an average computed using a fixed but unknown sample size ne. Uponsimple algebraic rearrangement one may then solve for ne, the so-called effectivesample size. This informally corresponds to the size of iid Monte Carlo sample onewould expect to need to attain the same variance achieved via this importance sample,so that a low value indicates poor weight behaviour (since that corresponds to fewiid samples).Var(∑ni=1 f (xi )wi∑ni=1 wi)= σ2ne=⇒ ne = nw̄2w2,wherew̄2 =⎛⎝1nn∑j=1wj⎞⎠2and w2 = 1nn∑j=1w2j .The reason to use such a diagnostic and not simply rely on the empirical varianceestimates above is that they are themselves based on the sampling procedure andtherefore may be poor estimates too.Finally, it is critical to note that although small ne does diagnose a problem withimportance sampling, it is not necessarily true that large ne means everything is ok:2 Sampling from Complex Probability Distributions: A Monte … 31it is, for example, entirely possible that the sampler has missed whole regions of highprobability.2.3.3.3 ExampleConsider the toy problem of computingP(X > 3.09)when X is a univariate standardNormal random variable. The R-language [26] computes the distribution functionof the standard Normal to at least 18 significant digits using a rational Chebyshevapproximation [6] (see /src/nmath/pnorm.c) and we know the true answerto be 0.001001 (4 sf). However, if the reader suspends disbelief and imagines thatwe cannot accurately compute the distribution function, but can only compute theNormal density and draw random realisations from it, then evaluation of the aboveprobability might be approximated using Monte Carlo methods instead (given theunbounded support and extreme tail location this may be preferred to numericalintegration).Since we are assuming the ability to generate random realisations from theNormal distribution, a standard Monte Carlo approach would draw many samplesxi ∼ N(0, 1) and computeP(X > 3.09) = E [I[3.09,∞)(X)] = 1nn∑j=1I[3.09,∞)(x j )However, this will require many samples to achieve an accurate estimate of this tailprobability.In contrast, still only using simulations from a Normal distribution, we may electto use importance sampling with a proposal N(m, 1) for some choice m. We knowthe fully normalised density of a Normal distribution and therefore will be using theestimator (2.9) with associated single sample variance which can be approximatedusing (2.10). Therefore, to select m, we perform a small grid search over possibleproposals, computing σ̂ 2π̃each time, to find a good choice. This results in Fig. 2.4,showing that a proposal N(3.25, 1) is a good choice.A further final run of n = 100,000 samples renders an estimate μ̂ = 0.001002 (4sf). The same pseudo-random number stream using standardMonte Carlo renders anestimate 0.001140 (4 sf), which is a relative error 163× larger than the importancesampling estimate. To demonstrate this is not a ‘fluke’ result, we continue to repeatboth importance sampling and standard Monte Carlo estimation with runs of sizen = 100,000 and plot the density of estimates of P(X > 3.09) in Fig. 2.5.Note that Fig. 2.5 demonstrates how much more accurate importance samplingis for the same sample size n = 100,000 when computing this event probabilitycompared to standardMonte Carlo. One may reasonably object that we have ignoredthe 25 pilot runs of n =100,000 importance samples used to select m = 3.25, sothat the total computational effort expended on importance sampling was at least26× that of standard Monte Carlo. However, it is a simple calculation to determine32 L. J. M. AslettFig. 2.4 The estimate of σ̂ 2π̃using (2.10) for 25 different values ofm in theNormal proposalN(m, 1)used in an importance sampling estimator of P(X > 3.09), where X ∼ N(0, 1). Each estimate ofσ̂ 2π̃is based on n = 100, 000 samples. m varies on an equally spaced grid from 1.5 to 4.5. Theminimum is at m = 3.25Fig. 2.5 A total of 10, 000 runs of both importance sampling and standardMonteCarlo, each of sizen = 100,000. Each run was used to compute the estimate of P(X > 3.09) where X ∼ N(0, 1) anda kernel density plot of these estimates produced (importance sampling = dashed, standard MonteCarlo = solid). The vertical line is the ground truth computed usingpnorm(3.09, lower.tail= FALSE). The same pseudo-random number stream was used for each method to ensure a faircomparisonthat based on the standard deviation of the samples used to generate Fig. 2.5 andthe√n convergence of Monte Carlo, that it would require a standard Monte Carlosample of size n = 295×100,000 to achieve the same accuracy profile as importancesampling. Therefore, even accounting for the pilot computational effort to select aproposal distribution, there is a substantial benefit to using importance sampling.2.4 Further ReadingA textbook length introduction with a solid emphasis on implementation details inR can be found in [28]. The same authors have a more advanced textbook going intothe theoretical aspects more deeply [27]. Both these books also introduce Markov2 Sampling from Complex Probability Distributions: A Monte … 33Chain Monte Carlo methods, which are often used in practice in high dimensionalproblems. A nice tutorial paper introduction to MCMC is [1] and [4] is an excellentcollection of chapters on the topic.A classic Monte Carlo text is [10], which is now freely (and legally) availableonline and contains many results not easily found elsewhere.Although standard Monte Carlo and Markov Chain Monte Carlo arguably rep-resent the mainstay of most practical uses of Monte Carlo, there are an array ofadvancedmethodswhich are particularlywell suited to different settings. Someexcel-lent review texts as jumping off points to explore some of these include [8] (Sequen-tial Monte Carlo), [19, 20] (Approximate Bayesian Computation), [14] (Multi-LevelMonte Carlo), and [2] (MLMC in engineering reliability).Excellent software choices for practically performing inference in complex mod-els via sampling methods includes Stan [5] and Birch [23].DedicationThis chapter is dedicated to the memory of Brett Houlding (1982–2019). The tragicnews of Brett’s passing was received while I was in the act of writing this chapter.He will be sorely missed.References1. C. Andrieu, N. de Freitas, A. Doucet, andM. I. Jordan. An introduction toMCMC for MachineLearning. Machine Learning, 50:5–43, 2003.2. L. J.M.Aslett, T. Nagapetyan, and S. J. Vollmer.MultilevelMonte Carlo for Reliability Theory.Reliability Engineering & System Safety, 165:188–196, 2017.3. G. E. P. Box and M. E. Muller. A note on the generation of Normal deviates. Annals ofMathematical Statistics, 29(2):610–611, 1958.4. S. Brooks, A. Gelman, G. Jones, and X.-L. Meng, editors. Handbook of Markov Chain MonteCarlo. Chapman and Hall/CRC, 2011.5. B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker,J. Guo, P. Li, andA. Riddell. Stan: A probabilistic programming language. Journal of StatisticalSoftware, 76(1):1–32, 2017. https://mc-stan.org/.6. W. J. Cody. Rational Chebyshev approximations for the error function. Mathematics of Com-putation, 23(107):631–637, 1969.7. D. R. Cox. Applied statistics: a review. The Annals of Applied Statistics, 1(1):1–16, 2007.8. C. Dai, J. Heng, P. E. Jacob, andN.Whiteley. An invitation to sequentialMonte Carlo samplers,2020. arXiv:2007.11936 [stat.CO].9. H. Dai, M. Pollock, and G. O. Roberts. Monte Carlo Fusion. arXiv preprint arXiv:1901.00139,2019.10. L. Devroye. Non-Uniform Random Variate Generation. Springer-Verlag, New York, 1986.Available at http://luc.devroye.org/rnbookindex.html.11. G. Feng, E. Patelli, M. Beer, and F. P. A. Coolen. Imprecise system reliability and componentimportance based on survival signature.Reliability Engineering&SystemSafety, 150:116–125,2016.34 L. J. M. Aslett12. N. Friel, R. Rastelli, J. Wyse, and A. E. Raftery. Interlocking directorates in Irish companiesusing a latent space model for bipartite networks. Proceedings of the National Academy ofSciences, 113(24):6629–6634, 2016.13. J. E. Gentle. Random number generation and Monte Carlo methods. Springer Science &Business Media, 2006.14. M. B. Giles. Multilevel Monte Carlo methods. Acta Numerica, 24:259–328, 2015.15. G. Goertzel. Quota sampling and importance functions in stochastic solution of particle prob-lems, volume 2793. US Atomic Energy Commission, Technical Information Division, 1950.16. J. Haslett, M. Whiley, S. Bhattacharya, M. Salter-Townshend, S. P. Wilson, J. R. M. Allen,B. Huntley, and F. J. G. Mitchell. Bayesian palaeoclimate reconstruction. Journal of the RoyalStatistical Society: Series A (Statistics in Society), 169(3):395–438, 2006.17. B. Houlding and S. P. Wilson. Considerations on the UK re-arrest hazard data analysis. Law,Probability & Risk, 10(4):303–327, 2011.18. A. J. Kinderman and J. G. Ramage. Computer generation of Normal random variables. Journalof the American Statistical Association, 71(356):893–896, 1976.19. J. M. Marin, P. Pudlo, C. P. Robert, and R. J. Ryder. Approximate Bayesian computationalmethods. Statistics and Computing, 22(6):1167–1180, 2012.20. P. Marjoram, J. Molitor, V. Plagnol, and S. Tavaré. Markov chain Monte Carlo without likeli-hoods. PNAS, 100(26):15324–15328, 2003.21. G. Marsaglia. Random number generators. Journal of Modern Applied Statistical Methods,2(1):2, 2003.22. M. Matsumoto and T. Nishimura. Mersenne twister: a 623-dimensionally equidistributed uni-form pseudo-random number generator. ACM Transactions on Modeling and Computer Sim-ulation (TOMACS), 8(1):3–30, 1998.23. L.M.Murray and T. B. Schön. Automated learningwith a probabilistic programming language:Birch. Annual Reviews in Control, 46:29–43, 2018. https://www.birch.sh/.24. A. B. Owen. Monte Carlo theory, methods and examples. 2013.25. M. Pollock, A.M. Johansen, and G. O. Roberts. On the exact and ε-strong simulation of (jump)diffusions. Bernoulli, 22(2):794–856, 2016.26. R Core Team. R: A Language and Environment for Statistical Computing. R Foundation forStatistical Computing, Vienna, Austria, 2019.27. C. P. Robert and G. Casella. Monte Carlo Statistical Methods. Springer-Verlag New York,2004.28. C. P. Robert and G. Casella. Introducing Monte Carlo Methods with R. Springer, New York,NY, 2010.29. S. P. Wilson and M. J. Costello. Predicting future discoveries of European marine species byusing a non-homogeneous renewal process. Journal of the Royal Statistical Society: Series C,54(5):897–918, 2005.30. E. Zio, P. Baraldi, and E. Patelli. Assessment of the availability of an offshore installation byMonte Carlo simulation. International Journal of Pressure Vessels and Piping, 83(4):312–320,2006.2 Sampling from Complex Probability Distributions: A Monte … 35Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 3Introduction to the Theoryof Imprecise ProbabilityErik QuaeghebeurAbstract The theory of imprecise probability is a generalization of classical ‘pre-cise’ probability theory that allows modeling imprecision and indecision. This is apractical advantage in situations where a unique precise uncertainty model cannotbe justified. This arises, for example, when there is a relatively small amount ofdata available to learn the uncertainty model or when the model’s structure cannotbe defined uniquely. The tools the theory provides make it possible to draw conclu-sions and make decisions that correctly reflect the limited information or knowledgeavailable for the uncertainty modeling task. This extra expressivity however oftenimplies a higher computational burden. The goal of this chapter is to primarily giveyou the necessary knowledge to be able to read literature that makes use of the theoryof imprecise probability. A secondary goal is to provide the insight needed to useimprecise probabilities in your own research. To achieve the goals, we present theessential concepts and techniques from the theory, as well as give a less in-depthoverview of the various specific uncertainty models used. Throughout, examples areused to make things concrete. We build on the assumed basic knowledge of classicalprobability theory.3.1 IntroductionThe theory of imprecise probability is a generalization of classical ‘precise’ proba-bility theory that allows modeling imprecision and indecision. Why is such a theorynecessary? Because inmany practical applications a lack of information—e.g., aboutmodel parameters—and paucity of data—especially if we also consider conditionalmodels—make it impossible to create a reliable model.For example, consider a Bayesian context where a so-called prior probabilitydistribution must be chosen as part of the modeling effort. The lack of information(This chapter was written while at Delft University of Technology)E. Quaeghebeur (B)Eindhoven University of Technology, Eindhoven, The Netherlands© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_33738 E. Quaeghebeurmay make it difficult to determine the type of the prior distribution, let alone itsparameters. Then, even if we assume some prior has been chosen—e.g., a normalone—in a somewhat arbitrary way, a paucity of data will make the parameters of theposterior—updated—distribution depend to a large degree on the prior’s somewhatarbitrary parameters. The consequence is that conclusions drawn from the posteriorare unreliable and decisions based on it somewhat arbitrary.The theory of imprecise probability provides us with a set of tools for dealing withthe problem described above. For the example above, instead of choosing a singleprior distribution, a whole set of priors is used, one that is large enough to sufficientlyreduce or even eliminate the arbitrariness of this modeling step. The consequence isthat conclusions drawn from an imprecise probabilistic model are more reliable bybeing less committal—more vague, if youwish; somewould say ‘more honest’—andthat decisions based on it allow for indecision.In this chapter, we will go over the basic concepts of the theory of imprecise prob-ability theory. Therefore, we will consider ‘small’ problems, with finite possibilityspaces. However, the theory can be applied to infinite—countable and uncountable—possibility spaces as well. Also, only the basics of more advanced topics such asconditioning will be touched upon. But, and this is the chapter’s goal, after havingunderstood the material we do treat, the imprecise probability literature should havebecome substantially more accessible. Good extensive general treatments are avail-able [2, 16, 20] and the proceedings of the ISIPTA conferences provide an extensiveselection of papers developing imprecise probability theory or applying it [1, 3, 4,6–12].Concretely, we start with a discussion of the fundamental concepts in Sect. 3.2.This is done in terms of themore basic notion of sets of acceptable gambles. Probabil-ities only appear thereafter, in Sect. 3.3, together with the related notion of prevision(expectation). The connection with sets of probabilities is made next, in Sect. 3.4.Then we touch upon conditioning, in Sect. 3.5, and before closing add some remarksabout continuous possibility spaces, in Sect. 3.6. Throughout we will spend ampletime on a running example to illustrate the theory that is introduced.3.2 Fundamental ConceptsIn this section, we introduce the fundamental concepts of the theory of impreciseprobability [18] [20, §3.7]. First, in Sect. 3.2.1, we get started with some basicconcepts. Then, in Sect. 3.2.2, we list and discuss the coherence criteria on whichthe whole theory is built.3 Introduction to the Theory of Imprecise Probability 393.2.1 Basic ConceptsConsider an agent reasoning about an experiment with an uncertain outcome. Thisexperiment is modeled using a possibility space—a set—X of outcomes x . Nowconsider the linear space L = X → R of real-valued functions over the outcomes.We view these functions as gambles because they give a value, seen as a payoff, foreach outcome and because the outcome is uncertain and therefore the payoff is aswell.A special class of gambles are the outcome indicators 1x or subset indicators 1B ,which take the value one on that outcome or subset and zero elsewhere.The agent can then express her uncertainty by specifying a set of gambles, calledan assessment A, that she considers acceptable. Starting from such an assessment,she can reason about other gambles and decide whether she should also accept themor not. If she were to do this for all gambles, then the natural extension E of herassessment would be the set of all acceptable gambles. To reason in a principled way,she needs some guiding criteria; these are the next section’s topic.Let us now introduce our running example:Wiske and Yoko Tsuno want to bet on Belgium vs. JapanGiven a sports match between Belgium and Japan, there is uncertainty about whichcountry’s team will win. So we consider the possibility space{be, jp}. There are toagents—gamblers—:Wiske and Yoko Tsuno, two comic book heroines. Each has anassessment consisting of a single gamble that they find acceptable:• Wiske accepts losing 5 coins if Japan wins for the opportunity to win 1 coin ifBelgium wins; soAW = {1be − 5 · 1jp}.• Yoko Tsuno accepts losing 4 coins if belgium wins for the opportunity to win1 coin if Japan wins; soAY = {−4 · 1be + 1jp}.The heroines are also discussing joining forces and forming a betting pool. Thepools they consider are• ‘Simple’, formed by combining their assessments; soASP = {1be − 5 · 1jp,−4 · 1be + 1jp} .• ‘Empty’ in case of disagreement, without any acceptable gambles; so AEP = ∅.3.2.2 CoherenceIn the theory of imprecise probabilities, the classical rationality criteria used forreasoning about assessments are called coherence criteria. These are typically for-mulated as four rules that should apply to any gambles f and g. (There are different40 E. Quaeghebeurvariants in the literature, but the differences are not relevant in this introductory text.)We divide the criteria into two classes.Constructive State how to generate acceptable gambles from the assessment:Positive scaling If f is acceptable and λ > 0, then λ · f is acceptable.Addition If f and g are acceptable, then f + g is acceptable.Background State which gambles are always or never acceptable:Accepting gain If f is nonnegative for all outcomes, then f is acceptable.Avoiding sure loss If g is negative for all outcomes, then g is not acceptable.These criteria are quite broadly seen as reasonable, under the assumption that thepayoffs are ‘not too large’.The last criterion, ‘Avoiding sure loss’, puts a constraint on what is consideredcoherent; if it is violated, we say that an assessment incurs sure loss. The first threerules can be used to create an explicit expression for the natural extension:E ={∑f ∈K λ f · f : K A ∪ { f ∈ L : f ≥ 0} and (∀ f ∈ K : λ f ≥ 0)},where denotes the finite subset relation. Then E is the smallest convex cone ofgambles encompassing the assessment A and the nonnegative gambles—includingthe zero gamble.Let us apply the natural extension to our running example:The natural extensions of Wiske, Yoko Tsuno, and the betting poolsFor our finite possibility space,{ f ∈ L : f ≥ 0} = {∑x∈X μx · 1x : (∀x ∈ X : μx ≥ 0)}So, with λA, μx ≥ 0 for all outcomes x and agent identifiers A, we get the followingexpressions that characterize the natural extensions:WiskeλW · (1be − 5 · 1jp) + μbe · 1be + μjp · 1jp= (λW + μbe) · 1be + (−5 · λW + μjp) · 1jp,Yoko TsunoλY · (−4 · 1be + 1jp) + μbe · 1be + μjp · 1jp= (−4 · λY + μbe) · 1be + (λY + μjp) · 1jp,Simple pool (λW − 4 · λY + μbe) · 1be + (−5 · λW + λY + μjp) · 1jp,Empty pool μbe · 1be + μjp · 1jp.To check whether the natural extension incurs sure loss, we must check whether thecoefficients of 1be and 1jp can become negative at the same time. Only the simple poolincurs sure loss; e.g., fill in λW = λY = 1 and μbe = μjp = 0 to convince yourself.(Convince yourself as well that the others avoid sure loss indeed.)3 Introduction to the Theory of Imprecise Probability 413.3 Previsions and ProbabilitiesIn this section, we move from modeling uncertainty using sets of acceptable gam-bles to the more familiar language of expectation—or, synonymously, prevision—and probability [17]. We first transition from acceptable gambles to previsions inSect. 3.3.1 [18, §1.6.3] [17, §2.2] and in a second step, in Sect. 3.3.2, give the con-nection to probabilities [20, §2.6]. Next, in Sect. 3.3.3, we consider assessments interms of previsions and what the other fundamental concepts of Sect. 3.2 then looklike [17, §2.2.1, §2.2.4] [20, §2.4–5, §3.1]. Finally, in Sect. 3.3.4, we consider theimportant special case of assessments in terms of previsions defined on a linear spaceof gambles [17, §2.2.1] [20, §2.3.2–6].3.3.1 Previsions as Prices for GamblesBefore we start: ‘prevision’ is in much of the imprecise probability literature usedas a synonym for ‘expectation’; we here follow that tradition.Now, how do we get an agent’s previsions for a gamble—equivalently: expecta-tion of a random variable—given that we know the agent’s assessment as a set ofacceptable gambles A? We first define a price to be a constant gamble and iden-tify this constant gamble with its constant payoff value. Then we define the agent’sprevisions as specific types of acceptable prices:• The lower prevision P( f ) is the supremum acceptable buying price of f :P( f ) = sup {ν ∈ R : f − ν ∈ E} .• The upper prevision P( f ) is the infimum acceptable selling price of f :P( f ) = inf {κ ∈ R : κ − f ∈ E} .If E is coherent, then P and P are also called coherent. There is a conjugacy relationbetween coherent lower and upper previsions: P( f ) = −P(− f ). It allows us towork in terms of either type of prevision; we will mainly use the lower one.In case P( f ) = P( f ), then P( f ) = P( f ) is the called the (precise) prevision ofthe gamble f .3.3.2 Probabilities as Previsions of Indicator GamblesNow that we have definitions for lower and upper previsions, we can derive probabil-ities from those. For classical probability, we have that the probability of an event—a42 E. Quaeghebeursubset B of the possibility space X—is the prevision of the indicator for that event.For lower and upper previsions, we get:• The lower probability: P(B) = P(1B).• The upper probability: P(B) = P(1B).Notice that we reuse the same symbol for the prevision and probability functions,as is common in the literature. As long as the nature of the argument—gamble orevent—is clear, this does not cause ambiguity. If P and P are coherent as previsions,then so are they as probabilities. Also the conjugacy relationship can be translatedto coherent lower and upper probabilities; let Bc = X \ B, thenP(B) = P(1B) = P(1 − 1Bc ) = −P(−1 + 1Bc ) = 1 − P(1Bc ) = 1 − P(Bc).In case P(B) = P(B), then P(B) = P(B) is called the (precise) probability ofB.To make the definitions for lower and upper previsions and probabilities concrete,let us apply them to our running example:Lower and upper probabilities for all events and agentsWe work out the calculation of Wiske’s lower probability that Belgium will win.PW(be) =PW(1be) (def. lower probability)= sup {ν ∈ R : 1be − ν ∈ EW}(def. lower prevision)= sup{ν ∈ R :[1 − ν0 − ν]=[λW + μbe−5 · λW + μjp], λW ≥ 0, μbe ≥ 0, μjp ≥ 0}(write out natural extension EW of AW)= sup{5 · λW − μjp : 1 − 5 · λW + μjp = λW + μbe, λW ≥ 0, μbe ≥ 0, μjp ≥ 0}(eliminate ν)= sup{5 · λW − μjp : λW = 16 (1 + μjp − μbe), λW ≥ 0, μbe ≥ 0, μjp ≥ 0}(solve constraint for λW)= sup{56 − 16μjp − 56μbe : 1 + μjp ≥ μbe, μbe ≥ 0, μjp ≥ 0}(eliminate λW)= 56(feasible solution μbe = 0, μjp = 0 maximizes expression)Do the calculations also for the other agents and Japan. Then apply conjugacy to findthe following table of lower and upper probabilities:Agents P(be) , P(be) P(jp) , P(jp) NoteWiske 5/6 , 1 0 , 1/6 Will not bet against BelgiumYoko Tsuno 0 , 1/5 4/5 , 1 Will not bet against JapanSimple pool +∞ , −∞ +∞ , −∞ Sure loss, so absurd boundsEmpty pool 0 , 1 0 , 1 So-called vacuous model3 Introduction to the Theory of Imprecise Probability 43While in the above example the calculation of the lower prevision can be done byhand, in general it realistically requires a linear program solver.3.3.3 Assessments of Lower PrevisionsUp until now, we assumed a set of acceptable gamblesA—an agent’s assessment—to be given. But often the agent will directly specify lower and upper probabilitiesor previsions, e.g., as bounds on precise probabilities and previsions. However, thecoherence criteria and expression for the natural extension are based on having aset of acceptable gambles. In this section we will provide expressions based on anassessment specified as lower prevision values for gambles in a given set K .The approach is to derive an assessment as a set of acceptable gambles A fromthese lower prevision. Irrespective of what its natural extension E actually looks like,it follows from the definition of the lower prevision as a supremum acceptable buyingprice that0 ≤ sup {ν − P( f ) : ν ∈ R ∧ f − ν ∈ E ⊇ A}= sup {κ ∈ R : f − (κ + P( f )) ∈ E} = sup {κ ∈ R : ( f − P( f )) − κ ∈ E} .This implies that f − P( f ) + ε ∈ E for any ε > 0, because of coherence.We cannottake ε = 0, because the corresponding so-called marginal gamble f − P( f ) is notincluded in E in general, as the supremum value κ = 0 is not necessarily attainedinside the set. We therefore take A = ⋃ f ∈K{f − P( f ) + ε : ε > 0}.We can then apply the theory described above to this assessmentA. This leads tothe following nontrivial results for a lower prevision P defined on a set of gamblesK :• It avoids sure loss if and only if for all n ≥ 0 and fk ∈ K it holds thatsupx∈Xn∑k=1(fk(x) − P( fk)) ≥ 0.• It is coherent if and only if for all n,m ≥ 0 and fk ∈ K it holds thatsupx∈X(n∑k=1(fk(x) − P( fk)) − m · ( f0 − P( f0)))≥ 0.• Its natural extension to any gamble f in L isE( f ) = sup⎧⎨⎩infx∈X⎧⎨⎩f (x) −n∑k=1λk ·(fk(x) − P( fk))⎫⎬⎭: n ≥ 0, fk ∈ K, λk ≥ 0⎫⎬⎭.44 E. Quaeghebeur3.3.4 Working on Linear Spaces of GamblesThe coherence criterion for lower previsions on an arbitrary spaceK of gambles wegave in the preceding section is quite involved.However, in caseK is a linear space ofgambles, this criterion becomes considerably simpler. Namely, a lower prevision Pmust then satisfy the following criteria for all gambles f and g in K and λ > 0:Accepting sure gains P( f ) ≥ inf f ,Super-linearity P( f + g) ≥ P( f ) + P(g),Positive homogeneity P(λ f ) = λ · P( f ).Expressed for upper previsions P , these coherence criteria are very similar:Accepting sure gains P( f ) ≤ sup f ,Sub-linearity P( f + g) ≤ P( f ) + P(g),Positive homogeneity P(λ f ) = λ · P( f ).From the coherence criteria, many useful properties can be derived for a coherentlower prevision P and its conjugate upper prevision P . We provide a number of keyones, which hold for all gambles f and g inK and μ ∈ R; P denotes either P or P:Upper dominates lower P( f ) ≥ P( f ),Constants P(μ) = μ,Constant additivity P( f + μ) = P( f ) + μ,Gamble dominance if f ≥ g + μ then P( f ) ≥ P(g) + μ,Mixed sub/super-additivity P( f + g) ≤ P( f ) + P(g) ≤ P( f + g).3.4 Sets of ProbabilitiesIn Sect. 3.2 we modeled uncertainty using a set of acceptable gambles. In Sect. 3.3we showed how this can also be done in terms of lower or upper previsions (orprobabilities). In this section, we add a third representation, one using credal sets—sets of precise previsions [17, §2.2.2], [18, §1.6.2]. In Sect. 3.4.1 we show how toderive the credal set corresponding to a given lower prevision. In Sect. 3.4.2 we gothe other direction and show how to go from a credal set to lower prevision values[20, §3.3].3.4.1 From Lower Previsions to Credal SetsA credal set is a subset of the set of all precise previsionsP. (For possibility spaces Bdifferent fromX, we writePB .) This set is convex, meaning that any convex mixture3 Introduction to the Theory of Imprecise Probability 45of precise previsions is again a precise prevision.Because of this, a gamble’s previsionis a linear function over this space. A lower—and upper—prevision can be seenas providing a bound on the value of the precise prevision for that gamble andthereby represent a linear constraint on the precise previsions. So the credal set Mcorresponding to a lower prevision P defined on a set of gambles K is the subsetof P satisfying this constraint for all gambles in K :M =⋂f ∈K{P ∈ P : P( f ) ≥ P( f )} .Being defined as such an intersection, such credal sets are closed and convex.The rationality criteria for a lower prevision P we encountered before can alsobe expressed using its corresponding credal setM:• P incurs sure loss if and only if M is equal to the empty set.• P is coherent if and only if all constraints are ‘tight’, i.e., if there exists a P inMsuch that P( f ) = P( f ) for all f in K .Let us make the concept of a credal set concrete using our running example:Yoko Tsuno’s credal setFor a finite possibility space such as the one of our running example, a precise previ-sion P can be defined completely by the corresponding probability mass function pdefined by px = P({x}) for x in X ={be, jp}. The set of all precise previsionscan therefore be represented by the probability simplex—the set of all probabilitymass functions—on X. This set and the example probability mass function ( 12 , 12 )is shown below left. Below right, we illustrate how Yoko Tsuno’s lower previ-sion PY(jp) = 45 generates the credal set MY: The gamble 1jp as a linear functionover the simplex is shown as an inclined line. This linear relationship between p—equivalently, the corresponding prevision Pp—and Pp(1jp) = Pp(jp) transforms thebounds 45 ≤ Pp(jp) ≤ 1 intoMY.jpbep = (pbe, pjp)( 12 ,12 )(1, 0) (0, 1)jpbe01Pp(jp) = 0 · pbe + 1 · pjpPY(jp) = 45( 15 ,45 )MYThe set of extreme pointsM∗Y ofMY as probability mass functions is{( 15 ,45 ), (0, 1)}.46 E. Quaeghebeur3.4.2 From Credal Sets to Lower PrevisionsNowwe assume that the agent’s credal setM is given. Most generally this can be anyset of precise previsions, i.e., any subset of P. Often, to ensure equivalence betweencoherent lower previsions and non-empty credal sets, they are required to be closedand convex. In that case, a credal set is determined completely by its set of extremepoints M∗ in the sense that all other elements are convex mixtures of these.To determine the lower prevision corresponding to any credal set, we determineits value for each gamble f of interest using the lower envelope theorem:P( f ) = min {Pp( f ) : p ∈ M} = min {Pp( f ) : p ∈ M∗}.Let us again use the running example to provide a feeling for what this all means:A credal set for the empty pool facing penaltiesConsider the empty pool. Because its assessment is empty, its credal setMEP is thetrivial one corresponding to all probability mass functions onX = {be, jp}. Now weadd an extra element to the possibility space, ‘Penalties’. Below left we show MEPembedded in the corresponding larger probability simplex. Wiske and Yoko Tsunodecide to add the uniform probability mass function to it. Below right, you see theconvex hullMEUP of this extra probability mass function and the original credal set.P(enalties)jpbe(1, 0, 0) (0, 1, 0)( 13 ,13 ,13 )MEPP(enalties)jpbeMEUPPp(p) = 13If wewant to calculate lower and upper prevision values, we can here use the extremepoint version of the lower and—similar—upper envelope theorem. For example, forthe pool’s upper probability for Penalties:PEP(p) = PEP(1p) = max{p(0, 0, 1) : p ∈{(1, 0, 0), (0, 1, 0), ( 13 ,13 ,13 )}}= 13.To make it explicit where this maximum is achieved, we above right show the lineof probability mass functions p such that Pp(p) = 13 .3 Introduction to the Theory of Imprecise Probability 473.5 Basics of ConditioningConditioning an uncertainty model is the act of restricting attention to a subset Bof the possibility space. It is often used to update an uncertainty model after havingobserved the event B [20, §6.1].In the theory of imprecise probability, conditioning is a specific case of naturalextension [17, §2.3.3], [20, §6.4.1]. In terms of acceptable gambles, conditioningon B corresponds to restricting the space of gambles to those that are zero outside B[18, §1.3.3]. For lower previsions, this translates to the following conditioning rulefor all gambles f in L:E( f | B) ={inf x∈B f (x) if P(B) = 0,max{μ ∈ R : P(1B( f − μ) = 0)}if P(B) > 0.Conditioning a credal set M corresponds to taking the credal set M|B formed byconditioning each of the precise previsions inM:M|B ={PB if ∃P ∈ M : P(B) = 0,{P(· | B) : P ∈ M} if ∀P ∈ M : P(B) > 0.These rules based on natural extension give vacuous conditionals whenever thelower probability of the conditioning event is zero. Regular extension is a less impre-cise updating rule [17, §2.3.4], [18, §1.6.6], [20, App. J]: In credal set terms, itremoves those precise previsions P such that P(B) = 0 fromM.Let us apply the conditioning rules discussed here to our running example:Conditioning the empty-uniform pool’s credal setWe condition the empty-uniform pool’s credal set on {jp, p}, i.e., Belgium notwinning in regular time. Further down on the left, we show what happens if weapply natural extension: the conditional model is vacuous because P(1,0,0)({jp, p}) =P(1,0,0)(1{jp,p}) = (1, 0, 0)(0, 1, 1) = 0. Further down on the right, we apply natu-ral extension and therefore remove P(1,0,0) fromMEUP; this results in a non-vacuousconditional credal set.48 E. QuaeghebeurP(enalties)jpbe(1, 0, 0) (0, 1, 0)( 13 ,13 ,13 )MEUPMEUP| {jp, p}P(enalties)jpbe(0, 1, 0)( 13 ,13 ,13 )M′EUPM′EUP| {jp, p}3.6 Remarks About Infinite Possibility SpacesThe theory we presented is also applicable to denumerable and continuous pos-sibility spaces with some technical amendments to the coherence criteria and byconsidering only bounded gambles. However, the running example was based onfinite possibility spaces, so no feeling was created for applications with infinite pos-sibility spaces. Therefore we here give some remarks about imprecise probabilisticuncertainty models on continuous possibility spaces:• They are mostly defined using credal sets whose extreme points are parametricdistributionswhere the parameters vary in a set. A prime example are the impreciseDirichlet model [21] and its generalizations [19].• They are also commonly defined using probability mass assignments to subsets ofthe possibility space. This is in some way a reduction to the finite case. Examplesare belief functions [13, §5.2.1.1], some P-boxes [14, §4.6.4], and NPI models [5,§7.6].• Furthermore,modelswhich bound some specific description of a precise prevision,such as cumulative distribution functions andprobability density functions, are alsopopular in some domains. The extreme points of their credal set are, however, notknown. General P-boxes [15] and lower and upper density functions [20, §4.6.3]are examples of this class.• Calculating lower and upper previsions—i.e., performing natural extension—caneasily become difficult optimization problems, so this should be a key considera-tion when choosing a specific type of model.3 Introduction to the Theory of Imprecise Probability 493.7 ConclusionThis introduction to the theory of imprecise probability has prepared you for access-ing the broader literature on this topic and its applications. For those that wish toapply imprecise probabilistic techniques, this text only provides the first step: Youshould dive into the literature and contact experts to obtain the necessary knowledgeand feedback. The references of this chapter and their authors or editors provide astarting point for that.Acknowledgements I would like to thank Frank Coolen for giving me the opportunity to give alecture at the 2018 UTOPIAE Training School and all the students for their active participation.Thanks are also due to theWind Energy group at TUDelft, which gaveme the freedom to contributeto the Training School, despite the demands of my regular duties.References1. Alessandro Antonucci, Giorgio Corani, Inés Couso, and Sébastien Destercke, editors.ISIPTA ’17: Proceedings of the Tenth International Symposium on Imprecise Probability:Theories and Applications, volume 62 of Proceedings of Machine Learning Research, 2017.2. Thomas Augustin, Frank P. A. Coolen, Gert de Cooman, and Matthias C. M. Troffaes, editors.Introduction to Imprecise Probabilities. Wiley Series in Probability and Statistics.Wiley, 2014.3. Thomas Augustin, Frank P. A. Coolen, Serafín Moral, and Matthias C. M. Troffaes, editors.ISIPTA ’09: Proceedings of the Sixth International Symposium on Imprecise Probability: The-ories and Applications, 2009.4. ThomasAugustin, SerenaDoria, EnriqueMiranda, and Erik Quaeghebeur, editors. ISIPTA ’15:Proceedings of the Ninth International Symposium on Imprecise Probability: Theories andApplications, 2015.5. Thomas Augustin, GeroWalter, and Frank P. A. Coolen. Statistical inference. In Augustin et al.[2], chapter 7, pages 135–189.6. Jean-Marc Bernard, Teddy Seidenfeld, and Marco Zaffalon, editors. ISIPTA ’03: Proceedingsof the Third International Symposium on Imprecise Probabilities and Their Applications, vol-ume 18 of Proceedings in Informatics, Waterloo, Ontario, Canada, 2003. Carleton Scientific.7. Frank P. A. Coolen, Gert de Cooman, Thomas Fetz, and Michael Oberguggenberger, editors.ISIPTA ’11: Proceedings of the Seventh International Symposium on Imprecise Probability:Theories and Applications, 2011.8. Fabio Cozman, Thierry Denœux, Sébastien Destercke, and Teddy Seidenfeld, editors.ISIPTA ’13: Proceedings of the Eight International Symposium on Imprecise Probability:Theories and Applications, 2013.9. Fabio Gagliardi Cozman, Robert Nau, and Teddy Seidenfeld, editors. ISIPTA ’05: Proceedingsof the Fourth International Symposium on Imprecise Probabilities and Their Applications,2005.10. Gert de Cooman, Fabio Gagliardi Cozman, Serafín Moral, and Peter Walley, editors.ISIPTA ’99: Proceedings of the First International Symposium on Imprecise Probabilitiesand Their Applications, 1999.11. Gert de Cooman, Terrence L. Fine, and Teddy Seidenfeld, editors. ISIPTA ’01: Proceedingsof the Second International Symposium on Imprecise Probabilities and Their Applications,Maastricht, the Netherlands, 2001. Shaker Publishing.50 E. Quaeghebeur12. Gert de Cooman, Jiřina Vejnarová, and Marco Zaffalon, editors. Proceedings of the Fifth Inter-national Symposium on Imprecise Probabilities: Theories and Applications. SIPTA, Action MAgency for SIPTA, 2007.13. Sébastien Destercke and Didier Dubois. Other uncertainty theories based on capacities. InAugustin et al. [2], chapter 5, pages 93–113.14. Sébastien Destercke and Didier Dubois. Special cases. In Augustin et al. [2], chapter 4, pages79–92.15. Scott Ferson, Vladik Kreinovich, Lev Ginzburg, Davis S. Myers, and Kari Sentz. Constructingprobability boxes and dempster-shafer structures. Technical Report SAND2002-4015, SandiaNational Laboratories, 2002.16. Enrique Miranda. A survey of the theory of coherent lower previsions. International Journalof Approximate Reasoning, 48(2):628–658, 2008.17. Enrique Miranda and Gert de Cooman. Lower previsions. In Augustin et al. [2], chapter 2,pages 28–55.18. Erik Quaeghebeur. Desirability. In Augustin et al. [2], chapter 1, pages 1–27.19. Erik Quaeghebeur and Gert de Cooman. Imprecise probability models for inference in expo-nential families. In Cozman et al. [9], pages 287–296.20. Peter Walley. Statistical reasoning with imprecise probabilities, volume 42 ofMonographs onStatistics and Applied Probability. Chapman & Hall, 1991.21. Peter Walley. Inferences from multinomial data: learning about a bag of marbles. Journal ofthe Royal Statistical Society B: Methodological, 58(1):3–57, 1996.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 4Imprecise Discrete-Time Markov ChainsGert de CoomanAbstract I present a short and easy introduction to a number of basic definitionsand important results from the theory of imprecise Markov chains in discrete time,with a finite state space. The approach is intuitive and graphical.4.1 IntroductionAlthough imprecision and robustness in discrete-time Markov chains were alreadystudied in the 1990s [6–8], more significant progress [2, 3, 5, 11] could bemade afterthe graphical structure of imprecise probability trees underlying themwas uncoveredin 2008 [4]. Research has now moved firmly into the continuous-time domain, forwhich [1, 9] are good starting points.In this paper, I give a concise and elementary overview of a number of basicideas and results in discrete-time imprecise Markov chains, with an emphasis ontheir graphical representation. We begin with the basics of precise and impreciseprobability models in Sects. 4.2 and 4.3. When such models are used in a dynami-cal context, precise and imprecise probability trees arise naturally; they and the useof the fundamental Law of Iterated Expectations for making inferences about themconstitute the subjects of Sects. 4.4 and 4.5. Imprecise Markov chains correspondto special imprecise probability trees, and they and their basic inferences are dis-cussed in Sect. 4.6, followed by a number of examples in Sect. 4.7. These exampleshint at stationary distributions and ergodicity. These notions are briefly discussedin Sect. 4.8, which concludes the paper. Throughout, I have included a number ofsimple exercises to illustrate the arguments in the main text.G. de Cooman (B)Foundations Lab for imprecise probabilities, Ghent University, Ghent, Belgiume-mail: gert.decooman@ugent.be© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_45152 G. de Cooman4.2 Precise Probability ModelsAssume we are uncertain about the value that a variable X assumes in some finiteset of possible values X . This is usually modelled by a probability mass function mon X , satisfying (∀x ∈ X )m(x) ≥ 0 and ∑x∈X m(x) = 1.With m we can associate an expectation operator Em as followsEm( f ) :=∑x∈Xm(x) f (x) where f : X → R.If A ⊆ X is an event, then its probability is given by Pm(A) = ∑x∈A m(x) =Em(IA), where IA : X → R is the indicator of A and assumes the value 1 on Aand 0 elsewhere. This tells us that there are two equivalent mathematical languagesfor dealing with uncertainty: the language of probabilities and the language of expec-tations, and that we can go freely from one to the other.All possible (precise) probability models are gathered in the simplex X of allmass functions onX :X :={m ∈ RX : (∀x ∈ X )m(x) ≥ 0 and ∑x∈X m(x) = 1}.Any probability model for uncertainty about X is a point in that simplex, which indi-cated that mass functions have a geometrical interpretation. This is illustrated belowfor the case X = {a, b, c} and the uniform mass function mu.Expectation also has a geometrical interpretation: specifying a value E( f ) for theexpectation of a map f : X → R, namely, ∑x∈X m(x) f (x) = E( f ), imposes alinear constraint on the possible values for m in X . It corresponds to intersectingthe simplex X with a hyperplane, whose direction depends on f . This is alsoillustrated in the picture above; in this particular case two assessments turn out tocompletely determine a unique mass function.4.3 Imprecise Probability ModelsWe now turn to a generalisation of precise probability models, which we will callimprecise. To allow for more realistic and flexible assessments, we can envisageimposing linear inequality—rather than equality—constraints on the m in X :E( f ) ≤∑x∈Xm(x) f (x) or∑x∈Xm(x) f (x) ≤ E( f ).4 Imprecise Discrete-Time Markov Chains 53This corresponds to intersecting X with affine semi-spaces:Any such number of assessments leads to a credal set M, which is our first type ofimprecise probability model.Definition 4.1 A credal set M is a convex closed subset of X .Below, we show some more examples of such credal sets in the special case X ={a, b, c}. The credal set on the left corresponds to the assessment: ‘b is at least aslikely as c’; the one in the middle is a convex mixture of the uniform mass functionwith the entire simplex; and the one on the right represents a statement in classicalpropositional logic: ‘X = a or X = c’. This illustrates that the language of credalsets encompasses both precise probabilities and classical propositional logic.Lower and upper expectations are our second type of imprecise probabilitymodel.To see how they come about, consider the credal set in the figure below on the right.We can ask what we know about the probability of c, or the expectation of I{c},given this credal set: it is only known to belong to the closed interval [1/4, 4/7]. Thiscan be generalised from events to arbitrary elements of the set L(X ) = RX of allreal-valued maps f on X : As m ranges over the credal setM, Em( f ) will similarlyrange over a closed interval that is completely determined by its lower and upperbounds.This leads to the definition of the following two real functionals on L(X ):EM( f ) = min {Em( f ) : m ∈ M} lower expectationEM( f ) = max {Em( f ) : m ∈ M} upper expectationfor all f : X → R.Observe that these lower and upper expectations are mathematically equivalent mod-els, because54 G. de CoomanEM( f ) = −EM(− f ) for all f ∈ L(X ).We will in what follows focus on upper expectations.Exercise 4.1 What is the upper expectation EM when M = X ?This shows thatwe can go from the language of probabilities—and the use ofM—to the language of expectations—and the use of EM. To see that we can also go theother way, we need the following definition:Definition 4.2 Wecall a real functional E onL(X ) anupper expectation if it satisfiesthe following properties: for all f and g in L(X ) and all real λ ≥ 0:1. E( f ) ≤ max f [boundedness];2. E( f + g) ≤ E( f ) + E(g) [sub-additivity];3. E(λ f ) = λE( f ) [non-negative homogeneity].Upper expectations are also called coherent upper previsions [10, 12]. They constitutea model that is mathematically equivalent to credal sets, in very much the same wayas expectations are mathematically equivalent to probability mass functions:Theorem 4.1 A real functional E is an upper expectation if and only if it is theupper envelope of some credal set M.Proof UseM = {m ∈ X : (∀ f ∈ L(X ))(Em( f ) ≤ E( f ))}. Exercise 4.2 Consider any linear prevision Em and any ε ∈ [0, 1]. Verify that theso-called linear-vacuous mixture:is an upper expectation.Solution: Em andmax are upper expectations by Theorem 4.1, because they are upperenvelopes of the respective credal sets {m} and X—see Exercise 4.1. Now verifythat being an upper expectation is preserved by taking convex mixtures. The corre-sponding credal set (1 − ε){m} + εX := {(1 − ε)m + εq : q ∈ X } is indicatedin blue in the figure above. ♦4 Imprecise Discrete-Time Markov Chains 55Exercise 4.3 All upper expectations on a binary space X = {0, 1} are such linear-vacuous mixtures, and the corresponding credal set can be depicted as0 1p pLet p := m(1) and q := 1 − p = m(0). What is the relation between [p, p] andp, ε?Solution: p = E(I{1}) = (1 − ε)p = p − εp and p = E(I{1}) = (1 − ε)p + ε =p + εq. Hence, p − p = ε. ♦4.4 Discrete-Time Uncertain ProcessesWe now apply these ideas in a more dynamic context: the study of processes. Weconsider an uncertain process, which is a collection of uncertain variables X1, X2,…,Xn ,… assuming values in some finite set of statesX . This can be represented graphi-cally by a standard event treewith nodes (also called situations) s = (x1, x2, . . . , xn)for xk ∈ X and n ≥ 0. This is depicted below on the left for the special case thatX = {0, 1}, where we have limited ourselves to three variables X1, X2, and X3; butthe idea should be clear. Observe that we use the symbol for the initial situation,or root node, of the event tree.1(1, 1)(1, 1, 1)(1, 1, 0)(1, 0)(1, 0, 1)(1, 0, 0)0(0, 1)(0, 1, 1)(0, 1, 0)(0, 0)(0, 0, 1)(0, 0, 0)X 1X 2X 31(1, 1)(1, 1, 1)(1, 1, 0)(1, 0)(1, 0, 1)(1, 0, 0)0(0, 1)(0, 1, 1)(0, 1, 0)(0, 0)(0, 0, 1)(0, 0, 0)mm1m0m(1,1)m(0,0)m(0,1)m(1,0)The event tree becomes a probability tree as soon as we attach to each nodes = (x1, x2, . . . , xn) a local probability mass function ms on X with associatedexpectation operator Ems , expressing the uncertainty about the next variable Xn+1after observing the earlier variables X1 = x1, …, Xn = xn . This is depicted aboveon the right for the special case that X = {0, 1}.We now consider a very general inference problem in such a probability tree.Consider any function g : X n → R of the first n variables: g = g(X1, X2, . . . , Xn).We want to calculate its expectation E(g|s) in the situation s = (x1, . . . , xk), that is,after having observed the first k variables. Interestingly, this can be done efficientlyusing the following theorem,which is a reformulation of the Lawof Total Probability:56 G. de CoomanTheorem 4.2 (Law of Iterated Expectations) If we know E(g|s, x) for all x ∈ X ,then we can calculate E(g|s) by backwards recursion using the local model ms:E(g|s) = Ems︸︷︷︸local(E(g|s, ·)) =∑x∈Xms(x)E(g|s, x).This shows that expectations can be calculated recursively using a very basic step,illustrated below for the case X = {0, 1}:s(s, 0)(s, 1)msE(g|s) = ms(1)E(g|s, 1) + ms(0)E(g|s, 0)E(g|s, 1)E(g|s, 0)Hence, all expectations E(g|x1, . . . , xk) in the tree can be calculated from the localmodels ms as follows:1. start in the final cut X n and let E(g|x1, x2, . . . , xn) = g(x1, x2, . . . , xn);2. do backwards recursion using the Law of Iterated Expectations:E(g|x1, . . . , xk) = Em(x1 ,...,xk )︸︷︷︸local(E(g|x1, . . . , xk, ·))3. go on until you get to the root node , where we can identify E(g|) = E(g).Exercise 4.4 Consider flipping a coin twice independently, with probability p forheads—outcome 1—and q = 1 − p for tails—outcome 0. The corresponding prob-ability tree for this experiment is given below on the left, with, in red, in the nodes,the corresponding number of heads. What is the expected number of heads?12p1qp01p0qq12p1qp01p0qqp = p · 1 + q · 01 + p = 2p + q = p · 2 + q · 12p = p + (p · 1 + q · 0)Solution: Above on the right, we apply the Law of Iterated Expectations recursively,from leaves to root; the solution is the expectation 2p attached to the root. ♦Exercise 4.5 Extend the ideas in the solution toExercise 4.4 to calculate the expectednumber of heads when the coin is flipped n times independently.Solution:We apply the Law of Iterated Expectations recursively, from leaves to root.Below on the left, we consider starting from the leaves of the tree at depth n; applyingthe Law reduces to adding p to the number of heads in each of their parent nodes atdepth n − 1. On the right, we apply the Law to these nodes at depth n − 1, whichreduces to adding 2p to the number of heads in each of their parent nodes at depthn − 2.4 Imprecise Discrete-Time Markov Chains 57kk + 1pkqn − 1 nk + p = k + (p · 1 + q · 0)at time n − 1 in a situation with k headskk + p + 1pk + pqn − 2 n − 1k + 2p = k + p + (p · 1 + q · 0)at time n − 2 in a situation with k headsGoing on in this way, we see that the solution is the expectation np attached to theroot at depth 0. ♦Exercise 4.6 We now flip the same coin time and time again, independently, untilwe reach heads for the first time. Calculate the expected number of coin flips.Solution: Below is the (unbounded) probability tree associated with this experiment.0(0, 0)· · ·q(0, 0, 1)pq(0, 1)pq1pdepth 1αα + 1Call the unknown expectation α. We apply the Law of Iterated Expectations tothe situations at depth 1. In the situation 1, the expected number of heads is 1, theactual number of heads there. In the situation 0, we see a copy of the original treeextending to the right, but since we have already flipped the coin once here, theexpected number of heads in this situation is α + 1. In the parent node, the expectednumber of heads α is therefore also given by p · 1 + q · (α + 1) = 1 + qα, whenceα = 1/p. ♦4.5 Imprecise Probability TreesUntil now, we have assumed that we have sufficient information in order to specify,in each node s, a local probability mass function ms on the set X of possible valuesfor the next state.s(s, 0)(s, 1)ms −→ s(s, 0)(s, 1)MsWe now let go of this major restrictive assumption by allowing for more generaluncertaintymodels.Wewill consider credal sets as ourmore general local uncertaintymodels: closed convex subsets Ms of X . See the figure below for a special casewhen X = {0, 1}.Definition 4.3 An imprecise probability tree is an event tree where in each node sthe local uncertainty model is a credal set Ms , or equivalently, its associated upperexpectation Es , with Es( f ) := max {Em( f ) : m ∈ Ms} for all f ∈ L(X ).58 G. de CoomanAn imprecise probability tree can be interpreted as an infinity of compatible preciseprobability trees: choose in each node s a probability mass function ms from thesetMs .1(1, 1)(1, 1, 1)(1, 1, 0)(1, 0)(1, 0, 1)(1, 0, 0)0(0, 1)(0, 1, 1)(0, 1, 0)(0, 0)(0, 0, 1)(0, 0, 0)MM1M0M(1,1)M(0,0)M(0,1)M(1,0)For each real map g = g(X1, . . . , Xn), each node s = (x1, . . . , xk), and each suchcompatible precise probability tree, we can calculate the expectation E(g|x1, . . . , xk)using the backwards recursion method described before. By varying over each com-patible probability tree, we get a closed real interval, completely characterised bylower andupper expectations E(g|x1, . . . xk) and E(g|x1, . . . , xk): [E(g|x1, . . . , xk),E(g|x1, . . . , xk)]. The complexity of calculating these bounds in this way is clearlyexponential in the number of time steps n. But, there is a more efficient method tocalculate them:Theorem 4.3 (Law of Iterated Upper Expectations [4, 5]) If we know E(g|s, x) forall x ∈ X , then we can calculate E(g|s) by backwards recursion using the localmodel Es:E(g|s) = Es︸︷︷︸local(E(g|s, ·)) = maxms∈Ms∑x∈Xms(x) E(g|s, x).This shows that expectations can be calculated recursively using a very basic step,illustrated below for the case X = {0, 1}:s(s, 0)(s, 1)MsE(g|s) = Es(E(g|s, ·))E(g|s, 1)E(g|s, 0)The method for, and the complexity of, calculating the E(g|s), as a function of n, istherefore essentially the same as in the precise case!Exercise 4.7 Extend the ideas in the solution to Exercise 4.5 to calculate the upperexpected number of heads when the coin is flipped n times independently, but wherenowwe have an imprecise probabilitymodel for a coin flip, with a probability interval[p, p] for heads, and a corresponding interval [q, q] = [1 − p, 1 − p] for tails.Solution:We apply the Law of Iterated Upper Expectations recursively, from leavesto root. Below on the left, we consider starting from the leaves of the tree at depth n;4 Imprecise Discrete-Time Markov Chains 59applying the Law reduces to adding p to the number of heads in each of their parentnodes at depth n − 1.kk + 1[p, p]k[q, q]n − 1 nk + p = k + E(I{1})at time n − 1 in a situation with k headskk + p + 1[p, p]k + p[q, q]n − 2 n − 1k + 2p = k + p + E(I{1})at time n − 2 in a situation with k headsOn the right, we apply the Law to these nodes at depth n − 1, which reduces toadding 2p to the number of heads in each of their parent nodes at depth n − 2. Goingon in this way, we see that the solution is the expectation n p attached to the root atdepth 0. A similar result holds for the lower expectation. ♦Exercise 4.8 We now flip the same coin with the imprecise probability model timeand time again, independently, until we reach heads for the first time. Calculate theupper expected number of coin flips.Solution: Below is the (unbounded) probability tree associated with this experiment.Call the unknown upper expectation α. We apply the Law of Iterated Upper Expecta-tions to the situations at depth 1. In the situation 1, the upper expected number of headsis 1, the actual number of heads there. In the situation 0, we see a copy of the originaltree extending to the right, but since we have already flipped the coin once here, theupper expected number of heads in this situation is α + 1. In the parent node, theupper expected number of heads α is therefore also given by 1 + E(α I{0}) = 1 + αq ,whence α = 1/p. A similar result holds for the lower expectation. ♦The attentive reader will have observed that in all these simple exercises, we canalso obtain the ‘imprecise’ result from the ’precise’ one by optimising over thesingle parameter p. We have to warn against too much optimism: in more involvedexamples, this will no longer be the case.4.6 Imprecise Markov ChainsWe now look at a special instance of a probability tree, corresponding to a stationary(precise) Markov chain. This happens when the precise local models m(x1,...,xn) only60 G. de Coomandepend on the last observed state xn—this is theMarkov Condition—and also do notdepend explicitly on the time step n:m(x1,...,xn) = q(·|xn)for some family of transition mass functions q(·|x), x ∈ X .Definition 4.4 The uncertain process is a stationary precise Markov chain when allMs are singletons {ms} and M(x1,...,xn) = {q(·|xn)}, for some family of transitionmass functions q(·|x), x ∈ X .For each x ∈ X , the transition mass function q(·|x) corresponds to an expectationoperator, given by E( f |x) = ∑z∈X q(z|x) f (z) for all f ∈ L(X ).Definition 4.5 Consider the linear transformation T ofL(X ), called transition oper-ator: T : L(X ) → L(X ) : f → T f , where T f is the real map defined by:T f (x) := E( f |x) =∑z∈Xq(z|x) f (z) for all x ∈ X .In the parlance of linear algebra, or functional analysis, T is the dual of the lineartransformation with Markov matrix M with elements Mxy := q(y|x).Up to now, we have mainly been concerned with conditional expectations of thetype E(·|s).Wewill now look at particular unconditional expectations, where s = .For any n ≥ 0, we define the expectation for the (single) state Xn at time n byEn( f ) = E( f (Xn)) = E( f (Xn)|) for all f : X → Rand we denote the corresponding mass function bymn . Applying the Law of IteratedExpectations in Theorem 4.2 now yields, with also E1 = Em and m1 = m:En( f ) = E1(Tn−1 f ), and dually, mn = Mn−1m1,so the complexity of calculating En( f ) is linear in the number of time steps n.Exercise 4.9 Consider the stochastic process where we first flip a fair coin. Fromthen on, on heads, we select a biased coin with probability p for heads for the nextcoin flip, and on tails, a biased coin with probability q = 1 − p for heads, and keepon flipping one of the two biased coins, selected on the basis of the outcome of theprevious coin flip. This produces a Markov chain. Find T f , T2 f , and E1( f ), E2( f )and E3( f ) for f ∈ L({0, 1}).Solution: Clearly, E1( f ) = 1/2 f (1) + 1/2 f (0), T f (0) = E( f |0) = q f (1) + p f (0)and T f (1) = E( f |0) = p f (1) + q f (0), whenceE2( f ) = E1(T f ) = p + q2f (1) + p + q2f (0) = 12f (1) + 12f (0).4 Imprecise Discrete-Time Markov Chains 61Similarly,T2 f (0) = qT f (1) + pT f (0) = q[p f (1) + q f (0)] + p[q f (1) + p f (0)]= (p2 + q2) f (0) + 2pq f (1)T2 f (1) = (p2 + q2) f (1) + 2pq f (0),whenceE3( f ) = E1(T2 f ) = p2 + q2 + 2pq2f (1) + 2pq + p2 + q22f (0) = 12f (1) + 12f (0),and so on.We see that at the level of expectations of single state variables, the processcannot be distinguished from flipping a fair coin. ♦The generalisation from precise to imprecise Markov chains goes as follows:Definition 4.6 The uncertain process is a stationary imprecise Markov chain whenthe Markov Condition is satisfied with stationarity: M(x1,...,xn) = Q(·|xn) for somefamily of transition credal sets Q(·|x), x ∈ X .1(1, 1)(1, 1, 1)(1, 1, 0)(1, 0)(1, 0, 1)(1, 0, 0)0(0, 1)(0, 1, 1)(0, 1, 0)(0, 0)(0, 0, 1)(0, 0, 0)mq(·|1)q(·|0)q(·|1)q(·|0)q(·|1)q(·|0)−→1(1, 1)(1, 1, 1)(1, 1, 0)(1, 0)(1, 0, 1)(1, 0, 0)0(0, 1)(0, 1, 1)(0, 1, 0)(0, 0)(0, 0, 1)(0, 0, 0)MQ(·|1)Q(·|0)Q(·|1)Q(·|0)Q(·|1)Q(·|0)An imprecise Markov chain can be seen as an infinity of (precise) probability trees:choose a precise mass function from Ms in each situation s. It should be clear thatnot all of these satisfy the Markov property or stationarity. This implies that solvingthe optimisation problem in order to find the tight upper bounds E(g|s), as discussedin Sect. 4.5, is not (necessary always) simply an optimisation over a parametrisedcollection of stationary (or even non-stationary) Markov chains, although it can turnout be so simple in a number of special cases.For each x ∈ X , the local transition modelQ(·|x) corresponds to an upper expec-tation operator E(·|x), with E( f |x) = max {Ep( f ) : p ∈ Q(·|x)}for all f ∈ L(X ).This leads to the following definition, which generalises the definition of transitionoperators for precise Markov chains:Definition 4.7 Consider the non-linear transformation T of L(X ), called the uppertransition operator: T : L(X ) → L(X ) : f → T f where the real map T f is definedby T f (x) := E( f |x) = max {Ep( f ) : p ∈ Q(·|x)}for all x ∈ X .62 G. de CoomanFor any n ≥ 0, we define the upper expectation for the (single) state Xn at time n byEn( f ) = E( f (Xn)) = E( f (Xn)|) for all f : X → R.Then the Law of Iterated Upper Expectations of Theorem 4.3 yields, with also E1 =EM :En( f ) = E1(Tn−1 f ) for all n ≥ 1 and all f ∈ L(X ),so the complexity of calculating En( f ) is still linear in the number of time steps n.4.7 ExamplesConsider a two-element state space X = {1, 0}, with upper expectation E1 =EM for the first variable, and for each (x1, . . . , xn) ∈ {1, 0}n , with 0 < ε ≤1, M(x1,...,xn) = Mxn = (1 − ε){q(·|xn)} + ε{1,0}, or equivalently, for the uppertransition operator T = (1 − ε)T + ε max. In other words, each transition credalsetQ(·|x) is a linear-vacuous mixture (see Exercise 4.2, also for the notations used)centred on the transition mass function q(·|x), where the mixture coefficient ε is thesame in each state x .It is amatter of simple and direct verification that for n ≥ 1 and f ∈ L(X ): Tn f =(1 − ε)nTn f + ε ∑n−1k=0(1 − ε)k max Tk f , and therefore, using the Law of IteratedExpectations, En+1( f )=E1(Tn f ) = (1 − ε)n E1(Tn f ) + ε ∑n−1k=0(1 − ε)k max Tk f .If we now let n → ∞, it is not too hard to see that the limit exists and is independentof the initial upper expectation E1:limn→∞ En( f ) = ε∞∑k=0(1 − ε)k max Tk f for all f ∈ L(X ).We consider two special cases:1. Contaminated randomwalk:whenT f (1)=T f (0) = 1/2[ f (1) + f (0)], the under-lying precise Markov chain is actually like flipping a fair coin. We then find thatE∞( f ) = (1 − ε)1/2[ f (1) + f (0)] + ε max f for all f ∈ L(X ).2. Contaminated cycle: when T f (1) = f (0) and T f (0) = f (1), the underlyingprecise Markov chain is actually like deterministic cycle between the states 0and 1. We then find that E∞( f ) = max f for all f ∈ L(X ).The probability intervals for 1 corresponding to these two limit models are given by0 11/2εand0 1ε = 14 Imprecise Discrete-Time Markov Chains 63As another example, we considerX = {a, b, c} and the transitionmodels depictedbelow, which are imprecise models not very far from a simple cycle:Below, we depict the time evolution of the En (as credal sets) for three cases (red,yellow and blue). We see that, here too, regardless of the initial distribution E1, thedistribution En seems to converge to the same distribution.4.8 A Non-linear Perron–Frobenius Theorem,and ErgodicityThe convergence behaviour in the previous examples can also be observed in generalimprecise Markov chains under fairly weak conditions. The following theorems canbe derived from the more general discussions and results in [3, 5].Theorem 4.4 Consider a stationary imprecise Markov chain with finite state set Xand upper transition operator T. Suppose that T is regular, meaning that there issome n > 0 such that min Tn I{x} > 0 for all x ∈ X . Then for every initial upperexpectation E1, the upper expectation En = E1 ◦ Tn−1 for the state at time n con-verges point-wise to the same stationary upper expectation E∞: limn→∞ En(h) =limn→∞ E1(Tn−1h) := E∞(h) for all h in L(X ). The limit upper expectation E∞ isthe only T-invariant upper expectation on L(X ), meaning that E∞ = E∞ ◦ T.In that case we also have an interesting ergodicity result. For a detailed descriptionof the notion of ‘almost surely’, we refer to [3], but it roughly means ‘with upperprobability one’.64 G. de CoomanTheorem 4.5 Consider a stationary imprecise Markov chain with finite state set Xand upper transition operator T. Suppose that T is regular with stationary upperexpectation E∞. Then, almost surely, for all h in L(X ):E∞(h) ≤ lim infn→∞1nn∑k=1h(Xk) ≤ lim supn→∞1nn∑k=1h(Xk) ≤ E∞(h).4.9 ConclusionThe discussion in this paper lays bare a few interesting but quite basic aspects ofinference in imprecise probability trees and Markov chains in discrete time. A moregeneral and deeper treatment of these matters can be found in [3–5]. For recent workon impreciseMarkov chains in continuous time, I refer the interested reader to [1, 9].References1. Jasper De Bock. The Limit Behaviour of Imprecise Continuous-Time Markov Chains. Journalof Nonlinear Science, 27(1):159–196, 2017.2. Jasper De Bock and Gert de Cooman. An efficient algorithm for estimating state sequencesin imprecise hidden Markov models. Journal of Artificial Intelligence Research, 50:189–233,2014.3. Gert de Cooman, Jasper De Bock, and Stavros Lopatatzidis. Imprecise stochastic processes indiscrete time: global models, imprecise markov chains, and ergodic theorems. InternationalJournal Of Approximate Reasoning, 76(C):18–46, 2016.4. Gert de Cooman and Filip Hermans. Imprecise probability trees: Bridging two theories ofimprecise probability. Artificial Intelligence, 172(11):1400–1427, 2008.5. Gert de Cooman, Filip Hermans, and Erik Quaeghebeur. Imprecise Markov chains and theirlimit behaviour. Probability in the Engineering and Informational Sciences, 23(4):597–635,2009. arXiv:0801.0980.6. D. J. Hartfiel and E. Seneta. On the theory of Markov set-chains. Advances in Applied Proba-bility, 26:947–964, 1994.7. Darald J. Hartfiel. Sequential limits in Markov set-chains. Journal of Applied Probability,28(4):910–913, 1991.8. Darald J. Hartfiel. Markov Set-Chains. Number 1695 in Lecture Notes in Mathematics.Springer, Berlin, 1998.9. Thomas Krak, Jasper De Bock, and Arno Siebes. Imprecise continuous-time Markov chains.International Journal of Approximate Reasoning. 88:452–528, 2017.10. Matthias C. M. Troffaes and Gert de Cooman. Lower Previsions. Wiley, 2014.11. Damjan Škulj and Robert Hable. Coefficients of ergodicity for Markov chains with uncertainparameters. Metrika, 76(1):107–133, 2013.12. Peter Walley. Statistical Reasoning with Imprecise Probabilities. Chapman and Hall, London,1991.4 Imprecise Discrete-Time Markov Chains 65Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 5Statistics with ImpreciseProbabilities—A Short SurveyThomas AugustinAbstract This chapter aims at surveying and highlighting in an introductory waysome challenges and big opportunities a paradigmatic shift to imprecise probabili-ties could induce in statistical modelling. Working with an informal understandingof imprecise probabilities, we discuss the concepts of model imprecision and dataimprecision as the two main types of imprecision in statistical modelling. Then weprovide a short survey of some major developments, methodological questions andapplications of imprecise probabilistic models under model imprecision in the con-text of different inference schools and summarize some recent developments in thearea of data imprecision.5.1 IntroductionBy promising powerful solutions to some of the deepest foundational problems ofprobability and statistics, imprecise probabilities offer great opportunities also forthe area of statistical modelling. Still, in statistics, theories of imprecise probabili-ties live in the shadows, and admittedly the development of many of the impreciseprobability-based methods is often in a comparatively early stage. Nevertheless, inall areas of statistics, the desire for a more comprehensive modelling of complexuncertainty had popped up again and again. The rather scattered work covers a hugevariety of methods and topics, ranging from contributions to the methodological andphilosophical foundations of inference to very concrete questions of applications.T. Augustin (B)Department of Statistics, Ludwig-Maximilians-Universität München (LMU Munich),Munich, Germanye-mail: augustin@stat.uni-muenchen.de© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_56768 T. AugustinThe present chapter aims at providing a rough and informal survey of some of themajor questions and developments.1 It is structured as follows. Section5.2 collectsthe basic concepts that are needed later on. Then, Sect. 5.3 looks at the major sourcesof imprecision in statistical modelling. We distinguish there between several typesof data imprecision and of model imprecision. Section5.4 focuses on the issue ofmodel imprecision and discusses it from the angle of different inference schools. Weput some emphasis on the (generalized) frequentist and Bayesian setting, but alsobriefly adopt other perspectives. In Sect. 5.5, approaches to handle so-called onticand epistemic data imprecision, respectively, are surveyed. Section5.6 is reservedfor some concluding remarks.5.2 Some Elementary Background on ImpreciseProbabilitiesIn this section, we briefly summarize the concepts in the background. With respectto the basic notions and the technical framework for statistical inference, we refer tothe first chapter of this book [41]. We rely on the same basic setting, where we useobservations (data) on some underlying stochastic phenomenon2 to learn character-istics of that mechanism, mathematically described by an unknown parameter of theunderlaid probability model.With respect to imprecise probabilities, a very rough and eclectic understandingshall be sufficient to read this chapter.3 It is not necessary, for our aims here, to dis-tinguish different approaches with respect to many technical details. Thus, in a ratherinclusive manner, we subsume here under imprecise probabilities any approach thatreplaces in its modelling precise, traditional probabilities p(·) by non-empty sets Pof precise probabilities as the basic modelling entity, including also all approachesthat can be equivalently transferred into a set of precise probability. This comprisesapproaches directlyworkingwith sets of probabilities, like robust Bayes analysis (seeSect. 5.4.2), Kofler &Menges’ linear partial information (e.g. [42]), Levi’s approachto epistemology (e.g. [45]), as well as the whole bunch of corresponding approachesbased on non-linear functionals and non-additive set-functions, covering lower andupper previsions in tradition ofWalley’s book [66], interval probabilities building on1A longer survey discussing the state of the art of statistical inference with imprecise probabilitiesat that time is aimed at by [8].2 To avoid any commitment to a certain interpretation of probability, we use the term “stochasticphenomenon” as a kind of neutral, superordinate concept, in particular including approaches refer-ring to random phenomena only as well as approaches addressing generally situations of epistemicuncertainty.3An introduction into imprecise probabilities on an intermediate level is aimed at by [4]; see also[57] (in this volume) for an introduction, [15] for a survey by a philosopher and [9] for a reviewwithan engineering background. Current developments in research on imprecise probabilities are mostlydiscussed at the biannual ISIPTA (International Symposium on Imprecise Probabilities: Theoriesand Applications) meetings; see [1, 5, 27] for the most recent ones.5 Statistics with Imprecise Probabilities—A Short Survey 69Weichselberger ([72]4), probabilistically interpretable approaches based on capaci-ties including the suitable branch of Dempster-Shafer theory (e.g. [29]), random sets(e.g. [12, Chap. 3]) and p-boxes following [32].Moreover, there is a smooth transitionto several approaches propagating systematic sensitivity analysis (e.g. [52]).It is important that our basic entity, the set P, has to be understood and treated asan entity of its own. It is by no means possible to distinguish certain of its elementsas more likely than others or to mix its elements, eventually leading to a precisetraditional probability distribution. P may be assessed directly by collecting severalprecisemodels, so-to-say as possibleworlds / expert opinions, to be considered.Moreoften P is constructed, typically as the set of all probability distributions• that respect bounds on the probabilities of certain classes of events or, more gen-erally, the expectations of certain random variables,5, 6• or that are (in a topologically well-defined sense) close to a certain precise prob-ability distribution p0(·) (neighbourhood models), providing a formal frameworkfor expressing statements like “p0(·) is approximately true”,7• or that are described by a parameterϑ varying in an interval/cuboid (parametricallyconstructed models), like “imprecise versions” of a normal distribution.85.3 Types of Imprecision in Statistical ModellingThere are many situations in statistics where imprecision occurs, i.e. where a carefulmodelling should go beyond the idealized situation of perfect stochasticity and dataobservedwithout any error and in an ideal precision. To study these situations further,an ideal-typical distinction between model and data precision is helpful.Model Imprecision has to be taken into account whenever there is doubt in a con-crete application that the strong requirements (perfect precision and, by the additivityaxiom, absolute internal consistency) the traditional concept of probability calls forcan be realistically satisfied. This includes all situations of incomplete or conflictinginformation, robustness concerns, and repeated sampling from a population wherethe common assumption of i.i.d. repetitions is violated by some unobserved hetero-geneity or hidden dependence. In a Bayesian setting, in addition, quite often, if atall, a precise prior distribution is not honestly deducible from the knowledge actuallyavailable. Ellsberg’s [31] seminal thought experiments on urns with only partiallyknown compositions have made it crystal clear that the amount of ambiguity, i.e. the4 See also [73] for a summary in English language.5 For instance, in the approaches founded by [66, 72], respectively.6 It may be explicitly noted that the latter also includes comparative/ordinal probabilities (e.g. [49]),only ordering the probability of certain events like “A is more likely to occur than B, and B morelikely than C’.7 See Sect. 5.4.1.8 See, for instance, the sets of conjugated distributions in Sect. 5.4.2 for an example.70 T. Augustinuncertainty about the underlying stochastic mechanism, plays a constitutive role indecision making and thus has to modelled carefully.Data Imprecision comprises all situations where the observations are not avail-able in the granularity that is originally intended in the corresponding application.Following [26], for the modelling of data imprecision, it is crucial to distinguishtwo situations, the precise observation of something inherently imprecise (ontic dataimprecision) and the imprecise observation of something precise (epistemic dataimprecision).9 The difference between these two may be explained by a pre-electionstudy where some voters are still undecided which single party they will vote for.(Compare [55], and for more details [44, 53]). If we understand a voting preferencelike “I am still undecided between parties A and B” as a political position of its own,then we interpret the information as an instance of ontic imprecision. If we focuson the forthcoming election and take this information as an imprecise observationallowing us to predict that the voting decision on the election day will either beA or B, then we are coping with epistemic imprecision. In particular in engineer-ing, another frequent example of epistemic data imprecision occurs from insufficientmeasurement precision, where intervals instead of precise values are observed. Acontrasting example where (unions of) intervals are to be understood as an entity oftheir own arises when the time span of certain spells is characterized: “This machinewas under full load from November 10th to December 23rd.”Epistemic imprecision very naturally occurs in many studies on dynamic pro-cesses, from socio-economic over medical to technical studies. There, censoring isalways a big issue: the spells of some units are still unfinished when the study ends,providing only lower bounds on the spell duration. Typical examples include theduration of unemployment, the time to recurrence of a tumour or the lifetime of anelectronic component. In addition, also interval censoring is quite common, whereone only learns that the event of interest occurred within a certain time span. It shouldbe noted explicitly that missing data, a frequent problem for instance in almost everysocial survey, may be comprised under this setting, taking the whole sample spaceas the observation for those units missing.5.4 Statistical Modelling Under Model ImprecisionOf course, there are also strong reservations against imprecise probabilities in tra-ditional statistics. For a traditional statistician, imprecise probabilities are just asuperfluous complication, misunderstanding either the generality of the concept ofuncertainty or the reductionist essence of the modelling/abstraction process. Indeed,for an orthodoxBayesian, all kinds of not-knowing are simply situations under uncer-tainty, and any kind of uncertainty is eo ipso expressible by traditional probabilities.From the modelling perspective, imprecision is taken as part of the residual category9 See also the distinction between the disjunctive and conjunctive interpretation of random sets, asdiscussed, e.g. in [30, Sect. 1.4].5 Statistics with Imprecise Probabilities—A Short Survey 71that is naturally lost when abstracting and building models.10 Box (and Draper)’soften cited dictum “Essentially, all models are wrong, but some of them are useful”[14, p. 424] has generally been (mis)understood as a justification to base details ofthe model choice on mathematical convenience and as an irrevocable argument totake model imprecision as negligible.5.4.1 Probabilistic Assumptions on the Sampling ModelMatter: Frequentist Statistics and ImpreciseProbabilitiesIndeed, however, Box’s quotation could be continued by “…and some models aredangerous”, since the general neglecting of model imprecision implicitly presup-poses a continuity (in an informal sense) of the conclusions in the models that byno means can be taken as granted. A well-known example from robust statistics11is statistical inference from a “regular bell-shaped distribution”. The standard pro-ceeding would be to assume a normal distribution. But, for instance, the density ofa Cauchy distribution is phenomenologically de facto almost indistinguishable fromthe density function of a normal distribution, suggesting bothmodels to be equivalentfrom a practical point of view. However, statistical inference based on the samplemean shows fundamentally different behaviours.12 Under normality, the distributionof the sample mean behaves nicely, contracting itself around the correct value if thesample size n increases. In the case of the Cauchy distribution, however, the distri-bution of the sample mean stays the same irrespective of the sample size,13 makingany learning by the sample mean impossible.This shocking insight—optimal statistical procedures may behave disastrouslyeven under “tiny deviations” from the ideal model—demonstrates that imprecisionin the underlying model may matter substantially. In this context, the theory of (fre-quentist) robust statistics as the theory of approximately true models emerged, andimprecise probabilities provide a natural superstructure upon it (see, e.g. [38] forhistorical connections). In particular, neighbourhood models, also briefly mentionedabove, have become attractive.14 Building on an influential result by Huber andStrassen [39], a comprehensive theory of testing in situations where the hypotheses10 This is particular the case when the prescriptive character of models is emphasized over itsdescriptive role (cf., e.g. [13, p. 99]).11 See also [8, Sect. 7.5.1, in particular, Fig. 7.3] to illustrate this.12 Generally, if one looks at an i.i.d. sample X1, . . . , Xn from a normal distribution with parametersμ and σ 2, denoted by N(μ, σ 2), and an i.i.d. sample Z1, . . . , Zn from a Cauchy distribution withparameter α and β, denoted by C(α, β), respectively, one obtains for the sample means X̄ :=1n∑ni=1 Xi ∼ N(μ, σ2n ) and Z̄ := 1n∑ni=1 Zi ∼ C(α, β).13 The Law of Large Numbers is not applicable to the Cauchy distribution because it does not havemoments.14 See, e.g. [6] for a survey.72 T. Augustinare described by imprecise probabilities emerged (see [2], [8, Sect. 7.5.2] and thereferences in the corresponding review sections therein.) Further insights are pro-vided from interpreting frequentist statistics as a decision problemunder an imprecisesampling distribution, leading to the framework investigated in [34].Other frequentist approaches, starting from different angles, include Hampel’sfrequentist betting approach (e.g. [36]) and some work on minimum distance esti-mation under imprecise sampling models (e.g. [35]). The statistical consequences ofthe chaotic models in the genuine frequentist framework to imprecise probabilities,developed by Fine and followers (e.g. [33]), might be quite intriguing, but are stillalmost entirely unexplored.5.4.2 Model Imprecision and Generalized Bayesian InferencePriors and Sets of Priors, Generalized Bayes Rule. The centrepiece of Bayesianinference is the prior distribution. Apart from very large sample sizes, where theposterior is de facto determined by the sample, the prior naturally has a stronginfluence on the posterior and on all conclusions drawn from it. In the rare situationswhere very strong prior knowledge is available, it can be used actively, but mostoften the strong dependence on the prior has been intensively debated and criticized.Working with sets of prior probabilities (or interval-valued priors) opens newavenues here. This set can naturally be chosen to reflect the quality/determinacy ofprior knowledge: strong prior knowledge leads to “small” sets; weak prior knowl-edge to “large” sets. Typical model classes include neighbourhood models or sets ofparametric distributions which often are conjugate 15 to the sampling distribution,which typically still is assumed to be precise.16 In imprecise probability, is under-stood as naturally inducing the set x of all posteriors arising from a prior in .17Interpretations of andx vary to the extent they are understood as principled enti-ties. A pragmatic point of view sees an investigation of x just as a self-evident wayto perform a sensitivity analysis. On the other extreme, Walley’s [66] generalizedtheory of coherence, having initiated the most vivid branch of research on impreciseprobabilities, provides a rigorous justification of exactly this way to proceed as the“Generalized Bayes Rule (GBR)”. Important developments have also been achievedfor a variety of different model classes under the term “Robust Bayesian Analysis”;see, e.g. [58] for a review on this topic.Near Ignorance Models. One importantway to use sets of priors is that they allowfor quite a natural formulation of (rather) complete ignorance. A traditional Bayesian15 See, e.g. [41, Chap. 1.2] (in this volume).16 Imprecise sampling models are quite rarely studied, see, however, [66, Sect. 8.5] and [61].17 There are some approaches in the traditional Bayesian context that work with sets of priors, too,but merely understanding such a set as an interim step, based on which an ultimate precise prior isselected. This can be done in a data-driven way, as in empirical Bayes approaches, including theML-II approach (e.g. [11, Sect. 3.5.4]), or generally, as in maximum entropy approaches (e.g. [11,Sect. 3.4]).5 Statistics with Imprecise Probabilities—A Short Survey 73model eo ipso fails in expressing ignorance/non-informativeness. Assigning a pre-cise probability is never non-committal; every precise prior delivers probabilisticinformation about the parameter. The genuinely non-informative model is the setof all probability distributions. While this model would yield vacuous inferences,it motivates so-called near-ignorance models, where, informally spoken, the innercore of this set is used, excluding extreme probabilities that are immune to learning.Near-ignorance models still assign non-committal probabilities to standard eventsin the parameter space, but allow for learning. By far the most popular model is theImprecise Dirichlet Model (IDM) [67] for categorical data. Different extensions fol-lowed, including general near-ignorance models for exponential families (e.g. [10]).Another direction of enabling the formulation of near-ignorance uses all priors withbounded derivatives ([68]; for a general exposition of the concept of bounded influ-ence, see the book [69]).Prior-Data Conflict. In some sense, a complementary application of generalizedBayesian inference is the active modelling of prior-data conflict. In practice, general-ized Bayesian models are quite powerful in expressing substantial prior knowledge.In particular, in areas where data are scarce, it is important to use explicitly all priorknowledge available, for instance by borrowing strength from similar experiments.Then, however, it is crucial to have some kind of alert system warning the analystwhen the prior knowledge appears doubtful in the light of data. Indeed, sets of pri-ors can be designed to react naturally to potential prior-data conflict: If data andprior assumptions are in agreement, the set of posterior distributions contracts moreand more with increasing sample size. In contrast, in the case of prior-data conflictand intermediate sample size, the set of posterior distributions is inflated substan-tially, perfectly indicating that one should refrain from most decisions before havinggathered further information.185.4.3 Some Other ApproachesWith generalized Bayesian approaches, and less pronounced with generalized fre-quentist statistics, the major statistical inference school are also predominant in thearea of imprecise probabilities. Nevertheless, there has also been considerable suc-cess in other inference frameworks. Again and again, the desire to save Fisher’sfiducial argument, aiming at providing probability statements on parameters withouthaving to rely on a prior distribution, has been a driving force for developmentsin imprecise probabilities. Dempster’s concept of multivalued mappings (e.g. [28]),which become evenmore famous in artificial intelligence by Shafer’s reinterpretationfounding Dempster-Shafer Theory (see, for instance, again the survey by [29]), is tobe mentioned here, but also work by Hampel (e.g. [36]) and by Weichselberger (e.g.18 See [71] for a general treatment in one-parameter exponential families, also providing someillustrating plots, and, for instance, [70] for an application in the context of system reliability.74 T. Augustin[74]), see also [7] tracing back its roots.19 A generalized likelihood-based frameworkhas been introduced by Cattaneo (e.g. [16, 17]).Another direct inference approach is Nonparametric Predictive Inference (NPI),as introduced by Coolen originally for Bernoulli data [18]. Based on exchangeabilityarguments, NPI yields direct conditional probabilities of further real-valued randomquantities, relying on the low structure assumption that all elements of the naturalpartition produced by the already observed data are equally likely; see, for instance,[19] for a detailed discussion, [3] for a clear embedding into imprecise probabilitiesand [21] for a web-page documenting research within this framework. The basicapproach can be naturally extended to censored data / competing risks (e.g. [22]),and to categorical data [20]. NPI has been developed further in a huge variety offields; see, for instance, [23, 24] for recent applications in biometrics and finance,respectively.5.5 Statistical Modelling Under Data ImprecisionIn this section, we turn to statistical modelling under data imprecision. Keeping thedistinction from Sect. 5.3, we briefly discuss ontic data imprecision and then turn toepistemic data imprecision.Ontic data imprecision, where we understand the imprecise observation as anentity of its own, may be argued to be a border case between classical statisticsand its extensions. Technically, we change the sample space of each observation to(an appropriate subset of) the power set. For instance, recalling the election examplefromSect. 5.3, thismean that instead of {a, b, c, . . .} representing the vote for a singleparty, we now also allow for combinations {a, b}, {b, c}, . . ., {a, b, c}, . . ., represent-ing the indecision between several parties. As long as we are in amultinomial setting,nothing has changed from an abstract point of view, providing powerful opportunitiesfor complex statistical modelling. In the spirit of this idea, [44, 55] applymultinomialregression models, classification trees, regularized discrete choice models from elec-tion research, and spectral clustering methods to German pre-election survey data.The situation changes substantially when ordinal or continuous data are considered,because, after changing to the power set, the underlying ordering structure is onlypartially preserved.Epistemic data imprecision is, as the examples at the end of Sect. 5.3 show, ofgreat importance in many applications and is quite vividly addressed in classicalstatistics. Even here, traditional statistics keeps its focus on full identification, i.e.the selection of one single probability model fitting the observed data optimally. Onesearches for, and then implicitly relies on, conditions under which one gets handson the so-to-say deficiency process as a thought pattern, making ideal precise obser-vations imprecise. For that purpose, most classical approaches assume either some19It may also be argued that, although less explicitly aiming at fiducial reasoning, the work onnear-ignorance models discussed above fits well into this category.5 Statistics with Imprecise Probabilities—A Short Survey 75kind of uninformativeness of the deficiency process (independent censorship, coars-ening at random (CAR) or missingness (completely) at random (MCAR, MAR)) oran explicit modelling of the deficiency process; see the classical work by [37, 46].Both the uninformativeness as well as the existence of a precisely specifiable defi-ciency process are very strong assumptions. They are—eo ipso by making explicitstatements about unobservable processes—typically not empirically testable. When-ever these assumptions are just made for purely formal reasons, the price to pay forthe seemingly precise result of the estimation process is high. In terms of Manski’sLaw of Decreasing Credibility,20 results may suffer severely from a loss of theircredibility, and thus of their practical relevance.Against this background, in almost any area of application, the desire for lesscommittal, cautious handling of epistemic data imprecision arose. Mostly isolatedapproacheswere proposed that explicitly try to take all possibleworlds into account ina reliable way, aiming at the set of all models optimally compatible with potentiallytrue data. These approaches include, for instance, work from reliable computingand interval analysis in engineering, like [51], extensions of generalized Bayesianinference (e.g. [75]) to reliable statistics in social sciences (e.g. [56]); see also [8,Sect. 7.8.2], who try to characterize and unify these approaches by the concept ofcautious data completion, and the concept of collection regions in [60].There is a smooth transition to approaches that explicitly introduce cautious mod-elling into the construction of estimation procedures; see, for instance, for recentdifferent likelihood- and loss minimization-based approaches addressing epistemicdata imprecision, [25, 40, 43, 54]. Such approaches have the important advantage thattheir construction often also allows the incorporation of additional well-supportedsubject matter knowledge, too imprecise to be useful for the precision focused meth-ods from traditional statistics, but very valuable to reduce the set of compatiblemodels by a considerable extent.Congenial is work in the field of partial identification and systematic sensitivityanalysis, providing methodology for handling observationally equivalent models;see [48, 65], respectively, for classical work and [47, 62] for introductory surveys.The framework of partial identification is currently receiving considerable attentionin econometrics, where in particular the embedding of fundamental questions intothe framework of random sets is of particular importance [50].5.6 Concluding RemarksThe contribution provided a—necessarily painfully selective—survey of some devel-opments of statistical modelling with imprecise probabilities (in a wider sense, alsoincluding closely related concepts). Both in the area of model imprecision as well asunder data imprecision, imprecise probabilities prove to be powerful and particularly20 “The credibility of inferences decreases with the strength of the assumptions maintained.” [48,p. 1].76 T. Augustinpromising. Further developments urgently needed include a proper methodology forsimulations with imprecise probabilities (see [64] for recent results), a careful studyof the statistical consequences of the rather far developed probabilistic side of thetheory of stochastic processes with imprecise probabilities (e.g. [63]), a more fruit-ful exchange with recent research on uncertainty quantification in engineering (see,e.g. [59] (in this volume)), an open mind towards recent developments in machinelearning and more large scale applications. Not only for these topics it is important tocomplement the still recognizable focus on so-to-say defensive modelling by a moreactive modelling. Far beyond sensitivity and robustness aspects, imprecision canactively be used as a strong modelling tool. The proper handling of prior-data con-flict and the successful incorporation of substantive matter knowledge in statisticalanalysis under data imprecision are powerful examples of going in this direction.References1. A. Antonucci, G. Corani, I. Couso, and S. Destercke, editors. ISIPTA ’17, volume 62 of Pro-ceedings of Machine Learning Research. PMLR, 2017.2. T. Augustin. Neyman-Pearson testing under interval probability by globally least favorablepairs: ReviewingHuber-Strassen theory and extending it to general interval probability. Journalof Statistical Planning and Inference, 105:149–173, 2002.3. T.Augustin and F. Coolen.Nonparametric predictive inference and interval probability. Journalof Statistical Planning and Inference, 124:251–272, 2004.4. T. Augustin, F. P. A. Coolen, G. de Cooman, andM. Troffaes, editors. Introduction to ImpreciseProbabilities. Wiley, Chichester, 2014.5. T. Augustin, S. Doria, E. Miranda, and E. Quaeghebeur, editors. ISIPTA ’15. SIPTA, 2015.6. T. Augustin and R. Hable. On the impact of robust statistics on imprecise probability models:A review. Structural Safety, 32:358–365, 2010.7. T. Augustin and Rudolf Seising. Kurt Weichselberger’s contribution to imprecise probabilities.International Journal of Approximate Reasoning, 98:132–145, 2018.8. T. Augustin, G. Walter, and F. Coolen. Statistical inference. In T. Augustin, F. Coolen,G. de Cooman, and M. Troffaes, editors, Introduction to Imprecise Probabilities, pages 135–189. Wiley, 2014.9. M. Beer, S. Ferson, and V. Kreinovich. Imprecise probabilities in engineering analyses.Mechanical Systems and Signal Processing, 37:4–29, 2013.10. A. Benavoli andM. Zaffalon. Prior near ignorance for inferences in the k-parameter exponentialfamily. Statistics, 49:1104–1140, 2014.11. J. O. Berger. Statistical Decision Theory and Bayesian Analysis. Springer, 2nd edition, 1985.12. A. Bernardini and F. Tonon. Bounding Uncertainty in Civil Engineering – Theoretical Back-ground. Springer, Berlin, 2010.13. J. M. Bernardo and A. Smith. Bayesian Theory. Wiley, Chichester, 2000.14. G. Box and N. Draper. Empirical Model-building and Response Surface. Wiley, New York,1987.15. S. Bradley. Imprecise probabilities. In Edward N. Zalta, editor, The Stanford Encyclopediaof Philosophy (Spring 2019 Edition). Standford University, 2019. https://plato.stanford.edu/archives/spr2019/entries/imprecise-probabilities/.16. M. Cattaneo. Statistical Decisions Based Directly on the Likelihood Function. PhD thesis, ETHZurich, 2007.17. M. Cattaneo. Likelihood decision functions. Electronic Journal of Statistics, 7:2924–2946,2013.5 Statistics with Imprecise Probabilities—A Short Survey 7718. F. Coolen. Low structure imprecise predictive inference for Bayes’ problem. Statistics & Prob-ability Letters, 36:349–357, 1998.19. F. Coolen. On nonparametric predictive inference and objective Bayesianism. Journal of Logic,Language and Information, 15:21–47, 2006.20. F. Coolen and T. Augustin. A nonparametric predictive alternative to the imprecise Dirichletmodel: The case of a known number of categories. International Journal of ApproximateReasoning, 50:217–230, 2009.21. F. Coolen, P. Coolen-Schrijner, and T. Coolen-Maturi. Nonparametric Predictive Inference(npi), 2020. https://npi-statistics.com/.22. T. Coolen-Maturi. Nonparametric predictive pairwise comparison with competing risks. Reli-ability Engineering & System Safety, 132:146–153, 2014.23. T. Coolen-Maturi. Predictive inference for best linear combination of biomarkers subject tolimits of detection. Statistics in Medicine, 36:2844–2874, 2017.24. T. Coolen-Maturi and F. Coolen. Nonparametric predictive inference for the validation of creditrating systems. Journal of the Royal Statistical Society: Series A, 182:1189–1204, 2019.25. I. Couso and D. Dubois. A general framework for maximizing likelihood under incompletedata. International Journal of Approximate Reasoning, 93:238–260, 2018.26. I. Couso, D. Dubois, and L. Sánchez. Random Sets and Random Fuzzy Sets as Ill-PerceivedRandom Variables. Springer, 2014.27. J. De Bock, C.. de Campos, E. de Cooman, G.and Quaeghebeur, and G. Wheeler, editors.ISIPTA ’19, volume 103 of Proceedings of Machine Learning Research. PMLR, 2019.28. A. Dempster. Upper and lower probabilities induced by a multivalued mapping. The Annals ofMathematical Statistics, 38:325–339, 1967.29. T. Denoeux. 40 years of Dempster-Shafer theory. International Journal of Approximate Rea-soning, 79:1–6, 2016.30. D. Dubois. Belief structure, possibility theory and decomposable confidence measures on finitesets. Computers and Artificial Intelligence, 5:403–416, 1986.31. D. Ellsberg. Risk, ambiguity, and the Savage axioms.Quarterly Journal of Economics, 75:643–669, 1961.32. S. Ferson, L. Ginzburg, V. Kreinovich, D. Myers, and K. Sentz. Constructing probability boxesand Dempster-Shafer structures (Sandia report), 2003. https://doi.org/10.2172/809606.33. P. Fierens. An extension of chaotic probability models to real-valued variables. InternationalJournal of Approximate Reasoning, 50:627–641, 2009.34. R. Hable.Data-based Decisions under Complex Uncertainty. PhD thesis, LMUMunich , 2009.35. R. Hable. Minimum distance estimation in imprecise probability models. Journal of StatisticalPlanning and Inference, 140:461–479, 2010.36. F. Hampel. The proper fiducial argument. In R. Ahlswede, L. Bäumer, N. Cai, H. Aydinian,V. Blinovsky, C. Deppe, and H. Mashurian, editors, General Theory of Information Transferand Combinatorics, volume 4123, pages 512–526. Springer, 2006.37. D. Heitjan and D. Rubin. Ignorability and coarse data. The Annals of Statistics, 19:2244–2253,1991.38. P. Huber. The use of Choquet capacities in statistics. Proceedings of the 39th Session of theInternational Statistical Institute, 45:181–191, 1973.39. P. Huber and V. Strassen. Minimax tests and the Neyman-Pearson lemma for capacities. TheAnnals of Statistics, 1:251–263, 1973.40. E. Hüllermeier, S. Destercke, and I. Couso. Learning from imprecise data: Adjustments ofoptimistic and pessimistic variants. In N. Ben Amor, B. Quost, and M. Theobald, editors,Scalable Uncertainty Management, pages 266–279. Springer, 2019.41. G. Karagiannis. Introduction to Bayesian statistical inference. In L. Aslett, F. Coolen, and J. DeBock, editors, Uncertainty in Engineering – Introduction to Methods and Applications, pages1–14. Springer, 2022.42. E. Kofler and G.Menges. Entscheiden bei unvollständiger Information. Springer, Berlin, 1976.43. D. Kreiss and T. Augustin. Undecided voters as set-valued information – towards forecastsunder epistemic imprecision. In K. Tabia and J. Davis, editors, Scalable Uncertainty Manage-ment. Springer LNCS, 242–250, 2020.78 T. Augustin44. D. Kreiss, M. Nalenz, and T. Augustin. Undecided voters as set-valued information - machinelearning approaches under complex uncertainty. In Destercke S. and E. Hüllermeier, editors,ECML/PKDD 2020 Workshop on Uncertainty in Machine Learning, 2020. https://sites.google.com/view/wuml-2020/program.45. I. Levi.The Enterprise of Knowledge. An Essay on Knowledge, Credal Probability, and Chance.MIT Press, Cambridge, 1980.46. R. Little and D. Rubin. Statistical Analysis with Missing Data. Wiley, 2nd edition, 2014.47. C. Manski. Credible interval estimates for official statistics with survey nonresponse. Journalof Econometrics, 191:293–301, 2015.48. Charles Manski. Partial identification of probability distributions. Springer, 2003.49. E. Miranda and S. Destercke. Extreme points of the credal sets generated by comparativeprobabilities. Journal of Mathematical Psychology, 64:44–57, 2015.50. I. Molchanov and F. Molinari. Random Sets in Econometrics. Econometric Society Mono-graphs. Cambridge University Press, 2018.51. H.Nguyen, V.Kreinovich, B.Wu, andG.Xiang.Computing Statistics under Interval and FuzzyUncertainty: Applications to Computer Science and Engineering. Springer, Berlin, 2011.52. M. Oberguggenberger, J. King, and B. Schmelzer. Classical and imprecise probability methodsfor sensitivity analysis in engineering: A case study. International Journal of ApproximateReasoning, 50:680–693, 2009.53. J. Plass. Statistical modelling of categorical data under ontic and epistemic imprecision. PhDthesis, LMU Munich, 2018.54. J. Plass, M. Cattaneo, T. Augustin, G. Schollmeyer, and C. Heumann. Reliable inference in cat-egorical regression analysis for non-randomly coarsened observations. International StatisticalReview, 87:580–603, 2019.55. J. Plass, P. Fink, N. Schöning, and T. Augustin. Statisticalmodelling in surveyswithout neglect-ing ‘The undecided’. InT.Augustin, S.Doria,M.Miranda, andE.Quaeghebeur, editors, ISIPTA’15, pages 257–266. SIPTA, 2015.56. U. Pötter. Statistical Models of Incomplete Data and Their Use in Social Sciences. HabilitationThesis, Ruhr-Universität Bochum, 2008.57. E. Quaeghebeur. Introduction to the theory of imprecise probability. In L. Aslett, F. Coolen, andJ. De Bock, editors, Uncertainty in Engineering – Introduction to Methods and Applications,pages 37–50. Springer, 2022.58. F. Ruggeri, D. Ríos Insua, and J. Martín. Robust Bayesian analysis. In D. Dey and C. Rao,editors, Handbook of Statistics. Bayesian Thinking: Modeling and Computation, volume 25,pages 623–667. Elsevier, 2005.59. S.Bi andM. Beer. Overview of stochasticmodel updating in aerospace application under uncer-tainty treatment. In L. Aslett, F. Coolen, and J. De Bock, editors, Uncertainty in Engineering– Introduction to Methods and Applications, pages 115–130. Springer, 2022.60. G. Schollmeyer and T. Augustin. Statistical modeling under partial identification: Distinguish-ing three types of identification regions in regression analysis with interval data. InternationalJournal of Approximate Reasoning, 56, Part B:224–248, 2015.61. N. Shyamalkumar. Likelihood robustness. In D. Rios Insua and F. Ruggeri, editors, RobustBayesian Analysis, pages 127–143. Springer, New York, 2000.62. E. Tamer. Partial identification in econometrics. Annual Review of Economics, 2:167–195,2010.63. N. T’Joens, J. De Bock, and G. de Cooman. In search of a global belief model for discrete-time uncertain processes. In J. De Bock, C. de Campos, G. de Cooman, E. Quaeghebeur, andG. Wheeler, editors, ISIPTA ’19, volume 103 of Proceedings of Machine Learning Research,pages 377–385, Thagaste, Ghent, Belgium, 2019. PMLR.64. M. Troffaes. Imprecise Monte Carlo simulation and iterative importance sampling for theestimation of lower previsions. International Journal of Approximate Reasoning, 101:31–48,2018.65. S. Vansteelandt, E. Goetghebeur, M. Kenward, and G.Molenberghs. Ignorance and uncertaintyregions as inferential tools in a sensitivity analysis. Statistica Sinica, 16:953–979, 2006.5 Statistics with Imprecise Probabilities—A Short Survey 7966. P. Walley. Statistical Reasoning with Imprecise Probabilities. Chapman &Hall, London, 1991.67. P. Walley. Inferences from multinomial data: Learning about a bag of marbles. Journal of theRoyal Statistical Society. Series B (Methodological), 58:3–57, 1996.68. P. Walley. A bounded derivative model for prior ignorance about a real-valued parameter.Scandinavian Journal of Statistics, 24:463–483, 1997.69. P. Walley. BI Statistical Methods: Volume I. Lightning Source, 2015.70. G. Walter, L. Aslett, and F. Coolen. Bayesian nonparametric system reliability using sets ofpriors. International Journal of Approximate Reasoning, 80:67–88, 2017.71. G. Walter and T. Augustin. Imprecision and prior-data conflict in generalized Bayesian infer-ence. Journal of Statistical Theory and Practice, 3:255–271, 2009.72. K. Weichselberger. Elementare Grundbegriffe einer allgemeineren Wahrscheinlichkeitsrech-nung I: Intervallwahrscheinlichkeit als umfassendes Konzept. Physica, 2001.73. K. Weichselberger. The theory of interval probability as a unifying model for uncertainty.International Journal of Approximate Reasoning, 24:149–170, 2000.74. K. Weichselberger. Elementare Grundbegriffe einer allgemeineren Wahrscheinlichkeitsrech-nung II. Symmetrische Wahrscheinlichkeitstheorie. (Manuscript), 2013.75. M. Zaffalon and E. Miranda. Conservative inference rule for uncertain reasoning under incom-pleteness. Journal of Artificial Intelligence Research, 34:757–821, 2009.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 6ReliabilityLisa Jackson and Frank P. A. CoolenAbstract This chapter introduces key concepts for quantification of system relia-bility. In addition, basics of statistical inference for reliability data are explained, inparticular, the derivation of the likelihood function.6.1 IntroductionReliability of systems is of crucial importance in all aspects of human life. Systemsare understood to be groupings of components in specific structure, with the systemfunctioning depending on the functioning of the components and the system struc-ture. Uncertainty about the functioning of components leads to uncertainty aboutthe system reliability. To study how system reliability depends on the reliabilityof components, several leading methods from the engineering literature are brieflyintroduced in Sect. 6.2. In situations of uncertainty about reliability in engineering,appropriate statistical methods are required to deal with the specific nature of data.This chapter provides an overview of such basic methods. Section6.3 introduces thekey statistical concepts, basic statistical models are presented in Sect. 6.4. Through-out the emphasis is on explanation of the likelihood function, which is at the heartof most statistical inference approaches as commonly used in reliability applica-tions. Unknown model parameters can conveniently be estimated by maximisationof the likelihood function, with corresponding theory to assess the uncertainty ofthe estimates. Bayesian methods for statistical inference are based on the likelihoodfunction as well, hence understanding of the likelihood function is crucial for studyof system reliability under uncertainty. Stochastic process models are also crucial toL. Jackson (B)Aeronautical and Automotive Engineering, Loughborough University, Loughborough,United Kingdome-mail: l.m.jackson@lboro.ac.ukF. P. A. Coolen (B)Department of Mathematical Sciences, Durham University, Durham, United Kingdome-mail: frank.coolen@durham.ac.uk© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_68182 L. Jackson and F. P. A. Coolendescribe variable reliability over time, for example, to reflect the effects of mainte-nance on a system’s reliability. A short introduction to such models is presented inSect. 6.5, again with an emphasis on deriving the likelihood function to enable sta-tistical inference. As this chapter brings together generic introductory material, nospecific references are included throughout the text. Instead, the chapter is endedwitha brief list of useful resources, pointing to a number of important books which arehighly recommended for further reading, and some brief comments about journalsin the field.6.2 System Reliability MethodsWhen modelling systems there are a number of tools, ranging from combinatorialmethods including reliability block diagrams, fault tree analysis and event tree anal-ysis, to more complex methods that cater for a greater range of system characteristicsincluding dependencies, e.g.Markov orMonte Carlo simulationmethods. Notewhenmodelling systems where repair times are involved, these typically do not follow theExponential distribution, e.g. the Lognormal or Weibull distributions may be moresuitable. Also including maintenance teams can mean dependencies are introducedvia prioritising strategies. In such instances, more complex methods are required. Abrief introduction to some of these modelling methods is provided in this section.6.2.1 Fault Tree AnalysisOne of themost common quantification techniques is fault tree analysis. This methodprovides a diagrammatic description of the various causes of a specified systemfailure in terms of the failures of its components. The choice of the system failuremode often follows from a failure mode and effects analysis (FMEA). There isthe assumption of independence of failures of the components. Logical gates areused to link together events (intermediate events shown as rectangles and basicevents representing failure events as circles), where the more common gates includeAND or OR, shown in Fig. 6.1, with an example tree shown in Fig. 6.2. Evaluationof the tree using Boolean algebra yields minimal cut sets, denoted as Ci , whichare failure combinations of components that are necessary and sufficient to causefailure of the system. Application of kinetic tree theory and the inclusion-exclusionprinciple (Eq.6.1) enables the system unavailability (Qsys) performance measure tobe calculated, where P(Ci ) is the probability of failure of minimal cut set Ci .6 Reliability 83Fig. 6.1 Common fault tree gate typesFig. 6.2 Example fault tree structure84 L. Jackson and F. P. A. CoolenQsys =Nc∑i=1P(Ci ) −Nc∑i=2i−1∑j=1P(Ci ∩ C j ) +Nc∑i=3i−1∑j=2j−1∑k=1P(Ci ∩ C j ∩ Ck) − . . .+(−1)Nc+1P(C1 ∩ C2 ∩ . . . ∩ CNc). (6.1)If the unavailability of the systemdoes notmeet performance acceptability criteria,then the system must be redesigned. An indication of where to make changes in thesystem can be achieved by generating component importance measures. This is ameasure of the contribution that each component makes to the system failure. Onesuch measure is the Fussel–Vesely measure of importance (IFVi ), defined as theprobability of the union of the min cut sets containing each component given thatthe system has failed, as shown in Eq.6.2.IFVi =P(⋃{C j |i ∈ C j })Qsys(6.2)6.2.2 Fault Tree Extensions: Common Cause FailuresTypically, fault trees have only a limited capability to cater for dependencies withinsystems. One example is that of common cause failures. Safety systems, for example,often feature redundancy, incorporated such that they provide a high likelihood ofprotection. However, redundant sub-systems or components may not always failindependently. A single common cause can affect all redundant channels at the sametime, examples include ageing (all made at the same time from the same materials),system environment (e.g. pressure or stress related) and personnel (e.g. maintenanceincorrectly carried out by the same person). There are several methods to analysesuch occurrences, including beta factor, limiting factor, boundary method and alphamethod. The beta factor method assumes that the common cause effects can berepresented as a proportion of the failure probability of a single channel of themulti-redundant channel. Hence, it assumes that the total failure probability (QT ) ofeach component is divided into two contributions: (i) the probability of independentfailure, QI , and, (ii) the probability of common cause failures, QCCF . The parameterβ is defined as the ratio of the common cause failure probability to the total failureprobability, as shown in Eq.6.3.β = QCCFQCCF + QI =QCCFQT. (6.3)There are further extensions to traditional fault tree analysis to cater for a greaterrange of dependencies, e.g. dynamic fault trees. As the nature of the dependencybecomes more complex, other modelling methods will be more suitable. Other typesof dependency include standby redundancy, where the probability of failure of the6 Reliability 85redundant component may change when it starts to function and experience load,hence it is dependent on the operating component and hence failure of both is notstatistically independent. Another form is multiple state component failure modes,where a component can exist in more than one state, hence being mutually exclu-sive cannot be considered in the fault tree (not failing in the open mode does notmean that it works successfully). In such instances, using Markov modelling meth-ods may be desirable. Given the state-space explosion that can exist when usingMarkov approaches, if there are just small subsections of the system exhibiting thedependency it may be possible to analyse these subsections with Markov methodsand use the result embedded in the fault tree approach.6.2.3 Phased Mission AnalysisAs systems become more complex, they can be required to undertake a numberof tasks, typically in sequence. An example might include an aircraft flight, whereit is required to taxi from the stand, take off, climb to the required altitude, cruise,descend, land and taxi to newdestination stand. The collection of tasks can be referredto as a mission, where each task is denoted by a phase which has an associated timeperiod. Mission success requires successful completion of all phases and in additiontheremay be different consequences resulting from the failure in each phase. For eachphase, an appropriate modelling technique can be employed to assess its reliability oravailability. When the systems phases are non-repairable then fault tree analysis canbe used to assess mission and phase success. For non-repairable scenarios, Markovor Petri nets can be used. The parameter of interest is the mission reliability. It is notappropriate to analyse the reliability of each phase and multiply these together to getthe mission reliability because (i) the phases are not independent (i.e. failure in onephasemay influence failure in another phase); (ii) the assumption that all componentsare working at the start of the phase is not correct and (iii) the system can fail on aphase change. For the non-repairable case, the initial step is to construct a missionfault tree. The general form is as shown in Fig. 6.3. The top event is a mutuallyexclusive OR gate, indicating that only one of the input events must happen to causethe output event (mission failure). When considering failures in phase 2 onwards,you need to include in the fault tree that the mission has been successful up to thispoint, namely, that it has functioned in the preceding phases. For this reason, you cansee the introduction of the NOT gate, i.e. under the intermediate event ‘Functions inPhase 1’. Alongside this the failure of the component in earlier phases also needs to betaken into account, hence the component failure is represented as shown inFig. 6.4. Toperform the analysis, both qualitatively and quantitatively, new algebra is required asshown in Fig. 6.5.C j corresponds to the failure of componentC in phase j , where thebar above C j corresponds to the working state. Considering the success of previousphases i = 1, . . . , j − 1, for failure in phase j , makes the analysis non-coherent,yielding prime implicant sets from a qualitative analysis (necessary and sufficientcombinations of events (success and failure)). The size of the problem to be solved can86 L. Jackson and F. P. A. CoolenFig. 6.3 Mission failure fault treeFig. 6.4 Revised componentrepresentation(and usually does) become prohibitively expensive, so approximations are required;these are usually based on coherent approximations of the non-coherent phases (i.e.conversion of prime implicants to minimal cut sets). Approximate quantificationformulae (e.g. Rare Event, Minimal Cut Set Upper Bound) can then be used.When analysing repairable systems, there are two requirements for mission suc-cess: (1) the system must satisfy the success requirements throughout each phaseperiod and (2) at the phase change times the system must occupy a state which issuccessful for both phases involved. The second point implies that we will con-sider failures on phase transition when calculating mission reliability. Analysis withrepairable systems is very similar in terms of generating themissionmodel and phase6 Reliability 87Fig. 6.5 Additional algebraFig. 6.6 Petri Net representation of a phased mission with three phasesmodels. This can be achieved using Markov or Petri Net methods. An illustrationof a Petri Net example is shown in Fig. 6.6, where the circles represent places andthe rectangular boxes represent transitions. Mission and phase reliabilities can beobtained from analysis of the model.6.3 Basic Statistical Concepts and Methods for ReliabilityDataConsider a random quantity T > 0, often referred to as a ‘failure time’ in relia-bility theory, but it can denote any ‘time to event’. Relevant notation for charac-teristics of its probability distribution includes the cumulative distribution function(CDF) F(t) = P(T ≤ t), the survival function S(t) = P(T > t) = 1 − F(t) (alsocalled the reliability function and denoted by R(t)), the probability density function(PDF) f (t) = F ′(t) = −S′(t) and the hazard rate h(t) = f (t)/S(t). The hazardrate has a possible interpretation with conditioning on surviving time t , for smallδt > 0, h(t)δt ∼ P(T ≤ t + δt | T > t). Harder to interpret, but also of use, is thecumulative hazard function (CHF) H(t) = ∫ t0 h(x)dx . Assuming R(0) = 1, we getH(t) = ∫ t0 f (x)S(x) dx = − ln S(t), so S(t) = exp{−H(t)} = exp{−∫ t0 h(x)dx}.88 L. Jackson and F. P. A. CoolenA constant hazard rate, h(t) = λ > 0 for all t > 0, gives S(t) = e−λt , the Expo-nential distribution. This can be interpreted as modelling ‘no ageing’, that is, if anitem functions, then its remaining time until failure is independent of its age. Thisproperty is unique to the Exponential distribution. An increasing hazard rate mod-els ‘wear-out’, roughly speaking this implies that an older unit has shorter expectedresidual life, and decreasing hazard rate models ‘wear-in’, implying that an olderunit has greater expected residual life. It is often suggested that ‘wear-out’ is appro-priate to model time to failure of many mechanical units, whereas electronic units’times to failures may be modelled by ‘wear-in’. In a human-life analogy, we canperhaps think about ‘wear-in’ as modelling time to death at very young age (‘infantmortality’) and ‘wear-out’ as modelling time to death at older age, with a periodin between where death is mostly ‘really random’, e.g. caused by accidents. This‘human-life analogy’ should only be used for general insight, and is included hereas engineers often claim that ‘typical hazard rates’ for components over their entirepossible lifetime are decreasing early on, then remain about constant for a reasonableperiod, and then become increasing (‘bath-tub shaped’).A popular parametric probability distribution for T is defined by the hazard rateh(t) = αβ(αt)β−1, for α, β > 0. This leads to S(t) = exp {−(αt)β}, and is called aWeibull distribution with scale parameter α and shape parameter β. This distributionis often used in reliability, due to the simple form for the hazard rate. For example,with β = 2 it models ‘linear wear-out’ (‘twice as old, twice as bad’).An interesting aspect of reliability data is that these are often affected by censoring,in particular, so-called right-censoring. Thismeans that, instead of actually observinga time at which a failure occurs, the information in the data is a survival of a certainperiod of time without failing. Clearly, such information must be taken into account,as neglecting it would lead to underestimation of expected failure times.Two main statistical methodologies use the likelihood function, namely, Bayesianmethods and maximum likelihood estimation. Hence, derivation of the likelihoodfunction is an important topic in reliability inference. Let t1, . . . , tn be observed fail-ure times, and c1, . . . , cm right-censored observations. For inference on a parameterθ of an assumed parametric model, the likelihood function based on these data isL(θ |t1, . . . , tn; c1, . . . , cm) = ∏nj=1 f (t j |θ)∏mi=1 S(ci |θ). This actually requires theassumption that the censoring mechanism is independent of the data distribution, ifthat is not the case the dependence would need to be modelled. It is also possible toconsider the likelihood over all possible probability distributions, so not restrictingto a chosen parametric model. In this case, the maximum likelihood estimator is theso-called Product-Limit (PL) estimator, presented by Kaplan andMeier in 1958. Thetheory of counting processes also provides a powerful framework for nonparametricanalysis of failure time data, based on stochastic processes and martingale theory.A well-known result within this theory is the Nelson–Aalen estimator for the CHF,which can be regarded as an alternative to the PL estimator.6 Reliability 896.4 Statistical Models for Reliability DataThis is a very wide topic, we can only mention a few important models. We first con-sider regression models for reliability data. Regression models are generally popularin statistics, and also very useful in reliability applications, where often Weibullmodels are used, with the survival function depending on a vector of covariates x ,and given by S(t; x) = exp{−(tαx)ηx}. Some simple forms are often used for theshape and scale parameters as functions of x , e.g. the loglinear model for αx , spec-ified via ln αx = xTβ, with β a vector of parameters, and similar models for ηx .The statistical methodology is then pretty similar to general regression methods, andimplemented in statistical software packages. Suchmodels need to be fully specified,so are less flexible than nonparametric methods, but they allow information in theform of covariates to be taken into account.Semi-parametric models enable covariates to be taken into account, but do sowithout fully specifying a parametric model, keeping some more flexibility. Usually,a parametric form for the effect of the covariates on a nonparametric ‘baselinemodel’is assumed. Most famous are the Proportional Hazards (PH) models, presented byCox in 1972. Here, the hazard rate for covariates x is defined by h(t; x) = h0(t)ψx ,with h0(t) the baseline hazard rate (normally left unspecified, so nonparametric),and ψx some positive function of x , independent of time t (normally a fully para-metric form is assumed for ψx ). The name of such models results from the fact thath(t;x1)h(t;x2) =ψx1ψx2, independent of t , so the hazard rates corresponding to different covari-ates are in constant proportion. For these models, ln S(t; x) = − ∫ t0 h(u; x)du =−ψx∫ t0 h0(u)du = ψx ln S0(t), so S(t; x) = [S0(t)]ψx . These models are used mostfor survival data in medical applications, but are also common and useful in reli-ability. As there are no assumptions on the form of the baseline hazard rate, theyprovide a valuable method to compare the effect of the covariates. An often used PHmodel is the linear PH model, with ψx = exp{xTβ}, with β a vector of parameters.We now describe the analysis of this particular model, with no further assumptionson h0(t). The goal is to estimate β and R0(t), the baseline survival function relatedto h0(t). This is far from trivial, as it is not clear how the likelihood function can bederived, since this is neither uniquely defined by a fully specified parametric model,nor completely free as was the case for fully nonparametric models (leading to thePL estimate). Hence, we need to use a different concept.Suppose we have data on n items, consisting of r distinct event times, t(1) 0, for i = 1, . . . , n, be very small. The process observed canthen be described as consisting of: 0 failures in (0, t1), and 1 failure in [t1, t1 + δt1),and 0 failures in [t1 + δt1, t2), etc., until no failures in [tn + δtn, r ] (this last bit is justdeleted if r = tn , so when observation of the process is ended at the moment of then-th failure). To derive the corresponding likelihood function, we take the productof the probabilities for these individual events, soexp{−m(0, t1)} × m(t1, t1 + δt1) exp{−m(t1, t1 + δt1)} × exp{−m(t1 + δt1, t2)} × . . .. . . × exp{−m(tn + δtn, r)}={n∏i=1[∫ ti+δtitiv(t)dt]}× exp[−∫ r0v(t)dt].Now, use that for very small δti , we have that∫ ti+δtitiv(t)dt ≈ v(ti )δti .Now we divide through by∏ni=1 δti , and let all δti ↓ 0 (this is, exactly the same thatleads to the PDFs appearing in the likelihood function for iid data). This gives thelikelihood function, for this model and based on these n failure data and observationover the period [0, r ]:L ={n∏i=1v(ti )}exp[−∫ r0v(t)dt].For optimisation, it is easier to use the log-likelihood function, which is also neededfor related statistical inference, and which is equal to:l =n∑i=1ln v(ti ) −∫ r0v(t)dt.It is possible to work with this likelihood non-parametrically, but often one assumesa parametric form for the ROCOF, making maximum likelihood estimation againconceptually straightforward (although it normally requires numerical optimisation).Two simple, often used parametric ROCOFs arev1(t) = exp(β0 + β1t)6 Reliability 93andv2(t) = γ ηtη−1,with γ, η > 0.Many models that have been suggested, during about the last three decades, forsoftware reliability, are NHPPs which model the software testing process as a faultcounting process. A famous model was proposed by Jelinski and Moranda in 1972,which is based on the following assumptions: (1) software contains an unknownnumber of bugs, N ; (2) at each failure, one bug is detected and corrected; (3) theROCOF is proportional to the number of bugs present. So, they use a NHPP withfailure times Ti , i = 1, . . . , N and T0 = 0, defined byv(t) = (N − i + 1)λ, for t ∈ [Ti−1, Ti ),for some constant λ. Then N and λ are both considered unknown, and estimated fromdata, where, of course, inference for N tends to be of most interest, or, in particular,the number of remaining bugs. Many authors have contributed to such theory bychanging somemodel assumptions. For example, non-perfect repair of bugs has beenconsidered, and even the possibility of such repair introducing new bugs (possibly arandom number), for this last situation so-called ‘birth-death processes’ can be used.Also non-constant λ has been considered, e.g. with the idea that some bugs maytend to show earlier than others. Also Bayesian methods for such models, and evensoftware reliability models more naturally embedded in Bayesian theory, have beensuggested and studied. However, although there is an enormous amount of literaturein this area, as indeed mathematical opportunities appear to have no limit here, thepractical relevance of such models seems to be rather limited and few interestingapplications of such models have been reported in software reliability. Recently, theimportant topic of testing of reliability of systems including software has receivedincreasing attention, which is much needed to ensure reliable systems.Useful sourcesP.K. Andersen, O. Borgan, R.D. Gill and N. Keiding, Statistical Models Based onCounting Processes (Springer, 1993).T. Aven and U. Jensen, Statistical Models in Reliability (Springer, 1999).R.E. Barlow, Engineering Reliability (SIAM, 1998).R.E. Barlow and F. Proschan, Mathematical Theory of Reliability (Wiley, 1965).T. Bedford and R. Cooke, Probabilistic Risk Analysis: Foundations and Methods(Cambridge University Press, 2001).P. Hougaard, Analysis of Multivariate Survival Data (Springer, 2000).R.S. Kenett, F. Ruggeri and F.W. Faltin (Eds), Analytic Methods in Systems and Soft-ware Testing (Wiley, 2018).J.F. Lawless, Statistical Models and Methods for Lifetime Data (Wiley, 1982).R.D. Leitch, Reliability Analysis for Engineers (Oxford University Press, 1995).H.F. Martz and R.A. Waller, Bayesian Reliability Analysis (Wiley, 1982).94 L. Jackson and F. P. A. CoolenW.Q. Meeker and L.A. Escobar, Statistical Methods for Reliability Data (Wiley,1998).N.D. Singpurwalla and S.P. Wilson, Statistical Methods in Software Engineering:Reliability and Risk (Springer, 1999).Leading international journals in this field include Reliability Engineering andSystem Safety, IEEE Transactions on Reliability, Journal of Risk and Reliability,Quality and Reliability Engineering International. Statistical methods for reliabilitydata are presented in a wide variety of theoretical and applied Statistics journals, the-ory and methods for decision support are often published in the Operations Researchliterature.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 7Simulation Methods for the Analysisof Complex SystemsHindolo George-Williams, T. V. Santhosh, and Edoardo PatelliAbstract Everyday systems like communication, transportation, energy and indus-trial systems are an indispensable part of our daily lives. Several methods have beendeveloped for their reliability assessment—while analytical methods are computa-tionally more efficient and often yield exact solutions, they are unable to account forthe structural and functional complexities of these systems. These complexities oftenrequire the analyst tomakeunrealistic assumptions, sometimes at the expenseof accu-racy. Simulation-based methods, on the other hand, can account for these realisticoperational attributes but are computationally intensive and usually system-specific.This chapter introduces two novel simulation methods: load flow simulation andsurvival signature simulationwhich together address the limitations of the existinganalytical and simulation methods for the reliability analysis of large systems.7.1 IntroductionA system is classed as complex from one of two fronts—in terms of the functionalrelationships between its components and in terms of its structure. A structurallycomplex systemdoes not conform to a series, parallel, or series-parallel configuration.Most real-world systems are composed of components that can operate at multipleH. George-WilliamsInstitute of Energy and Sustainable Development, School of Engineering and SustainableDevelopment, De Montfort University, Oxford, United Kingdome-mail: hindolo.george-williams@dmu.ac.ukT. V. SanthoshInstitute for Risk and Uncertainty, University of Liverpool, Liverpool, United Kingdome-mail: s.santhosh@liverpool.ac.ukBhabha Atomic Research Centre, Mumbai, IndiaE. Patelli (B)Centre for Intelligent Infrastructure, Civil and Environmental Engineering,University of Strathclyde, Glasgow, United Kingdome-mail: edoardo.patelli@strath.ac.uk© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_79596 H. George-Williams et al.performance levels or states and components with a functional coupling with othercomponents. Such systems are deemed functionally complex, since their states cannotbe directly deduced from their traditional two-state structure functions. They arecharacterised bymultiple states, with the number of states determined by the diversityin the states of their components, structure and the functional relationships betweentheir components [21]. In these systems, the number of performance levels mayor may not be finite, depending on the performance measure under considerationand the type of system [21]. For instance, the power generated by a power plantmay take any value between zero and its maximum achievable value, dependingon the performance levels of its component and the demand on the grid. Complexsystems may be standalone or form an indispensable part of some critical system likehealthcare, safety-critical and industrial control systems. It is, therefore, importantto be able to assess their susceptibility to failures, as well as quantify and predict theensuing consequences, for effective planning of restoration andmitigation measures.7.2 Reliability Modelling of Systems and NetworksIn system reliability evaluation, the analyst has numerous techniques at their disposal,which can be classified as heuristic-, analytical- or simulation-based [1] and furtheras static or dynamic. In particular, dynamic techniques not only model the systembased on the functional and structural relationships between its components, but alsosupport dynamic relationships like inter-component and inter-system dependencies.7.2.1 Traditional ApproachesReliability Block Diagrams and Fault Trees have been extensively used in the reli-ability evaluation of binary-state systems. Both techniques have proven particularlyuseful for moderately sized systems with series-parallel configurations. However,they become difficult to apply with large or complex systems and often requireadditional techniques to decompose the system. The Reliability Graph [40] was,therefore, developed to overcome this difficulty and proved very efficient in mod-elling structural complexities. Reliability block diagrams, fault trees and reliabilitygraphs, however, assume components to be statistically independent, which rendersthem inadequate for systems susceptible to restrictivemaintenance policies and inter-component dependencies. However, techniques including but not limited to dynamicreliability block diagrams [10], dynamic fault trees [6], condition-based fault trees[35], dynamic flow graphs [2], Petri Nets [26] and other combinatorial techniques[38] have been developed to model these dynamic relationships. They have foundapplication in a wide range of reliability engineering problems, including repairablesystems with restrictive maintenance policies.7 Simulation Methods for the Analysis of Complex Systems 97Though the earliest forms of these techniques including binary decision diagramswere applicable only to binary-state systems, numerous instances of their recentextension tomulti-state systems exist, see, e.g. [39]. However, these extensions eitherrequire state enumeration or the derivation of the minimal path or cut sets of thesystem, which is an NP-hard problem [41].The extended block diagram technique and graph-based algorithms sharetwo common limitations. First, they define reliability with respect to the maximumflow through the system. Therefore, they are limited to systems with single outputnodes and those with multiple output nodes where only the presence of flow at thesenodes is relevant and not the relative magnitude of the flow. The second limitationarises from the assumption that there are no flow losses in the system, making theminapplicable to certain practical systems like energy systems and pipe networks,susceptible to losses in some failure modes. More recently, various researchers havemade invaluable contributions to multi-state system reliability analysis, developingtechniques applicable to a wide range of systems [22]. These techniques have mainlybeen based on either the structure function approach, stochastic process, simulationor the Universal Generating Function approach [21, 25].The most popular stochastic process employed in system reliability analysis isthe Markov Chain (MC), which involves enumerating all the possible states of thesystem and evaluating the associated state probabilities [25]. This technique is onlyeasily applicable to exponential transitions or distributions with simple cumulativedistribution functions, requires complicated mathematics and becomes complex forlarge systems. For an M component binary-state system, the number of states inthe model ranges from M + 1 for series systems, to 2M for parallel systems. Forlarge multi-state systems, the number of states increases exponentially, rendering themodel difficult to construct and expensive to compute.The Universal Generating Function was introduced to address the state explo-sion problem of the MC. It allows the algebraic derivation of a system’s perfor-mance from the performance distribution of its components [21, 24]. However, boththe Universal Generating Function and Markov Chain are limited in the numberof reliability indices they can quantify. Also, like all multi-state system reliabilityevaluation techniques, they are maximum-flow-based and assume flow conservationacross components. The Universal Generating Function, though straightforward forseries/parallel systems, it requires a substantial effort for complex topologies.Simulation methods are the most suitable for multi-state system reliability andperformance evaluation, since they mimic the actual operation of systems. Theiradvantage over their analytical counterpart is due to the fact that they support anytransition distribution, allow the effects of external factors on system performance tobe investigated [43] and are easily integrated with other methods [36]. In particular,they allow the explicit consideration of the effects of uncertainty and imprecisionon the system, providing a powerful tool for risk analysis and by extension, ratio-nal decision-making under uncertainty. They are, therefore, mostly used to analysesystems for which analytical approaches are inadequate. However, even some of theexisting simulation methods [23, 43] require prior knowledge of the system’s pathset, cut set or structure function and are mostly limited to binary-state systems [42].98 H. George-Williams et al.7.2.2 Interdependencies in Complex SystemsEngineers and system designers are under immense pressure to build systems robustand adequate enough to meet the ever-increasing human demand and expectation.Unavoidably, the resulting systems are complex and highly interconnected, whichironically constitute a threat to their resilience and sustainability [18]. Two systemsare interdependent if at least a pair of components (one from each system) are cou-pled by some phenomena, such that a malfunction of one affects the other. In suchsystems, an undesirable glitch in one system could cascade and cause disruptionsin the coupled system. The cascade could be fed back into the initiating systemand the overall consequences may be catastrophic [5]. To minimise the effects offailures, some interdependent systems are equipped with reconfiguration provisions.This normally entails transferring operation to another component, rerouting flowthrough alternative paths, or shutting down parts of the system.The achievement of maximum overall system performance is, in general, desir-able. However, in many applications (nuclear power plants, for instance), it is moreimportant to guarantee system availability and recovery in the shortest possible time,following component failure [16]. Interdependencies are manifested in engineeringsystems at two levels: between components (inter-component), which can be func-tional or induced and between systems/subsystems (inter-system) [15].Functional dependencies are due to the topological and/or functional relationshipsbetween components. Induced dependencies, on the other hand, are due to a statechange in one component (the initiator) triggering a corresponding state change inanother (the induced), such that evenwhen the initiator is reinstated, the induced doesnot reinstate, unless manually made to do so. Functional dependencies in standalonesystems are intrinsically accounted for by the innate attributes of the system reliabilitymodelling and evaluation technique while induced dependencies require explicitmodelling. Inter-system dependencies, on the other hand, are due to functional orinduced couplings between multiple systems. The functional dependencies in thesesystems, however, may require explicit modelling. This is the case for componentsrelying on material generated by another system. For instance, an electric pump in awater distribution system relies on the availability of the electricity network.Induced dependencies are further divided into Common-Cause Failures (CCF)[27] and cascading events, as summarised in Fig. 7.1. Common-cause failures arethe simultaneous failure of multiple similar components due to the same root cause.Their origin is traceable to a coupling that normally is external to the system. Notableinstances are shared manufacturing lines, shared maintenance teams, shared envi-ronments and human error. A group of components susceptible to the same CCFevent is called a Common-Cause Group (CCG). An important point to note aboutcommon-cause failures is that, on occurrence of the failure event, there is a probabil-ity associated with multiple component failure and that the affected components failin the same mode. Consequently, the number of components involved in the eventranges from 1 to the total number of components in the CCG. CCF events may affectan entire system or only a few of its components and, therefore, pose a consider-7 Simulation Methods for the Analysis of Complex Systems 99InterdependenciesInter-component Inter-systemFunctional InducedCascading FailuresCommon-Cause FailuresFig. 7.1 Types of interdependencies in complex systems. Functional dependencies—such as whenthe failure of power supply forces the unavailability of connected components. Common-CauseFailures—due to earthquake excitation, vibration, environmental conditions (temperature, humidity,contaminants), shared maintenance. Cascading events such as the failure of one component mightoverload other componentsable threat to the reliability of systems. CCF modelling and quantification attractskeen interest from system reliability and safety researchers, as well as practitioners.Examples of the work that has been done in this field can be found in [28, 33, 37].Most of the methods presented in these publications, however, are built on reliabilityevaluation techniques that do not segregate the topological from the probabilisticattributes of the system. As such, they are computationally expensive for problemsinvolving multiple reliability analysis of the same system. They also have yet to beapplied to multi-state systems, as well as systems susceptible to both cascading andcommon-cause failures.Cascading failures are those with the capacity to trigger the instantaneous failureof one ormore components of a system. They can originate from a component or froma phenomenon outside the system boundary. The likelihood of the initiating eventoriginating from within the system distinguishes them from CCF. Another point ofdichotomy is that the affected components do not necessarily have to be similar orfail in the same mode. In addition, at the occurrence of the initiating event, the prob-ability of all the coupled components failing is unity, same for the case when they arein a state rendering them immune [15, 18]. A few prominent examples of initiatingevents external to the system are extreme environmental events, natural disasters,external shocks, erroneous human-system interactions and terrorist acts. Variousmodels have been developed to study the effects of cascading failures on complexsystems [29]. However, a good number of these models only assess their response100 H. George-Williams et al.to targeted attacks, variation in some coupling factor or the relative importance ofsystem components. When faced with the additional situation of random componentfailures, a complete reliability and availability analysis should be performed [18].Even methods that fulfill this requirement have their applicability hampered by com-ponents that undergo non-Markovian transitions, components susceptible to delayedtransitions, and reconfigurable systems.7.3 Load Flow SimulationThe load flow simulation is a recently proposed technique for the reliability andperformance analysis of multi-state systems [17]. It is based on the fact that if theperformance levels of a system’s components are known, the performance levelsof the system can be directly derived from its network model. In this formalism,each component is modelled as a semi-Markov stochastic process and the systemas a directed graph whose nodes are the components of the system. The approachis intuitive and applicable to any system architecture and easily programmable ona computer. It outperforms other multi-state system reliability analysis approaches,since it does not require state enumeration or cut set definition. Efficient algorithmsfor manipulating the adjacency matrix of this directed graph to obtain the flow equa-tions of the system are available in OpenCossan [31].The operation of the system is simulated using Kinetic Monte Carlo method byinitially sampling the state and time to the next transition (hereafter referred to astransition parameters) of each component. The simulation jumps to the smallest sam-pled transition time tmin , at which time the states of the components undergoing thetransition are updated. Using the updated performance levels of the components ofthe system, the virtual flow across the system is computed via a linear programmingprocedure that employs the interior-point algorithm. The new transition parametersof the components undergoing a transition are then sampled and the simulation jumpsto the next smallest transition time. This cycle of component transition parametersampling, transition forcing and system performance computing continues until themission time T is reached. The system performance computed at every componenttransition is captured and saved in counters, from which the performance indices ofthe system can be deduced. A component shutdown and restart procedure is incor-porated to replicate the actual operating principles of most practical systems. In thisprocedure, the availability of each system component is tested against its predefinedreference minimum input load level at every transition and the effects of functionalinterdependence on the failure probability of the components are accounted for.Figure7.2 provides a high-level illustration of the load flow simulation procedure.Ageing and component performance degradation is common in most systems.For such systems, techniques built around the flow conservation principle becomeobsolete, as the flow generated by sources can be dissipated in intermediate com-ponents in certain failure modes. For instance, consider a 100MW power generatorsupplying a 95MW load through a 125MW transformer. If there are no power losses7 Simulation Methods for the Analysis of Complex Systems 101StopStartInitialisationSet time t=0Compute systemperformanceSample componentstransition timeFind smallest timeSet t=tmin t > TMissionCompute systemreliability and performanceSample component status update system configurationFig. 7.2 Flowchart of the load flow simulationin the transformer, 95MWwill be drawn from the generator and delivered to the load.However, if the efficiency of the transformer deteriorates to say 75%, it now takesall 100MW from the generator but delivers only 75MW to the load. In both cases,the apparent difference between the generation capacity and demand is the same butthe power drawn from the generator increases while the effective power suppliedto the load deteriorates. For this example, the demand would have to be slashed to75MW or less, to preserve the operational integrity of the generator. Other scenarioswhere component inefficiency affects system reliability are: a power transmissionline prone to losses and an oil pipeline where a failure mode is a hole in a pipe orgasket failure at some flange [17].The load flow simulation approach has been successfully applied to the availabil-ity assessment of a reconfigurable offshore installation [18], dynamic maintenancestrategy optimization of power systems [19] and the probabilistic risk assessment ofstation blackout accidents in nuclear power plants [16].Advantages Over Existing Techniques:1. Inherits all the advantages of simulation approaches used for system reli-ability and performance evaluation.2. Implements any system structurewith relative ease, since it doesn’t requireknowledge of the minimal path or cut sets prior to system analysis.3. Calculates the actual flow across every node of the system.4. Models systemsmadeupofmultiple source and sinknodeswith competingstatic or dynamic demand.5. Models losses in components and across links.6. Models component restart and shutdown.7. Not limited to integer-valued node capacities and system demand, asrequired by other graph-based algorithms.7.3.1 Simulation of Interdependent and ReconfigurableSystemsLoad flow simulation allows the modelling of inter-component and inter-systemdependencies, thereby supporting the reliability assessment of realistic engineeringsystems [18]. Components and external events that influence the operation of the102 H. George-Williams et al.Fig. 7.3 Illustration of decoupling procedure for interdependent systemssystem are identified and numbered, followed by the identification and modelling ofall the inter-component dependencies in the system. The strategy is to decouple theinterdependent system into its constituent systems (subsystems) as shown in [18].The nodes associated with each subsystem are then identified and its directed graphobtained (i.e. only nodes with actual commodity flow are considered). The states ofeach node are then identified and modelled as described in [17].For illustrative purposes, consider the original system in Fig. 7.3 (left panel). It isan interdependent four commodity system—each solid line transports a commodityand the broken line depicts an induced dependency in the direction of the arrow.Node 2 is part of subsystem S2 and relies the commodity from subsystem S3 todrive its operation. One would say it is functionally dependent on subsystem S3 andexhibits a dual operation mode, operating both as a sink and an intermediate node.Its sink mode directly influences flow in S3, while its transmission mode directlyinfluences flow in S2. It is, therefore, logical to separate node 2 into its constituentnodes, each representing a mode of operation. Virtual nodes representing the sinkmodes of dual nodes are created and assigned new IDs, creating a decoupled system(see Fig. 7.3 (right panel)). A load-source functional dependency exists between thedecoupled nodes, since the transmission node is incapacitated if flow into the sinknode is inadequate. Therefore, they make a load-source pair, with the transmissionnode being the load and the sink node, the local source node.Local sources, otherwise known as support nodes in load-source pairs, are mod-elled as binary-state objects: state 1 (active) has capacity l, depicting the availabilityof the dependent node; State 2 (inactive) has capacity 0 and depicts its unavailabil-ity. l is the minimum level of support required to operate the dependent/sink nodeand in practical cases represents the load rating of that component. By applying thedecoupling procedure described to all load dependency relationships in the system,the following load-source pairs; {2, 14}, {3, 16}, {1, 18}, {13, 15} and {9, 17} areobtained. Li = { j, l} signifies that node i requires a minimum of l units of a certaincommodity from node j to operate. If i has a load dependency relationship with7 Simulation Methods for the Analysis of Complex Systems 103multiple nodes, Li takes the form of a two-column matrix, with each row definingthe node’s relationship with another node.Induceddependencies are definedby theparameter Di = {d j1, d j2, d j3, d j4}u×4 |j = 1, 2, ..., u − 1, u, which defines the state change induced in other nodes as aresult of a state change in node i . d j1 is the state of i triggering the cascading event,d j2; the affected node, d j3; the state the node has to be in to be affected, and d j4; itstarget state on occurrence of the event. Each row of Di defines the behaviour of anaffected node, and u, the number of relationships. If node i and the affected noded j2 belong to different subsystems, the subsystem the latter belongs to is dependenton the subsystem of the former. For example, suppose state 2 of node 5 in Fig. 7.3forces node 7 into state 3 if it is in state 1 at the time node 5 makes the transition tostate 2. The induced dependency of node 7 on node 5 is defined by D5 asD5 =(2 7 1 3)(7.1)Once the system has been decoupled, the dependency tree depicting the relation-ships between its subsystems and their ranking is derived. The rank of a subsystemdepends on its position on the tree relative to the reference subsystem. The indepen-dent subsystem, which is also the reference subsystem, is assigned rank 1 and theremainder ranked in ascending order of their longest distance from this reference.See [18] for the details of the ranking, reconfiguration and simulation procedures.7.3.2 Maintenance Strategy OptimizationThe load flow simulation approach can be exploited to optimise the maintenancestrategies of complex systems. The multi-state semi-Markov models of componentsare extended to represent their behaviour under various maintenance strategies. Theoperation of the system is then simulated using a slightly modified version of thesimulation procedure depicted in Fig. 7.2 and detailed in [19]. Non-Markovian com-ponent transitions associatedwith the operational dynamics imposed bymaintenancestrategies are implemented. For example, the maintenance of a failed component canonly be initiated if there is an idle maintenance team, making the transition of thecomponent from its failed to working state non-Markovian, since it is conditional onthe availability of a maintenance team. Additional component states such as preven-tive maintenance, corrective maintenance, shutdown, diagnostics, idle and awaitingmaintenance are included to model different maintenance activities.To illustrate the derivation of the multi-state model of a component under variousmaintenance strategies, consider a binary-state component. The component is subjectto both preventive and corrective maintenance and maintained by a limited numberof maintenance teams. In addition, its corrective maintenance consists of two stages:a diagnosis stage and a restoration stage. Following diagnosis, the maintenance teamcould proceed with the actual repairs if spares are not required or make a sparesrequest. There is a known probability associatedwith spares being needed for a repair104 H. George-Williams et al.Fig. 7.4 Multi-state models of binary-state component under maintenance delaysand while the maintenance team awaits the spares, it could be assigned to anothercomponent. Similarly, there is a probability associated with spares being needed tocomplete the preventive maintenance of the component, which could be interruptedif these spares are not immediately available. The resulting multi-state models ofthe component under two contrasting maintenance strategies are shown in Fig. 7.4,with the component’s state assignments and possible transitions. Transitions areeither normal, forced or conditional. Normal transitions occur randomly and dependonly on their associated time-to-occurrence distributions. Forced transitions occurpurely as a consequence of events outside the component boundary, and their time-to-occurrence distributions are unknown. Conditional transitions, on the other hand,have a known time-to-occurrence distribution but are assigned a lower priority andonly occur on fulfilment of a predefined probabilistic condition or set of conditions[19]. Unlike normal transitions in which the next state of the component dependsonly on its current state, the next state of the component under forced transitionsmay also depend on its previous state. As such, the multi-state component transitionparameter sampling procedure presented in [17] cannot be used to determine thetransition parameters of the component. For this, the set of procedures presented in[19] are required. The binary-state component models in Fig. 7.4 can be generalisedfor multi-state components by defining one ‘Idle’ state (if components are kept outof operation during spares delay), a ‘Diagnosis’ state (where necessary) and one‘Corrective Maintenance’ state for each repairable failure mode.With this approach, multiple contrasting complex maintenance strategies can besimulated without the need to modify the simulation algorithm, as the maintenancestrategy is implemented at the component level. See, for instance, the optimal main-tenance strategies for a hydroelectric power plant derived in [19].7 Simulation Methods for the Analysis of Complex Systems 1057.3.3 Case Study: Station Blackout Risk AssessmentThe complete lack of AC power at a nuclear power plant is critical to its safety, sinceAC power is required for its decay heat removal. Though designed to cope withthese incidents, nuclear power plants can only do so for a limited time. The impact ofstation blackouts on a nuclear power plant’s safety is determined by their frequencyand duration. These quantities, however, are traditionally computed via a static faulttree analysis that deteriorates in applicability with increasing system complexity.The load flow simulation approach was used to quantify the probability and dura-tion of possible station blackouts at the Maanshan Nuclear Power Plant in Taiwan,accounting for interdependencies between system components, maintenance, systemreconfiguration, operator response strategies and human errors [16].TheMaanshan Plant is powered through two physically independent safety buses,which themselves are powered by six offsite power sources through two independentswitchyards. Each safety bus has a dedicated backup diesel generator and both busesshare a third diesel generator. Two gas turbine generators connected through thesecond switchyard power the plant’s safety systems if all three diesel generators areunavailable. The gas turbine generators, however, take about 30 min to become fullyoperational, when powered on. The goal in this case study was to quantify the riskto the plant, of station blackouts initiated by the failure of the grid sources, as wellas the switchyards and identify the best recovery strategy, to minimise this risk.The load flow simulation approach was used to model the structural/functionalrelationships between the components of the system as described in Sect. 7.3 andthe formalism described in Sect. 7.3.1 to model both the interdependencies betweencomponents and their dynamic behaviour under various recovery strategies. The fulldetails of the solution approach and results are available in [16].7.4 Survival Signature SimulationFor very large-scale systems and networks, the full system structure information (orstructure function, minimal paths sets, etc.) might not be available or may be difficultto obtain. Having a compact representation of the system, therefore, is advantageous.Survival signature [7] has been proposed as a generalisation of system signature[11, 12] to quantify the reliability of complex systems consisting of independentand identically distributed (i id) or exchangeable components, with respect to theirrandom failure time. It has been shown in [8] how the survival signature can bederived from the signatures of two subsystems in both series and parallel config-uration. The authors developed a non-parametric-predictive inference for systemreliability using the survival signature. Aslett et al. [3] demonstrated the applicabil-ity of the survival signature to system reliability quantification via a parametric, aswell as non-parametric approach.An efficient computational approach for computingapproximate and exact system and survival signatures has been recently presentedin [20, 34]. Feng et al. [13] developed an analytical method to calculate the sur-106 H. George-Williams et al.Fig. 7.5 Example of a bridge network composed of six-component of two typesTable 7.1 Survival signature for the system shown in Fig. 7.5l1 l2 (l1, l2) l1 l2 (l1, l2)0 [0, 1, 2, 3] 0 2 [0, 1] 01 [0, 1] 0 2 2 1/31 2 1/9 2 3 2/31 3 1/3 3 [0, 1, 2, 3] 1vival function of systems with uncertainty in the parameters of component failuretime distributions. These methods are all useful but less practical for larger complexsystems and not applicable to non-exponential transitions.As an illustration, consider a six-component bridge network with two componenttypes (Fig. 7.5), the survival function is given by Table7.1.Considering 2 working components of type 1; l1 = 2 and 3 of type 2; l2 = 3,there are three possible combinations in total but only two combinations lead tosuccess (the survival of the system) of the system. Hence, the survival signature ofthe system is 23 , as shown in Table7.1. Similarly, for l1 = 3 and l2 = [0, 1, 2, 3],there are eight possible combinations in total, all of which result in success. Hence,the survival signature of the system in this case is equal to 1.0. Thus, knowing thesuccess paths from the combinations of multiple types of active components, it ispossible to compute the survival function of a complex system.Exact analytical solutions are restricted to particular cases (e.g. systems withcomponent failure times following the exponential distribution and non-repairablecomponents). The survival function of a system with K component types is given byP(Ts > t) =m1∑l1=0...mK∑lK =0φ(l1, . . . , lK )P(K⋂k=1{Ck(t) = lk}) (7.2)whereP(K⋂k=1{Ck(t) = lk}) =K∏k=1(mklk)[Fk(t)]mk−lk [1 − Fk(t)]lk (7.3)Here, Ck(t) ∈ {0, 1, . . . ,mk} denotes the number of components of type k in thesystemwhich function at time t , and Fk(t) represents theCDF of the random failuretimes of components of the different types. In this approach, we have a strong i id7 Simulation Methods for the Analysis of Complex Systems 107Fig. 7.6 Flow chart of the Monte Carlo simulation algorithm for complex systems with repairablecomponents based on survival signature. Details of the simulation method are available in [30]assumption of failure times within same components types. With this assumption,all state vectors [7] are equally likely to occur.However, simulation methods can be applied to study and analyse any system,without introducing simplifications or unjustified assumptions. AMonteCarlo-basedapproach can be combined with survival signature, to estimate the reliability ofa system in a simple and efficient way. A possible system evolution is simulatedby generating random events (i.e. the random transition such as failure times ofthe system components) and then estimating the status of the system based on thesurvival signature (Eq. (7.2)). By counting the number of occurrences of a specificcondition (e.g. the number of times the system is in working status), it is possible toestimate the survival function and reliability of the system.The most generally applicable Monte Carlo simulation methods adopting the sur-vival signature for multi-state component and repairable systems have been proposedin [30]. Its procedural steps are presented in Fig. 7.6.7.4.1 Systems with ImprecisionThe reliability analysis of complex systems requires the probabilistic characterisationof all the possible component transitions. This usually requires a large dataset that isnot always available. To avoid the inclusion of subjective assumptions, imprecisionand vagueness of the data can be treated by using imprecise probabilities that combineprobabilistic and set theoretical components in a unified construct (see, e.g. [4, 9]).Randomness and imprecision are considered simultaneously but viewed separatelyat any time during the analysis and in the results [32].108 H. George-Williams et al.Imprecision can occur at component level, where the exact failure distribution isnot known or at system level, in the form of an imprecise survival signature. The latteroccurs when part of the system can be unknown or not disclosed. Such imprecisionleads to bounds on the survival function of the system, providing confidence in theanalysis, in the sense that it does not make any additional hypothesis regarding tothe available information. When the imprecision is at the component level, a naïveapproach, employing a double loop sampling approach where the outer loop is usedto sample realisations of component parameters, can be used. In other words, eachrealisation defines a new probabilistic model that needs to be solved adopting thesimulationmethods proposed above, fromwhich the envelop of the system reliabilityis identified. However, since almost all systems are coherent (a system is coherent ifeach component is relevant, and the structure function is nondecreasing), it is onlynecessary to compute the system reliability twice, using the lower and upper boundsfor all the parameters, respectively. If the imprecision is at the system level (i.e. inthe survival signature), the simulation strategy proposed in Fig. 7.6 can be adoptedwithout additional computational cost by collecting, in two separate counters, theupper and lower bounds of the survival signature at each component transition, asillustrated in [30]. Hence, imprecision at the component and system levels can beconsidered concurrently, without additional computational costs.7.4.2 Case Study: Industrial Water Supply SystemAn industrial water supply system consisting of 13 components, as shown in Fig. 7.7,is chosen as a case study, to demonstrate the capability of the survival signaturemethod. The system is expected to deliver water to at least one of the two tanksT 2 or T 3 from tank T 1, through a set of motor-operated pumps and valves. Thecomponent failure datawith the corresponding distributions are provided inTable7.2.The survival signature method is employed to compute the reliability of the system.Fig. 7.7 Industrial water supply system7 Simulation Methods for the Analysis of Complex Systems 109Table 7.2 Reliability parameters of the components of the water supply systemComponent Failure rate(h−1)MTTR (h) Repair rate(h−1)DistributiontypeT1, T2, T3 λ1 = 5 · 10−5 24 μ1 = 0.0417 ExponentialP1, P2, P3 λ2 = 3 · 10−3 17.4 μ2 = 0.0575 ExponentialV1, V2, V3, V4, V5, V6, V7 λ3 = 2 · 10−4 9 μ3 = 0.111 ExponentialTable 7.3 Survival signature (selected parts only) for the system shown in Fig. 7.7 computed withapproach proposed in [20]l1 l2 l3 l1 l2 l3 [0, 1] ∀ ∀ 02 1 2 1/63 3 1 2 1/212 1 5 8/63 3 1 5 8/212 1 7 2/9 3 1 7 2/32 2 6 22/63 3 2 6 6/72 3 5 8/21 3 3 5 6/72 3 7 2/3 3 3 7 1The components of the system are categorised into three types, namely, pumps,tanks and valves. The survival signature is given in Table7.3. The survival functionof the water system is then calculated analytically as shown below:P(TS > t) =3∑l1=03∑l2=07∑l1=0(l1, l2, l3)(3l1)[1 − e−λ1t ]3−l1 [e−λ1t]l1 ×(3l2) [1 − e−λ2t]3−l2 [e−λ2t]l2 ×(7l3) [1 − e−λ3t]7−l3 [e−λ3t]l3 (7.4)The resulting survival functions without repair and with repairable components areshown in Fig. 7.8.As shown in Fig. 7.8, the results of the simulation method are in agreement withthe analytical solution for both repairable and non-repairable components. The pro-posed simulation method is applicable to any distribution type, intervals or evenprobability boxes. It not only separates the system structure from its component fail-ure time distributions, but also doesn’t require the i id assumption between differentcomponent types, as illustrated in [14].110 H. George-Williams et al.Fig. 7.8 Survival function without repairable (left panel) and with repairable components (rightpanel) computed using 10000 samples and verified by the analytical solutions7.5 Final RemarksSystem topological complexity, component interdependencies, multi-state compo-nent attributes and complex maintenance strategies inhibit the application of sim-ple reliability engineering reasoning to systems. For systems characterised by theseattributes, simulation-based approaches allow the realistic analysis of their reliability,despite the relatively higher computational costs of these approaches. This, however,is not a problem, with recent advancement in computing.The load flow simulation approach is an intuitive simulation framework that isapplicable to binary and multi-state systems of any topology. It does not require theprior definition of the structure function, minimal cut sets or the minimal path sets ofthe system. Instead, it employs a linear programming algorithm and the principles offlow conservation to compute the flow through the system. Thus, it can model flowlosses and implement reconfiguration requirements relatively easily. It can modelall forms of interdependencies in realistic systems, using intuitive representations.These attributes render the framework intuitive and generally applicable.While the load flow simulation approach is optimised for multi-state systems, itmay not be the best for binary-state systems with identical components. Since thesurvival signature is a function of the system topology only, it can be calculated onlyonce and reused in multiple reliability analyses. This feature reduces the reliabilityevaluation of the system to the analysis of the failure probabilities of its components,which is computationally cheap. Efficient simulation methods based on system sur-vival signature allow the reliability analysis of complex systems without resorting tosimplifications or approximations.The load flow and survival signature simulation approaches are not alternative toeach other; instead, they can be coupled to take advantage of their unique features,especially for systems with multiple outputs and potentially, multi-state systems.7 Simulation Methods for the Analysis of Complex Systems 111The algorithms and examples presented are available at: https://github.com/cossan-working-group/SystemReliabilityBookChapter.References1. M. Abd-El-Barr. Design and analysis of reliable and fault-tolerant computer systems. Hack-ensack, NJ, USA: World Scientific Publishing Co., 2006.2. Ahmad W. Al-Dabbagh and Lixuan Lu. Reliability modeling of networked control sys-tems using dynamic flowgraph methodology. Reliability Engineering and System Safety,95(11):1202 – 1209, 2010.3. L.J.M Aslett, F.P. Coolen, and S.P. Wilson. Bayesian inference for reliability of systems andnetworks using the survival signature. Risk Analysis, 2014.4. Michael Beer and Edoardo Patelli. Editorial: Engineering analysis with vague and impreciseinformation. Structural Safety, 52:143, 2015.5. Sergey V. Buldyrev, Roni Parshani, Gerald Paul, H. Eugene Stanley, and Shlomo Havlin.Catastrophic cascade of failures in interdependent networks. Nature, 464(7291):1025–1028,April 2010.6. Sungil Byun, Inseok Yang,MooGeun Song, and Dongik Lee. Reliability evaluation of steeringsystem using dynamic fault tree. In IEEE Intelligent Vehicles Symposium, pages 1416 – 1420,2013.7. F.P. Coolen and T. Coolen-Maturi. Generalizing the signature to systems with multiple typesof components. Complex Systems and Dependability, pages 115 – 30, 2012.8. Frank P.A. Coolen and Tahani Coolen-Maturi. Predictive inference for system reliability aftercommon-cause component failures. Reliability Engineering & System Safety, 135:27 – 33,2015.9. Alvarez DA. Infinite random sets and applications in uncertainty analysis. Arbeitsbereich fürTechnischeMathematik am Institut für Grundlagen der Bauingenieurwissenschaften. Leopold-Franzens-Universität Innsbruck, 2007.10. Salvatore Distefano and Antonio Puliafito. Reliability and availability analysis of dependent-dynamic systems with drbds. Reliability Engineering and System Safety, 94(9):1381–1393,2009.11. S. Eryilmaz. Review of recent advances in reliability of consecutive k-out-of-n and relatedsystems. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk andReliability, 224 (3):225 - 237, 2010.12. Samaniego FJ. System signatures and their applications in engineering reliability. SpringerScience & Business Media, 110.13. M.Beer F. P.CoolenG. Feng, E. Patelli. Imprecise system reliability and component importancebased on survival signature. Reliability Engineering & System Safety, 150:116 - 125, 2016.14. Hindolo George-Williams, Geng Feng, Frank Coolen, Michael Beer, and Edoardo Patelli.Extending The Survival Signature ParadigmToComplex SystemsWithNon-repairable Depen-dent Failures. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Riskand Reliability, 233(4):505–519, August 2019.15. Hindolo George-Williams, Geng Feng, Frank PA Coolen, Michael Beer, and Edoardo Patelli.Extending the survival signature paradigm to complex systems with non-repairable dependentfailures. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk andReliability, 233(4):505–519, 2019.112 H. George-Williams et al.16. Hindolo George-Williams, Min Lee, and Edoardo Patelli. Probabilistic risk assessment ofstation blackouts in nuclear power plants. IEEE Transactions on Reliability, 67(2):494–512,June 2018.17. Hindolo George-Williams and Edoardo Patelli. A hybrid load flow and event driven simulationapproach to multi-state system reliability evaluation. Reliability Engineering & System Safety,152:351–367, August 2016.18. Hindolo George-Williams and Edoardo Patelli. Efficient Availability Assessment of Reconfig-urable Multi-State Systems with Interdependencies. Reliability Engineering & System Safety,165:431–444, September 2017.19. Hindolo George-Williams and Edoardo Patelli. Maintenance Strategy Optimization forComplex Power Systems Susceptible to Maintenance Delays and Operational Dynamics.66(4):1309–1330, 2017.20. M. Beer J. Behrensdorf, M. Broggi. Imprecise reliability analysis of complex interconnectednetworks. Safety and Reliability - Safe Societies in a Changing World - Haugen et al. (Eds),Taylor & Francis Group, London, ISBN 978-0-8153-8682-7 year = 2018.21. G. Levitin and A. Lisnianski. Multi-state System Reliability Analysis and Optimization, in:Handbook of Reliability Engineering, chapter 4, pages 61–90. Springer, 2003.22. G. Levitin, L. Xing, H. Ben-Haim, and Y. Dai. Multi-state systems with selective propagatedfailures and imperfect individual and group protections. Reliability Engineering & SystemSafety, 96(12):1657–1666, 12 2011.23. Jing-An Li, YueWu, KinKeung Lai, andKe Liu. Reliability estimation and prediction ofmulti-state components and coherent systems. Reliability Engineering & System Safety, 88(1):93 –98, 2005.24. A. Lisnianski and Y. Ding. Inverse lz-transform for a discrete-state continuous-time markovprocess and its application to multi-state system reliability. In Applied Reliability Engineeringand Risk Analysis. Wiley, 2014.25. Anatoly Lisnianski, Ilia Frenkel, and Yi Ding. Multi-State System Reliability Analysis andOptimization for Engineers and Industrial Managers. Springer-Verlag London Limited, 2010.26. ManishMalhotra and Kishor S. Trivedi. Dependability modeling using petri-nets. IEEE Trans-actions on Reliability, 44(3):428 – 440, 1995.27. A Mosleh, D M Rasmuson, and F M Marshall. Guidelines on modeling Commom-CauseFailures in probabilistic risk assessment. Technical Report NUREG/CR-5485, U.S. NuclearRegulatory Commission, 1998.28. AndrewO’Connor andAliMosleh.Ageneral cause basedmethodology for analysis of commoncause and dependent failures in system risk and reliability assessments.Reliability Engineering& System Safety, 145:341 – 350, 2016.29. Min Ouyang. Review on modeling and simulation of interdependent critical infrastructuresystems. Reliability Engineering & System Safety, 121:43 – 60, 2014.30. Edoardo Patelli. Handbook of Uncertainty Quantification, chapter COSSAN: A Multidisci-plinary Software Suite for Uncertainty Quantification and Risk Management, pages 1–69.Springer International Publishing, 2017.31. Edoardo Patelli, Hindolo George-Williams, Jonathan Sadeghi, Roberto Rocchetta, MatteoBroggi, and Marco de Angelis. OpenCossan 2.0: An efficient computational toolbox for risk,reliability and resilience analysis. In André T. Beck, Gilberto F. M. de Souza, and Marcelo A.Trindade, editors, Proceedings of the Joint ICVRAM ISUMA UNCERTAINTIES Conference,2018.32. Broggi M de Angelis M. Patelli E, Alvarez DA. Uncertainty management in multidisciplinarydesign of critical safety systems. J Aerosp Inf Syst, 12:140-69, https://doi.org/10.2514/1.I010273,2015.33. Jose E. Ramirez-Marquez andDavidW. Coit. Optimization of system reliability in the presenceof common cause failures. Reliability Engineering & System Safety, 92(10):1421 – 1434, 2007.34. S. Reed. An efficient algorithm for exact computation of system and survival signatures usingbinary decision diagrams. Reliability Engineering & System Safety, 165:257 - 267, https://doi.org/10.1016/j.ress.2017.03.036.,2017.7 Simulation Methods for the Analysis of Complex Systems 11335. DanM. Shalev and Joseph Tiran. Condition-based fault tree analysis (cbfta): A newmethod forimproved fault tree analysis (fta), reliability and safety calculations. Reliability Engineeringand System Safety, 92:1231 – 1241, 2007.36. Masoud Taheriyoun and Saber Moradinejad. Reliability analysis of a wastewater treatmentplant using fault tree analysis and monte carlo simulation. Environmental Monitoring andAssessment, 187(1), 2015.37. Matthias CM Troffaes, Gero Walter, and Dana Kelly. A robust Bayesian approach to model-ing epistemic uncertainty in common-cause failure models. Reliability Engineering & SystemSafety, 125:13–21, 2014.38. M. Veeraraghavan and K.S. Trivedi. A combinatorial algorithm for performance and reliabilityanalysis using multistate models. IEEE Transactions on Computers, 43(2):229–234, Feb 1994.39. Liudong Xing and Yuanshun Dai. A new decision-diagram-based method for efficient analysison multistate systems. Dependable and Secure Computing, IEEE Transactions on, 6(3):161–174, July 2009.40. Wei-Chang Yeh. An improved sum-of-disjoint-products technique for the symbolic networkreliability analysis with known minimal paths. Reliability Engineering & System Safety,92(2):260 – 268, 2007.41. Wei-Chang Yeh. A fast algorithm for quickest path reliability evaluations in multi-state flownetworks. Reliability, IEEE Transactions on, 64(4):1175–1184, Dec 2015.42. Wei-Chang Yeh, Yi-Cheng Lin, Y.Y. Chung, andMingchang Chih. A particle swarm optimiza-tion approach based on monte carlo simulation for solving the complex network reliabilityproblem. Reliability, IEEE Transactions on, 59(1):212–221, March 2010.43. Enrico Zio, Piero Baraldi, and Edoardo Patelli. Assessment of the availability of an offshoreinstallation by monte carlo simulation. International Journal of Pressure Vessels and Piping,83(4):312 – 320, 2006.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 8Overview of Stochastic Model Updatingin Aerospace Application UnderUncertainty TreatmentSifeng Bi and Michael BeerAbstract This chapter presents the technique route of model updating in the pres-ence of imprecise probabilities. The emphasis is put on the inevitable uncertainties,in both numerical simulations and experimental measurements, leading the updat-ing methodology to be significantly extended from deterministic sense to stochas-tic sense. This extension requires that the model parameters are not regarded asunknown-but-fixed values, but random variables with uncertain distributions, i.e. theimprecise probabilities. The final objective of stochastic model updating is no longera single model prediction with maximal fidelity to a single experiment, but ratherthe calibrated distribution coefficients allowing the model predictions to fit withthe experimental measurements in a probabilistic point of view. The involvementof uncertainty within a Bayesian updating framework is achieved by developing anovel uncertainty quantification metric, i.e. the Bhattacharyya distance, instead ofthe typical Euclidian distance. The overall approach is demonstrated by solving themodel updating sub-problem of the NASA uncertainty quantification challenge. Thedemonstration provides a clear comparison between performances of the Euclidiandistance and the Bhattacharyya distance, and thus promotes a better understand-ing of the principle of stochastic model updating, as no longer to determine theunknown-but-fixed parameters, but rather to reduce the uncertainty bounds of themodel prediction and meanwhile to guarantee the existing experimental data to bestill enveloped within the updated uncertainty space.S. Bi (B)School of Aerospace Engineering, Beijing Institute of Technology, Beijing, Chinae-mail: sifeng.bi@bit.edu.cnM. BeerInstitute for Risk and Reliability, Leibniz Universität Hannover, Hanover, Germanye-mail: beer@irz.uni-hannover.de© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_8115116 S. Bi and M. Beer8.1 IntroductionComputational models of large-scale structural systems with acceptable precision,robustness and efficiency are critical, especially for applicationswhere a large amountof experimental data is hard to be obtained such as the aerospace engineering. Modelupdating has been developed as a typical technique to reduce the discrepancy betweenthe numerical simulations and the experimental measurements. Recently, it is a ten-dency to consider the inevitable uncertainties involved in both simulations and exper-iments. A better understanding of the discrepancy between them, in the backgroundof uncertainty, would achieve a better outcome of model updating.In structural analysis, the source of uncertainties, i.e. the reason of the discrep-ancy between simulations and measurements, can be classified into the followingcategories:• Parameter uncertainty: Imprecisely known input parameters of the numericalmodel, such asmaterial properties of novel composites, geometry sizes of complexcomponents and random boundary conditions;• Modelling uncertainty: Unavoidable simplifications and idealisations, such as lin-earised representations of nonlinear behaviours and frictionless joint approxima-tions;• Experimental uncertainty: Hard-to-control random effects, such as environmentnoise, measurement system errors and human subjective judgments.The above uncertainties motivate the trend to extend model updating from thedeterministic sense to the stochastic sense. The stochastic updating techniques drawmassive attention in the literature, in which the majority is based on the framework ofimprecise probability [4]. Considering the very typical categorisation of uncertain-ties, the term “imprecise probability” can be understood separately as “probability”corresponding to the aleatory uncertainty, and “imprecise” corresponding to the epis-temic uncertainty. The epistemic uncertainty is caused by the lack of knowledge. Asa better understanding of the problem is achieved, this part of uncertainty can bereduced by model updating. The aleatory uncertainty represents the natural random-ness of the system, such as the random wind load on launch vehicles, manufactureand measurement system errors. This part of uncertainty is irreducible, however, anappropriate quantification of the aleatory uncertainty is still required in stochasticmodel updating.The involvement of aleatory and/or epistemic uncertainties provides a clear logicfor the categorisation of model input parameters:I Parameters without any aleatory or epistemic uncertainty, appearing as explicitlydetermined constants;II Parameters with only epistemic uncertainty, which are still constants, however,with an undetermined position in an interval;III Parameters with only aleatory uncertainty, which are no longer constants, butpresented as random variables with exactly known distribution characteristics;8 Stochastic Model Updating with Uncertainty Treatment 117IV Parameters with both aleatory and epistemic uncertainties, which are randomvariables with undetermined distribution characteristics.The above parameters with various uncertainty characteristics lead the modelpredictions into the fourth category, i.e. the outputs with imprecise probabilities.Consequently, the objective of stochastic model updating is no longer a single modelpredictionwithmaximal fidelity to a single experiment, but rather aminimised uncer-tainty space of the outputs, whose bounds should still encompass the existing exper-imental data. In order to achieve this objective, a novel Uncertainty Quantification(UQ) metric based on the Bhattacharyya distance is introduced in this chapter. TheUQmetric refers to an explicit value quantifying the discrepancy between simulationsand measurements. Clearly, this metric is expected to be as comprehensive as pos-sible to capture all sources of uncertainty information simultaneously. Furthermore,the overall updating procedure is committed to being simple enough, by employinga uniform framework applicable to both of the classical Euclidian distance and thenovel Bhattacharyya distance. Within this uniform framework, comparison betweenthese two distance-based metrics can be performed conveniently, and therefore a bet-ter understanding of the difference between the deterministic and stochastic updatingis achieved with significantly reduced calculation cost.The following parts of this chapter are organised as follows. Section 8.2 givesan overview of the state of the art of deterministic and stochastic model updatingwhere key literature is provided. Section 8.3 presents the overall technique route withdescription of the key aspects, which can be helpful to generate a preliminary figureof the overall model updating campaign. As the emphasis, the parameter calibra-tion procedure with uncertainty treatment is explained in Sect. 8.4, where the Bhat-tacharyya distance-based UQ metric is presented along with the Bayesian updatingframework. The NASA UQ challenge problem is solved by the proposed approachand some interesting results are compared and analysed in Sect. 8.5. Section 8.6presents the conclusions and prospects.8.2 Overview of the State of the Art: Deterministicor Stochastic?Although this chapter focuses on the stochastic model updating, the deterministicupdating is still the footstone of its stochastic extension. A comprehensive review ofthe deterministicmodel updating techniques can be found fromRef. [14]. The readersare also suggested to refer to the fundamental book by Friswell and Mottershead[15] on this subject covering key aspects including model preparation, vibration dataacquisition, sensitivity analysis, error localisation, parameter calibration, etc.Among the plentiful techniques for deterministic parameter calibration, thesensitivity-based method is one of the most popular approaches based on the lin-earisation of the generally nonlinear relation between the model inputs and outputs.Mottershead et al. [20] provide tutorial literature for the sensitivity-based updat-118 S. Bi and M. Beering procedure of finite element models with both demonstrative and industry-scaleexamples. However, the sensitivity-based method is only valid for typical outputs,e.g. natural frequencies, modal shapes and displacement, of modal analysis or staticanalysis. Other more complex applications such as strong nonlinear dynamics ortransient analysis lead the analytically solved sensitivity to be unpractical. Conse-quently, the random sampling method, more specifically the Monte Carlo method,attracts more and more interest by providing a direct connection between the modelparameters and any output features via multiple deterministic model evaluations.The rapid development of computational hardware makes it possible for large-sizesamples, fromwhich the statistical information of the inputs/outputs can be preciselyestimated [17, 18]. SuchMonte Carlo-basedmethods have been successfully appliedin large-scale structures, see e.g. Refs. [10, 16].The widely used Monte Carlo methods obviously benefit to the research ofstochastic model updating, meanwhile, its conjunction with the Bayesian theoryfurther promotes this topic. Beck and Katafygiotis [3] proposed the fundamentalframework of Bayesian updating, whichwas further developed via theMarkov ChainMonte Carlo (MCMC) algorithm by Beck and Au [2]. The Bayesian updating frame-work in conjunction with the MCMC algorithm possesses the advantage to capturethe uncertainty information presented by rare experimental data. This approach hasbeen developed as a standard solution of stochastic model updating for differentapplications such as uncertainty characterisation [22].Besides the Bayesian interface, other imprecise probability techniques also haveconsiderable potential to be applied in stochastic model updating, such as intervalprobabilities [13], evidence theory [21], info-gap theory [5] and fuzzy probabilities[19]. For the background of imprecise probability, the comprehensive review by Beeret al. [4] is suggested for an overall understanding of this topic.8.3 Overall Technique Route of Stochastic Model UpdatingThe implement of stochastic model updating requires a complete theoretic systemincluding various key steps from the originally developed model, with a series oftechniques to define the features, to select and calibrate the parameters, to locateand reduce the modelling errors, and finally to validate the model with indepen-dent measurements. Special treatment of uncertainty propagation and quantificationpromotes the extension of model updating from deterministic domain to stochasticdomain. This extension is specifically implemented to key steps such as parametercalibration, model adjustment and validation. Overview of all related steps in modelupdating is provided in the following subsections.8 Stochastic Model Updating with Uncertainty Treatment 1198.3.1 Feature ExtractionA numerical model cannot be universally feasible for all scenarios with differentoutput features. Here, the feature is defined as the quantity that the engineer wants topredict with the model and it is also dependent on the capacity of the practical exper-imental set-up. Clearly, different features lead to different sensitivities of parameters,and require different strategies for uncertainty quantification. And therefore the firststep of model updating is to define a suitable feature according to the existing exper-imental setup and the practical application. The most typical features in structuralupdating are themodal quantities, e.g. the natural frequencies andmodal shapes.Mul-tiple orders of frequencies constitute a vector and the absolute mean error betweentwo vectors can be easily utilised to quantify the discrepancy between the simulatedand measured data. For modal shapes, the Modal Assurance Criterion (MAC) [1]is the most popular tool to quantify the correlation between two sets of eigenvec-tors. The continuous quantities are also commonly utilised as features, such as thedisplacement response in time domain and Frequency Response Functions (FRFs)in frequency domain. Classical techniques to evaluate the difference between twocomplex and continuous quantities are the Signature Assurance Criterion (SAC) andCross Signature scale Factor (CSF). Reference [8] provides an integrated applicationof SAC and CSF for a comprehensive comparison between two FRFs.8.3.2 Parameter SelectionA sophisticated model of a large structure system always contains a massive numberof parameters, which lead to a huge calculation burden and even failure of the updat-ing procedure. Parameter selection is therefore a key step to select or filter parametersaccording to their significances to the features defined in the first step. The core tech-nique for parameter selection is the sensitivity analysis focusing on a quantitativemeasurement of the parameter significance. The classical sensitivity analysis tech-nique is the Sobol variance-based method [25]. For a comprehensive knowledge ofthe global sensitivity analysis inspired by Sobol’s method, the well-written book bySaltelli et al. [24] is suggested to the readers. Another extension of Sobol methodincludes, e.g. the Analysis of Variance (ANOVA) based on the hypothesis testing inprobabilistic theory. Reference [7] proposes an integrated application of ANOVAandDesign of Experiment (DoE), which provides a significant coefficient matrix con-taining the complete sensitivity information of amultiple parameter–multiple featuresystem. When uncertainties are involved, the sensitivity analysis requires extensionfrom the deterministic procedure to the stochastic procedure. The sensitivity of acertain parameter, in the background of uncertainty treatment, can be represented asthe degree of how much the uncertainty space of the features is reduced, when theepistemic uncertainty space of the parameter is completely reduced. This requiresadditional techniques for uncertainty propagation and quantification, which will be120 S. Bi and M. Beeraddressed in the following parameter calibration step. A comprehensive literaturereview on the subject of sensitivity analysis can be found in Ref. [26].8.3.3 Surrogate ModellingThe employment of surrogate models, known as mate-models in some literature,becomes increasingly important along with the sampling-based updating method-ologies where a large number of model evaluations are generally required. A surro-gate model is a fast-running script between inputs and outputs, which can replacethe time-consuming model, e.g. the large finite element model, in the updating pro-cedure. An original input/output sample, i.e. a training sample, generated from thecomplex model is required to train the surrogate model. Since the surrogate model isproposed to handle the conflict between efficiency and precision, the training sampleis expected to be as small as possible, while the precision of the surrogate modelaccording to the original model should be high enough. The typical types of surro-gate models include the polynomial function, radial basis function, support vectormachine, Kriging function neural network, etc. The selection of a suitable surrogatemodel type is determined by various aspects such as its efficiency, generality, andnonlinearity. Another technique, in conjunction with the surrogate modelling, is theDoE with the purpose to efficiently and uniformly configure a spatial distribution ofthe training sample within the whole parameter space. A comprehensive review ofthe existing techniques of surrogate modelling and DoE can be found in Ref. [27].8.3.4 Test Analysis Correlation: Uncertainty QuantificationMetricsTest Analysis Correlation (TAC) is the core step of the overall updating procedure,not only because it significantly influences the updating outcome but also becauseit is the part mostly extended by the uncertainty treatment. TAC refers to the pro-cess to quantitatively measure the agreement (or lack thereof) between test mea-surements and analytical simulations, taking uncertainties into account. It thereforerequires a comprehensive metric which is capable of capturing multiple uncertaintysources simultaneously. This chapter proposes UQ metrics under various distanceconcepts. The Euclidian distance, i.e. the absolute geometry distance between twosingle points, is probably the most common metric used in deterministic updatingapproaches. However, it becomes insufficient for stochastic updating where mul-tiple simulations and multiple tests are presented. The Mahalanobis distance is aweighted distance considering the covariance between two datasets. And the alter-native Bhattacharyya distance is a statistical distance measuring the overlap betweentwo random distributions. A comprehensive comparison among the three distances8 Stochastic Model Updating with Uncertainty Treatment 121in model updating and validation can be found in Ref. [9], where the Bhattacharyyadistance is found to be more comprehensive to capture more sources of uncertaintyinformation. In the example section, the Bhattacharyya distance-based UQ metricis applied within a Bayesian updating framework, and the result is compared withthe one using the typical Euclidian distance. This work does not address the Maha-lanobis distance, since it has been found to be infeasible for parameter calibration,although it contains the covariance information among multiple outputs. Neverthe-less, the Mahalanobis distance has the potential to be utilised in model validation asdemonstrated in Ref. [9].8.3.5 Model Adjustment and ValidationAfter the model parameters are calibrated, the model still needs to pass the validationprocedure before it can be utilised in a practical application. The validation procedurecontains a series of criteria with increasing requirements: (1) The updated modelshould predict the existing measurements; (2) The updated model should predictanother set of measurement data which is different from the ones used for updating;(3) The updated model should predict any modification of the physical system bymaking the samemodification on themodel; (4) The updatedmodel, when utilised asa component of a whole system, should improve the prediction of the whole systemmodel. In the background of uncertainty treatment, the fit between the predictionand the measurement should be assessed by not only the precision but also thestochastic characteristics, e.g. probabilities, intervals and probability boxes. Modeladjustment is a procedure to deal with the modelling uncertainty. When the updatedmodel fails to fulfil the validation criteria, or some updated parameters are found tobe unphysical, e.g. a minus density value, the modelling uncertainty is too severeto be compensated by calibrating the parameters. The model can be adjusted byincreasing the resolution, changing the element type and adding a more detailedgeometry description, etc. Another round of parameter calibration is performed forthe adjusted model until the validation criterion is found to be fulfilled.8.4 Uncertainty Treatment in Parameter Calibration8.4.1 The Bayesian Updating FrameworkThe typical Bayesian model updating methodology is based on the following Bayes’equation:P(θ |Xexp) = PL(Xexp|θ)P(θ)P(Xexp)(8.1)122 S. Bi and M. Beerwith the key elements described as follows:• P(θ) is the prior distribution of the parameters, representing the prior knowledgebefore model updating;• P(θ |Xexp) is the posterior distribution of the parameters conditional to the existingmeasurement, i.e. P(θ |Xexp) is the outcome of Bayesian updating;• P(Xexp) is the so-called “normalisation factor” guaranteeing the integration of theposterior distribution P(θ |Xexp) equal to one;• PL(Xexp|θ) is the likelihood function defined as the probability of the existingmeasurements conditional to an instance of the parameters.The likelihood represents the probability of the measurement data under eachinstance of the updating parameters θ . And thus the objective ofmodel updating in theBayesian background is expressed as: to find the specific instance of the parametersallowing the experimental measurement to possess the largest probability, in otherwords, allowing the likelihood PL(Xexp|θ) reach the maximum. See Chap. 1 foradditional background on Bayesian inference.However, one of the difficulties in Bayesian updating is relative to the normal-isation fact P(Xexp). The direct integration of the posterior distribution over thewhole parameter space is quite difficult in practical application, especially when thenumber of parameters is large and the distribution format is complex, leading thedirect evaluation of the normalisation factor impractical. The well-known MCMCalgorithm is popular to solve this difficulty by replacing Eq. (8.1) withP(θ |Xexp) ≈ PL(Xexp|θ)β P(θ) (8.2)where β is the weighting coefficient fallen within the interval [0, 1]. When β equalsto zero, the right part of Eq. (8.2) is the prior distribution; when β equals to one, theright part of Eq. (8.2) converges to the posterior distribution. In the j-th iteration of theMCMCalgorithm, random samples are generated from the intermediate distributionswith weighting coefficient β j ∈ [1, 0]. In the (j + 1)-th iteration, parameter pointswhich lead to higher likelihood are selected from the random samples in the j-thiteration. New Markov chains are generated using the selected parameter points,and thus β j+1 is updated. The intermediate distribution converges to the posteriordistribution when β j = 1. The MCMC algorithm has been developed as a standardtool to stepwise generate samples from very complex target distributions. For thedetailed description of the MCMC algorithm, Refs. [2, 11] are suggested to thereaders.More applications of this algorithm can be found in the fields from stochasticmodel updating [6, 22] to structural healthmonitoring [23]. SeeChap. 2 for additionalbackground on Monte Carlo methods.8 Stochastic Model Updating with Uncertainty Treatment 1238.4.2 A Novel Uncertainty Quantification MetricThe MCMC algorithm makes it possible to avoid the evaluation of the normalisationfact in Bayes’ equation. However, the evaluation of the likelihood is inevitable andit becomes even more critical along with the uncertainty treatment. In the presenceof multiple sets of measurement data, the theoretical definition of the likelihood isPL(Xexp|θ) =Nexp∏k=1P(xk |θ) (8.3)where Nexp is the number of existing experiments. This equation requires accurateknowledge of the distribution of each measurement data point P(xk |θ). An accurateestimation of the distribution requires a large number of samples, which means alarge number of model evaluations.Clearly, a direct evaluation of Eq. (8.3) leads to a huge calculation cost. This iswhy the Approximate Bayesian Computation (ABC) becomes increasingly popularin Bayesian applications. Considering the principle of the likelihood function, itis natural to propose an approximate function to replace Eq. (8.3), as long as thisapproximate function still contains the information of the existing measurement dataand an instance of updating parameters. In this work, an approximate likelihoodbased on the Gaussian function is proposed asPL(Xexp|θ) ≈ 1σ√2πexp{−d(Xexp,Xsim)22σ 2}(8.4)where d(Xexp,Xsim) is the distance between the experimental and simulated featuredata. Equation (8.4) serves as an elegant connection between the Bayesian updatingframework and the distance concepts. It provides a uniform framework for variousdistance concepts. The approximate likelihood is applicable for not only theEuclidiandistance but also for the Mahalanobis and Bhattacharyya distances, in a uniformupdating framework. As explained in Sect. 8.3.4, the Mahalanobis distance is notutilised in this work. The Euclidian distance is evaluated asdE (Xexp,Xsim) =[(Xexp − Xsim)(Xexp − Xsim)T]1/2(8.5)whereXexp andXsim are the experimental and simulated feature datamatrices, respec-tively; X denotes the mean vector of the matrix. Clearly, the Euclidian distance onlyhandles the mean of the data. The Bhattacharyya distance has the definition asdB(Xexp, Xsim) = − log[∫Y√Pexp (x) Psim (x)dx](8.6)124 S. Bi and M. Beerwhere P(x) is the PDF of the feature; X is the feature space, implying∫Xis the inte-gration performed over the whole feature variable space. More detailed informationabout the evaluation method of the Bhattacharyya distance can be found in Ref. [6].8.5 Example: The NASA UQ ChallengeThe NASA UQ challenge [12] has been developed as a benchmark problem foruncertainty treatment techniques in multidisciplinary engineering. This challengeproblem contains a series of subproblems such as uncertainty characterisation, sen-sitivity analysis, uncertainty propagation, etc. In this work, only the SubproblemA (Uncertainty Characterisation) is investigated since it is equivalent to the task ofmodel updating herein. Figure 8.1 illustrates the key components of this problem,which contains one output evaluated via a black-box model using five input param-eters. According to the categorisation strategy in Sect. 8.1, the five parameters areclassified into three categories, as shown in Table 8.1.Only the parameters involving epistemic uncertainties, i.e. Categories II and IVparameters, require to be updated in this context. The uncertainty characteristics,including distribution types and distribution coefficients, are predefined as listedin Table 8.1. The predefined intervals of the distribution coefficients represent theepistemic uncertainty of the parameters. An output sample with 50 data points isavailable in the problem, which is generated by assigning a set of “true” values of thedistribution coefficients from the predefined intervals. The objective of this problemis, based on the existing 50 output points, to reduce the epistemic uncertainty space ofthe parameters, i.e. to reduce the predefined intervals of the distribution coefficients.Although the number of the model parameter is five, there are totally eight updatingcoefficients controlling the epistemic uncertainty space of the parameters, as shownin the last a column of Table 8.1.To solve the problem, the Euclidian and Bhattacharyya distances are utilised toconstruct the approximate likelihood functions, respectively. The Bayesian updatingFig. 8.1 Key components of the Subproblem A in the NASA UQ challenge8 Stochastic Model Updating with Uncertainty Treatment 125Table 8.1 The uncertainty characteristics of the parametersParameters Categories UncertaintycharacteristicsUpdating coefficientsp1 IV Unimodal Beta,μ1 ∈[0.6, 0.8],σ 21 ∈[0.02, 0.04]θ1 = μ1, θ2 = σ 21p2 II Constant, p2 ∈ [0.0,1.0]θ3 = p2p3 III Uniform, μ3 = 0.5,σ 23 = 1/12No updating requiredp4, p5 IV Joint Gaussian,μi ∈[−5.0, 5.0],σ 2i ∈[0.0025, 4.0],ρ∈[-1.0, 1.0], i = 4, 5θ4 = μ4, θ5 = σ 24θ6 = μ5, θ7 = σ 25 ,θ8 = ρframework employs these two distances as UQ metrics and generates two indepen-dent sets of results. This treatment is intended to make a clear distinction betweenthe deterministic updating and the stochastic updating, and furthermore to revealthe merits and demerits of these two distances. In practical applications, however, acombined application of these two distances is suggested by first using the Euclidiandistance for the mean updating, and second using the Bhattacharyya distance forthe variance updating. This two-step strategy has been demonstrated as a success insolving this problem in Ref. [6].The posterior distributions of the eight updating coefficients using respectively theEuclidian and Bhattacharyya distances as the UQ metrics are illustrated in Fig. 8.2.In the figure, the objects (“Samples” or “PDFs”) with the suffix “_ED” denote theFig. 8.2 Posterior distributions updated using the Euclidian and Bhattacharyya distances126 S. Bi and M. Beerresults updated using the Euclidian distance metric; and the ones with the suffix“_BD” denote the results using the Bhattacharyya distance metric. Except the distri-bution of θ1, most of the posterior distributions with the Euclidian distance metric arestill close to uniform, implying the deterministic updating procedure employing theEuclidian distance is incapable of solving this problem. As an obvious comparison,the updating procedure employing theBhattacharyya distance performswell formostof the updating coefficients by providing very peaked posterior distributions, such asθ2, θ3, θ6, and θ7. The vertical line in each subfigure represents the true value of theupdating coefficient, which were used to generate the existing 50 output samples.By comparing the position of the vertical line with the peak of the posterior distribu-tion, the updating precision is assessed. For θ1, the peak of the distribution with theEuclidian distance is apart from the vertical line, while the peak of the distributionwith the Bhattacharyya distance is quite close to the vertical line. This means theBhattacharyya distance performs better than the Euclidian distance when updatingθ1. The same conclusion is also achieved for θ3, θ4, and θ7. However, for θ2, and θ5,although the distributions with the Bhattacharyya distance have clear poles, they donot converge to the vertical lines. The updating precision of these two coefficients ispossible to be further improved by the two-step strategy as described in Ref. [6]. Notethat, there are still two coefficients (θ5 and θ8), which cannot be calibrated by neitherthe Euclidian distance nor the Bhattacharyya distance. A potential explanation is thatthe sensitivity of these two coefficients to the output is extremely low, leading theinverse procedure impossible to locate the “true” value.The quantitative updating results are detailed in Table 8.2, where the true valuesof the updating coefficients and the updated ones are provided. The updated valuesin the last two columns of Table 8.2 are obtained by determining the exact posi-tion of the distribution peaks in Fig. 8.2. For the posterior distributions which arestill close to uniform distribution, the determining process is meaningless and thusomitted. This is why only two updated coefficients are provided in the column withthe Euclidian distance. Note that, although the true values of the coefficients arereleased in Table 8.2, they are not necessarily to be treated as the final target of thisupdating problem. In the background of uncertainty treatment, the objective of thisproblem, stipulated by the problem designer [12], is to reduce the epistemic uncer-tainty space of the parameters, while making sure that the existing output sample canstill be included in the output uncertainty space. Hence, it makes more sense to assesshow much the output uncertainty space has been reduced after the intervals of thecoefficients are reduced in the updating procedure. This is relative to the tasks uncer-tainty propagation and model validation, which are out of the scope of this examplefor parameter calibration. A complete validation procedure considering the uncer-tainty space of the output for this problem can be found in Ref. [6]. Nevertheless,the stochastic Bayesian updating framework employing the Bhattacharyya distancehas been demonstrated to be more comprehensive and feasible than the Euclidiandistance in solving this NASA UQ challenge problem.8 Stochastic Model Updating with Uncertainty Treatment 127Table 8.2 The updated results using the Euclidian and Bhattacharyya distancesUpdating coefficients Predefined intervals True valuesa Updated resultsEuclidian Bhattacharyyaθ1: μ1 [0.6, 0.8] 0.6364 0.6990 0.6194θ2: σ 21 [0.02, 0.04] 0.0356 – 0.0397θ3: p2 [0, 1] 1.0 – 0.9855θ4: μ4 [−5, 5] 4.0 – 3.9770θ5: σ 24 [0.0025, 4] 0.04 – –θ6: μ5 [−5 , 5] −1.5 – −4.4732θ7: σ 25 [0.0025, 4] 0.36 – 0.2818θ8: ρ [−1, 1] 0.5 −0.3471 0.1804aData available online at https://uqtools.larc.nasa.gov/nda-uq-challenge-problem-2014 [retrievedin 2017]8.6 Conclusions and ProspectsThe tendency of uncertainty analysis has been rendering the typical model updat-ing full of vitalities but also challenges. This chapter reviews the key techniquesand components of the overall model updating campaign. Main emphasis is put onthe involvement of uncertainty, which leads the transformation from the determinis-tic approach to the stochastic approach. The stochastic model updating is executedwithin the Bayesian model updating framework, where the Bhattacharyya distance isproposed as a novel UQmetric. The approximate likelihood is critical by providing auniform connection between the Bayesian framework and various types of distancemetrics. The example demonstrates that the Bhattacharyya distance is more compre-hensive and feasible than the Euclidian distance to calibrate distribution coefficientsof parameters with imprecise probabilities.The Bhattacharyya distance is designed as a universal tool of UQ, which can beconveniently embedded into a technique route similar as the deterministic approach,but provides stochastic outcomes by capturing more uncertainty information. Thetendency of uncertainty analysis will be further promoted by the novel UQ metricin not only the stochastic parameter calibration but also other procedures, e.g. thestochastic sensitivity analysis andmodel validation,whichwill establish the completescenario of the stochastic model updating.References1. R. Allemang. The modal assurance criterion: Twenty years of use and abuse. Sound and vibra-tion, 8:14–21, 2008.2. J. Beck and S.K. Au. Bayesian updating of structural models and reliability usingmarkov chainmonte carlo simulation. Journal of Engineering Mechanics, 128(4):380–391, 2002.128 S. Bi and M. Beer3. J. Beck and L. Katafygiotis. Updating models and their uncertainties. i: Bayesian statisticalframework. Journal of Engineering Mechanics, 124(4):455–461, 1998.4. M. Beer, S. Ferson, and V. Kreinovich. Imprecise probabilities in engineering analyses.Mechanical Systems and Signal Processing, 37(1-2):4–29, 2013.5. Y. Ben-Haim. Info-gap decision theory: Decisions under severe uncertainty. Elsevier, 2006.6. S. Bi, M. Broggi, and M. Beer. The role of the bhattacharyya distance in stochastic modelupdating. Mechanical Systems and Signal Processing, 117:437–452, 2019.7. S. Bi, Z. Deng, and Z. Chen. Stochastic validation of structural fe-models based on hierarchicalcluster analysis and advanced monte carlo simulation. Finite Elements in Analysis and Design,67:22–33, 2013.8. S. Bi, M. Ouisse, and E. Foltête. Probabilistic approach for damping identification consideringuncertainty in experimental modal analysis. AIAA Journal, 56(12):4953–4964, 2018.9. S. Bi, S. Prabhu, S. Cogan, and S. Atamturktur. Uncertainty quantification metrics with varyingstatistical information in model calibration and validation. AIAA Journal, 55(10):3570–3583,2017.10. A. Calvi. Uncertainty-based loads analysis for spacecraft: Finite element model validation anddynamic responses. Computers & Structures, 83(14):1103–1112, 2005.11. J. Ching and Y. Chen. Transitional markov chain monte carlo method for bayesian modelupdating, model class selection, and model averaging. Journal of Engineering Mechanics,133(7):816–832, 2007.12. L.Crespo, S.Kenny, andD.Giesy.Thenasa langleymultidisciplinary uncertainty quantificationchallenge. In 16th AIAA Non-Deterministic Approaches Conference.13. M. Faes, M. Broggi, E. Patelli, Y. Govers, J. Mottershead, M. Beer, and D. Moens. A multi-variate interval approach for inverse uncertainty quantification with limited experimental data.Mechanical Systems and Signal Processing, 118:534–548, 2019.14. M. Friswell and J. Mottershead. Model updating in structural dynamics: a survey. Journal ofSound and Vibration, 162(2):347–375, 1993.15. M. Friswell and J. Mottershead. Finite Element Model Updating in Structural Dynamics.Kluwer Academic Press, Dordrecht, Netherlands, 1995.16. B. Goller, M. Broggi, A. Calvi, and G.I. Schueller. A stochastic model updating technique forcomplex aerospace structures. Finite Elements in Analysis and Design, 47(7):739–752, 2011.17. Y. Govers and M. Link. Stochastic model updating - covariance matrix adjustment from uncer-tain experimental modal data. Mechanical Systems and Signal Processing, 24(3):696–706,2010.18. C. Mares, J. Mottershead, and M. Friswell. Stochastic model updating: Part 1 - theory andsimulated example.Mechanical Systems and Signal Processing, 20(7):1674–1695, 2006.19. B. Möller and M. Beer. Fuzzy Randomness: Uncertainty in Civil Engineering and Computa-tional Mechanics. Springer, 2004.20. J. Mottershead, M. Link, and M. Friswell. The sensitivity method in finite element modelupdating: A tutorial.Mechanical Systems and Signal Processing, 25:2275–2296, 2011.21. W. Oberkampf and J. Helton. Evidence theory for engineering applications. In Engineeringdesign reliability handbook, chapter 10, pages 197–226. CRC Press, 2004.22. E. Patelli, D. Alvarez, M. Broggi, and M. De Angelis. Uncertainty management in multidisci-plinary design of critical safety systems. Journal of Aerospace Information Systems, 12(1):140–169, 2015.23. R. Rocchetta, M. Broggi, Q. Huchet, and E. Patelli. On-line bayesian model updating forstructural health monitoring.Mechanical Systems and Signal Processing, 103:174–195, 2018.24. A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, andS. Tarantola. Global Sensitivity Analysis. The Primer. John Wiley & Sons, 2008.25. I. Sobol. Sensitivity estimates for nonlinear mathematical models.Mathematical Modeling andComputational experiment, 1(4):407–414, 1993.26. P. Wei, Z. Lu, and J. Song. Variable importance analysis: A comprehensive review. ReliabilityEngineering and System Safety, 142:399–432, 2015.8 Stochastic Model Updating with Uncertainty Treatment 12927. R. Yondo, E. Andrés, and E. Valero. A review on design of experiments and surrogate modelsin aircraft real-time and many-query aerodynamic analyses. Progress in Aerospace Sciences,96:23–61, 2018.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.Chapter 9Aerospace Flight Modelingand Experimental TestingOlivier ChazotAbstract Validation processes for aerospace flight modeling require to articulateuncertainty quantification methods with the experimental approach. On this note,the specific strategies for the reproduction of re-entry flow conditions in ground-based facilities are reviewed. It shows how it combines high-speed flow physics withthe hypersonic wind tunnel capabilities.9.1 IntroductionSpace missions are built upon massive technology knowledge and on the latestprogress in engineering. They fascinate as they represent for the public the mostadvanced knowledge as well as a typical dive into the unknown. For scientists andengineers, such missions are an occasion to push the scientific knowledge and toestablish better what we know to offer solid basis for further discoveries.In practice, space exploration leads to extreme challenges as it aims to investigateplanets, or asteroids, in the solar system and return samples for analysis across verysevere flight conditions. Such missions need to be designed using physical modelsand robust numerical methods. However, those tools used by aerospace engineersremain, for most of them, on the need for more research development for theirvalidation and consolidation to be able to plan successful and fruitful missions.Aerothermodynamic testing is one of the crucial points for the design of aerospacevehicles. At first place, it aims at establishing as much as knowledge possible oncritical flight phenomena. Ground-based facilities are operated to reproduce flight-relevant environment for the testing of the vehicle configuration and its ThermalProtection System (TPS) to allow for an accurate evaluation of their performances.Two types of facilities are classically used for the ground testing to support thepre-flight analysis. First, the required flow-field is reproduced for its analysis inO. Chazot (B)Aeronautics and Aerospace Department, von Karman Institute for Fluid Dynamics,Sint-Genesius-Rode, Belgiume-mail: chazot@vki.ac.be© The Author(s) 2022L. Aslett et al. (eds.), Uncertainty in Engineering,SpringerBriefs in Statistics,https://doi.org/10.1007/978-3-030-83640-5_9131132 O. Chazothigh-enthalpy facilities such as shock tunnels or expansion tubes. Then in a secondstep, the interaction between the dissociated gas and the vehicle’ surface is studiedto determine the thickness of Thermal Protection Material (TPM) from databasesbuilt in plasma wind tunnels. Accurate flight duplication is thus necessary in orderto properly address those stringent requirements without over-sizing the TPS.All the data produced are processed to define at best the modeling for all therequired physical phenomena. Very performant frameworks for such processing arethe Uncertainty Quantification (UQ) methods. Bayesian approach, in particular, areextremely powerful as they allow todetermine the required information froma limitedexperimental knowledge with a probabilistic treatment. However, such key results,for hypersonic flight, are not only due to powerful mathematical treatment but alsothanks to consistent experimental methodologies. This combination between mathe-matical and experimental approach is essential to generate useful knowledge on thephysical modeling to be determined. Then it is primordial to have a good understand-ing on how the ground testing are linked to the fundamental mathematical modelsused for simulating the high-speed flow physics. It serves to setup the best experi-mental environment to allow a fruitful processing of the acquired data.Therefore, this note intends to review the rationale of the ground testing method-ology for aerospace vehicle. It offers a synthesis on the high-speed ground testingunderlying the links to their scientific basis. It presents how similitude laws couldbe applied or need to be adapted as more and more physical phenomena have tobe considered for aerospace flights. It would help engineers and applied mathemati-cians for working together facing the challenges of high-speed flight investigationsfor aerospace development.9.2 Aerospace Flights and Planetary Re-entryAerospace flight can be considered for very different trajectories. It could follow aplanet orbit up to the limit of its escape velocity, but it could also correspond to aninterplanetary transfer with super-orbital velocities. Those situations involve a largevariety of kinetic energy and trajectories. Figure 9.1 illustrates those typical orbittrajectories, around Earth, with their corresponding velocities.Super-orbital atmospheric re-entry, also known as hyperbolic re-entry, is charac-terized by high velocities and is encountered when a probe is entering an atmospherefrom a hyperbolic orbit rather than from an elliptical orbit. It is the case for someinterplanetary probes or sample return missions. Typical velocities for probes enter-ing the Earth’s atmosphere from hyperbolic orbits scale from 11 to 14 km/s, whichcorrespond to specific enthalpies between 60 and 100 MJ/kg. This is considerablyhigher than the usual re-entry velocities from circular or elliptic orbits, e.g., 8.2 km/sfor the Space Shuttle. Up to now, the Stardust probe was the fastest artificial objectto perform a controlled re-entry in the Earth’s atmosphere, at 12.8 km/s.The environment of a high-speed re-entry flight is muchmore severe compare to aLowEarthOrbit (LEO) entrywith particularly highheat flux and important stagnation9 Aerospace Flight Modeling and Experimental Testing 133Fig. 9.1 Orbit trajectoriesand velocities. μ is the Earthgravitational parameter and ris the radius of the orbitpressure. In addition to the classical features of re-entry flows like non-equilibriumthermo-chemistry and complex gas–surface interactions some specific phenomenabecome important as shock layer radiative heating and ablation phenomena. Thosephysical process were not considered much at smaller velocity but they cannot beanymore neglected for super-orbital re-entry. In top of this, more complex physicalreality coupling phenomena appear in the flow between radiation and ablation thatlead to a very intricate situation. It is therefore not possible to extrapolate what hasbeen studied and learned for orbital re-entry to super-orbital conditions and specificground testing strategies are required.9.3 Similitude Approach for Hypersonic FlowsSimilitude in classical fluid dynamics establishes a correspondence between differentflows, based on the mathematical model representing these flows, without having tosolve the set of equations. This correspondence then can be used to relate two realphysical flow situations or to relate a family of solutions for the model.With two correlated applications:• The flow fields are similar (i.e., solution of the same set of equations) even if thedimensions and the temporal evolution of the phenomena are on different scales.• When applicable, the use of the similarity laws allows to replicate in wind tunnelsthe flow field occurring in flight around a re-entry vehicle.This second aspect of the similitude approach is mostly useful since it allows to studyon ground typical re-entry situations. Hypersonic flows are particularly interestingfor the application of similitude approach because it exists many different situa-tions in hypersonic regime which can be describes with a variety of models leadingto different similitude laws. On one hand, it gives opportunities to develop sim-plify solutions, but on the other hand, the physical nature of the high-speed flows,mostly due to the high-temperature effects, severely restrict exact similitude andimpose to study approximate similitude. This family of flows is essentially deter-mined by the mathematical model chosen to describe the flows and of which the134 O. Chazotflows are themselves solutions. It appears that the similitude strategy evolves as themodel integrates more and more physical aspects for the description of the flow. Inhypersonics going from inviscid regime to high-temperature effects, the similitudeapproach will be adapted to retain the most relevant flow parameters correspondingto each case. Similitude in hypersonic has been studied and discussed by differentauthors: in earlier time, by Hayes and Probstein for inviscid and viscous hypersonics[11, 12], by Freeman for dissociated gases [9] and in a more general form by Viviandfor CFD development [6], but it provides also very useful material for experimentalstudies. The most common approaches are briefly presented below, as well as theirlimitations.9.3.1 Inviscid HypersonicsAt high Mach number when considering a slender body (small thickness ratioτ << 1) at small angle of attack ( α << 1) in a perfect gas, the governing Eulerequations can be further simplified using the small disturbances theory. The similar-ity parameters are thenMτ, ατ , and γ . The slenderness ratio τ is defined as τ = d/ l,where d is the body’s radius and l its length. The parameter K = Mτ is called thehypersonic similarity parameter [11]. On these conditions, with a constant isentropicexponent γ , the results of the similitude may be expressed in dimensionless form.For family of affinely related bodies of thickness ratio, τ in two-dimensional flowsthe pressure coefficient Cp could be written:Cpτ 2= 2⎡⎣γ + 14+√(γ + 14)2+ 1K 2⎤⎦ = f (K , γ ) (9.1)The typical testing strategy in this inviscid framework could be represented as inthe figure (Fig. 9.2).That approach holds for inviscid hypersonic, over a wide range of K for slenderbodies as long as the Mach number in the flow is large enough and τ small. Never-theless, it does not correlate well as the body thickness is increased, curved shockand boundary layers start to be predominant in the flow. As most hypersonic vehiclesare blunt rather than slender for thermal considerations this hypersonic similarity haslimited applications.9.3.2 Viscous HypersonicsHypersonic flows present thick boundary layers that are also growing fast as Machnumber is increased. These boundary layers cannot be ignored generally in hyper-sonic problems as they determine the major feature of the flow physics. They also9 Aerospace Flight Modeling and Experimental Testing 135Fig. 9.2 Inviscid hypersonic similitudeneed to be taking into account for their interaction with the inviscid flow; ground-based facilities need to be designed toward those considerations to provide a way forstudying these effects in relevant hypersonic flight conditions. To this end, similitudeapproach is a precious guide to identify which combination of flow parameters needto be taken into account in experimental simulation. Viscous hypersonic similitudehas been presented first by Hayes and Probstein in a paper [12] where they reviewthe general features of viscous similitude at high speed. The inviscid flow need tobe represented, then the hypersonic similarity parameter K = Mτ and γ need to beinvariant. The interaction of the viscous part of the flow with the inviscid field willbe determined by the viscous-inviscid interaction parameter χ expressed asχ = M3∞√C∞√Rex ,∞ (9.2)with C∞ = μr T∞μ∞Tr and Rex,∞ = ρ∞Uxμ∞ .With these conditions the perfect gas model has to be assumed and the continuumhypothesis valid, i.e., the mean free path is at least one order of magnitude smallerthan the characteristic length of the flow. The viscous hypersonic similarity requiresreproducing the free-stream Reynolds (Re) and Mach number (M) and the tempera-ture ratio TwT∞ , where the subscriptw indicates the wall temperature and the∞ symbolin the subscript refers to the free-stream. If the gas mixture is not the same, the heatcapacity ratio γ also needs to be reproduced [13]. This eases facility development,as lowering the gas temperature lowers the speed of sound, and thereby increasesthe achievable free-stream Mach number. The condensation temperature of the test136 O. ChazotFig. 9.3 Schematic representation of cold hypersonic testing in ground facilitygas imposes the upper limit of what is achievable in a given facility. In a generalview it resorts a Mach–Reynold simulation. Most of the ground testing capabilitiesare designed following this principle working with scaled models (Fig. 9.3). Theiroperation envelope are commonly presented in a Ma–Re graph that indicates whichflight domain can be simulated (Fig. 9.4).1 One has to remark that the temperatureratio parameter TwT∞ is not often taking into account but it could be an important aspectto consider in ground testing as it could appear as a limitation in the study.9.3.3 High-Temperature HypersonicsWhen a real gas, as air, experiences a strong shock at high speed, it will increasetremendously its thermal energy. In such conditions, the molecular collisions areenergetic enough to cause dissociation and ionization and the gas to depart fromthe perfect gas model. For a blunt body at a velocity of 7 km/s, the temperatureimmediately after the shock is around 14,000 K, and around 8,000 K downstreamthe shock, where the flow may return to equilibrium. At such high temperatures,chemical effects have to be taken into account. For air at a pressure of 1 atm, vibra-tional excitation begins at 800 K, O2 begins to dissociate at 2,500 K and is fullydissociated for 4,000 K, point for which N2 begins to dissociate. At 9,000 K, N2is fully dissociated and ionization begins. One can easily understand that the flowdownstream the shock becomes a plasma: molecules are dissociated and atoms arepartially ionized. The parameters to be respected for a flow involving chemistry areV 2/2D, where V is the free-stream velocity and D the typical dissociation energyof the gas molecule considered, the Damköhler number Da, defined as Da = L/ lD ,1 Graphic extracted from AGARD AR 319.9 Aerospace Flight Modeling and Experimental Testing 137Fig. 9.4 Mach–Reynoldsmap for hypersonic facilities[20]where L is the characteristic length of the flow and lD the characteristic length asso-ciated to the dissociation reaction, and the temperature ratio Tw/T [13]. It shouldbe brought to the reader’s attention that the gas behind the shock is a chemicallyreacting mixture of perfect gases, and not a real gas as it is sometime incorrectlyreferred to in literature.The Damköhler number for the gas appears considering the mass conservationequation in a non-dimensional form as it is expressed in the equation below:∂ρi∂t+ .(ρi V + Ji ) =Nr∑r=1(ν′′ir − ν′ir ){Da f rNs∏j=1ρν′irj − DabwNs∏j=1ρν′′irj}(9.3)Da f r = k f rρ∞LU∞Dabw = kbwρ2∞LU∞(9.4)Chemical reactions take a very short but finite time to happen. Assuming thatthe flow is composed of a single species, the dissociation rate is proportional tothe density, while the recombination rate is proportional to the square of density.Hence, the characteristic lengths associated to the dissociation and recombinationreactions, respectively, scale as lD ∝ 1/ρ and lR ∝ 1/ρ2, where ρ is the flow den-sity [13]. The recombination length is usually larger than the dissociation length.Under certain conditions, at very high altitude, one can assume that lD ∼ L , and138 O. ChazotFig. 9.5 High-enthalpyfacility map [15]therefore Da = o(1), while lR is much larger. The flow is frozen; it is too fast forrecombination reactions to take place. In that case, the binary scaling parameter ρLmust be reproduced in order to obtain the correct Da [16]. One could note thatthe diffusion phenomena is also considered in this analysis, as the diffusion termin the equation scale with ρL [5]. With this approach, the high-enthalpy facilitiesare designed to reproduce the real flight velocities (V∞) and to allow for an opera-tion considering the ρL parameter as they can be represented in the graph of Fig. 9.5.This leads to some complications. Firstly, the same air mixture as in real flight iscommonly used, as the chemistry processes are too complicated to reproduce to useanother one. The typical dissociation energy D is therefore conserved, and the actualfree-streamvelocitymust be reproduced to duplicate the group V 2/2D. Secondly, therequired density to achieve in the wind tunnel becomes large in order to maintain theproper value of the binary scaling parameter for duplication of flight at lower altitude.Thirdly, the binary scaling approach, strictly speaking, is built upon the hypothesisof a single species mixture. That approach has therefore limited application for morecomplex mixtures such as air [7]. Finally, as altitude decreases, density increases andboth lD and lR become smaller. The binary scaling parameter does not hold anymore,as both ρL and ρ2L should be reproduced at the same time. The same holds whenflow velocity increases, and therefore also temperature. This prevents from usingthe similarity laws, and flow duplication can thus only be performed on a full-scalemodel, with the actual flow velocity.9 Aerospace Flight Modeling and Experimental Testing 1399.4 Duplication of Dissociated Boundary Layerwith Surface ReactionFollowing the development presented in the previous paragraph, it is observed thatin spite of all the possibilities offer by the ground testing, considering the high-temperatures effects in hypersonic, only the dissociation in the shock layer couldbe simulated to some extend on scaled models. If one wanted to further take intoaccount, the recombination in the gas he would have to consider full-scale models.Looking into more details, one could underline that all the important phenomena,concerning the heat transfer typically, are laying in the boundary layer. Figure 9.6gives an illustration of real flight situation in front of an Aerospace vehicle.The study could then be reduced to this confined layer in particular for the heat-transfer problems in hypersonic. In this situation, the full-scale environment willbe reproduced by duplicating all the characteristics of the boundary layer in theground testing facilities. This statement was already made in earlier publicationsfrom researchers working on dissociated boundary layers using shock tube facilities[19].The stagnation point is of particular interest for this duplication because the flowconditions, returning to zero velocity, are more easily reproduced in a laboratory.Before identifying the parameters that need to be retained for the testing conditionsit will be useful to give a physical description of reacting boundary layers, as theymanifest themselves along surfaces in hypersonic flows.Dissociated non-equilibriumboundary layers with surface reaction have already been presented extensively bymajor authors [2, 19] and the reader is encouraged to consult these references. TheFig. 9.6 Re-entry environment along stagnation line140 O. Chazotpurpose of this section is not to expose in details the theory specific to such boundarylayers but to recall their main characteristics and better understand how they could besimulated in ground-based facilities.Boundary layers are the locationwhere the diffu-sion phenomena are dominant, their relative importance are characterized by typicalnon-dimensional numbers: Prandtl (Pr), Lewis (Le), and Schmidt (Sc). Those willnot be commented here to rather focus on the chemical non-equilibrium feature of theflow. The chemical reactions taking place in the gas phase are called homogeneousreactions while those happening between the gas and the solid surface are calledheterogeneous reactions. The chemical non-equilibrium in the gas phase is charac-terized by the gas Damköhler number (Dag). It corresponds to the ratio betweenthe typical time of the flow, to cross the boundary layer, to a typical reaction timefor the gas chemistry: Dag = τ f /τc. When Dag → ∞ the boundary layer reachesa local thermodynamic equilibrium (LTE). On the contrary when Dag → 0 it leadsto a frozen boundary layer, where no chemistry happen. These different conditionshave significant consequences on the wall heat flux as it has been shown in reference[8]. The reactions at the wall are usually considered as first-order reactions, they arerepresented by a reaction rate (kw) for each dissociated specie. kw could be related toa recombination coefficient γi (or catalycity parameter), interpreted as a probabilityof recombination at the wall, by the Hertz–Knudsen formula:kwi = γi√k · Tw2π · Mi (9.5)The wall reaction rates are also characterized by a surface Damkohler number(Daw). It compares the time of diffusion for the species across the boundary layer tothe time of reaction at the wall: Daw = τDi f f /τReact . When Daw → ∞, it does notnecessarily imply that the reaction rate at the wall tend to infinity (kw → ∞), butsimply that the surface reaction are much more faster than the diffusion process. Itis said that the GSI phenomena are “limited by diffusion”. In the other extreme case,when Daw → 0 the diffusion ismuch faster than the reaction and theGSI phenomenaare “limited by reaction” or “reaction controlled”. From this description, it appearsthat the diffusion and reaction processes should be accurately reproduced. To thisextend, it could be understood that the dimension and the environment of the reactionboundary layer must be duplicated (Fig. 9.7). Two situations could be distinguishedfor the boundary layer to be duplicated in the laboratory: stagnation point region andoff-stagnation point when the boundary is developing along the surface.If one considers only the stagnation point region, it has been shown that a completeduplication of real flight condition is possible in ground facility, if the total enthalpy(He), the total pressure (Pe), and the velocity gradient (β = du/dx), of the flightconditions, can be matched locally on the test sample [4, 14]. In this case, the testingis realized in plasma wind tunnels (Arcjet or Plasmatron facilities) that are able toproduce dissociated flows for a long time base which is suitable for tests involvingaerothermochemistry. The theoretical frame for the testing with ICP wind tunnelshas been adapted by Russian scientists, in a methodology called Local Heat Transfer9 Aerospace Flight Modeling and Experimental Testing 141Fig. 9.7 Typical features of the dissociated stagnation point boundary layerSimulation (LHTS) [14] and an assessment of this methodology has been conductedat the VKI Plasmatron facility [1, 4]. The duplication of the flight condition atstagnationpoint is strictly reduced to theboundary layerwith its appropriate treatment[1]. In the case of a subsonic plasma facility, like the VKI Plasmatron, the testingconfiguration could be presented as in Fig. 9.8.The experimental assessment of the LHTS methodology has been conducted atVKI by Chazot et al. [4] and Barbante and Chazot [1]. The results are presentedon the Figs. 9.9 and 9.10. It could be seen that the boundary layers profiles fromthe hypersonic flow is very well duplicated with the subsonic plasma flow, whenmatching the characteristics parameters at the edge of the boundary layer.Fig. 9.8 TPS testing in plasma wind tunnel in LHTS conditions142 O. ChazotFig. 9.9 Temperatureprofiles comparison betweenflight and ground testingconditions [1]Fig. 9.10 Mass fractionprofiles comparison betweenflight and ground testingconditions [1]9.5 Considering Flow RadiationShock layer radiative heating appears around 9 km/s in Earth’s atmosphere and 7km/s inMars’ atmosphere [17]. It reaches 10%of the total heat flux for probes havinga diameter smaller than 1 m and entry velocities approaching 13 km/s in the Earth’satmosphere [21]. It is then a physical phenomena to be considered for super-orbitalre-entry.9 Aerospace Flight Modeling and Experimental Testing 143However on ground testing for scaledmodels flow radiation is not often addressedin the literature. It is usually mentioned to explain that this feature of the flow cannotbe properly reproduced considering similarity rules for amodel of reduced dimension[18]. But even if this conclusion is fully valid some more details could be added tobetter understand the issues with radiative heat-transfer in hypersonic facilities.The radiative flux at a certain point P is the integral of the radiative intensityin all directions and over the entire frequency spectrum. Considering a point S ona surface, the radiative flux reaching S from a point P of the surrounding is thedifference between the energy emitted and that absorbed integrated along the opticalpath from P to S. In the case of a scaled model in a high-enthalpy facility, respectingthe same condition for the gas and the flow velocity, it could be shown that the opticalthickness in the radiating shock layer scales with the product ρL , considering theabsorption coefficient remaining the same on the two situations.One is then brought back to the same approach as the binary scaling exposedbefore. In these conditions, the radiative heat-flux on the surface of a scaled modelcould be reproduced in a ground-based facility.However, if one considers the flow passing through the radiating environmentaround the model, it appears that the scaling does not hold any longer. The amountof energy E radiated by a control volume is proportional to the mass contained in thatvolume: E ∝ m = ρ.L3. The amount of flow ṁ ingested in the shock layer could beexpressed as ṁ ∝ ρ.U∞.L2. Therefore, when the flow-radiation coupling need to betaken into account the heat radiated per unit mass passing in a control volume scalesas: E/ṁ ∝ L . In conclusion, even if the radiative flux on the surface of a scaledmodel could be reproduced, the radiative heating of the flow around the same modelis not respected.This problem arises if the radiative heating is important enough to have an influ-ence on the rest of the flowfield. This coupling is quantified by the radiative coolingparameter , also referred to as Goulard number, defined as = 2 · Qrad1/2 · ρ∞ ·U 3∞(9.6)where Qrad is the radiative heating for an adiabatic flow, that is without radiativecooling, ρ∞ the free-stream density, and v the shock velocity. This parameter servesas an approximate measure of the coupling between radiation and the flow [10]. It isthe ratio of the amount of radiation generated by the shock, assumed to be twice theradiative heating of the surface, by the kinetic energy heat flux entering the shocklayer. If > 0.01, radiation is coupled to the rest of the flowfields and the effect ofimproper radiation scaling extends to the rest of the flowfield.Furthermore, there is a coupling between radiation and gas–surface interactionprocesses. Indeed, the use of ablative material is compulsory for high-speed re-entry. The ablation processes are very efficient to prevent the hot shock layer gasfrom reaching the wall and absorb part of the shock layer radiative heat flux. Anaccurate estimation of the absorbed radiation is complex since the thickness andthermochemical state of the ablation gas layer are difficult to predict up to now [17].144 O. ChazotOther surface phenomena such as catalycity and oxidation should be also taken intoaccount, it makes the problem even more intricate when one is aware that such amodels are not yet fully established, especially in high-enthalpy flows [3].9.6 Ground Testing Strategy for High-Speed Re-entryGround testing facilities commonly used to study re-entry flows and that are mainlyconcerned with high-enthalpy flows can be divided in two categories:• Impulse facilities, such as shock tubes and ballistic ranges, are only able to produceflows that last typically a fraction of a second. It is usually assumed long enoughto let the steady flow establish itself, but too short compared to the thermal inertiaof the material surface. They are thus mainly used to investigate the aerothermo-dynamic effects, gas kinetic, and radiation processes. In this category, expansiontubes are able to reach free-stream enthalpies characteristic of high-speed re-entry.• Plasma wind tunnels, such as inductively coupled plasma facilities or arc-jets, areable to operate for long test durations, in the order of minutes. However, they havenot been designed to reproduce the flow radiation.Similar flights conditions or direct flight duplication are possible to reproduce forsub-orbital re-entry velocity in those facilities for a limited time (impulse facilities)or in a limited region (stagnation point in plasma wind tunnels). Each category offacility is addressing a specific aspect of the flowfield as it is sketched in Fig. 9.11.The velocity of the flow is much higher in the case of high-speed re-entry whichconcern super-orbital conditions. This results in considerably higher free-streamFig. 9.11 Different types of facilities for flow radiation or gas–surface interaction studies9 Aerospace Flight Modeling and Experimental Testing 145Fig. 9.12 Testing methodology for high-speed re-entry TPS sizingenthalpy and pressure in the flow, and leads to coupling phenomena that cannot bereproduced on scaled down models.High-speed re-entry requires therefore a new approach of ground testing. Groundtesting facilities should be used to investigate separately different phenomenon play-ing a role in the aerothermodynamic of the flow, rather than to duplicate flight, as itis the case for lower velocity re-entry.Those investigations should be performed on a panel of different facilities in orderto develop models and databases. Models of shock layer radiation can be developedbased on measurements performed in impulse facilities, under conditions similarto those encountered in high-speed re-entry. However, material processes such asablation and radiative heating cannot be performed in the same facility due to theshort test duration involved. In particular, gas–surface interaction have to be studiedin plasma wind tunnels. As it is known, those facilities could produce the requiredheat flux level, but have not been designed to take into account the correct couplingphenomena including the radiation processes present below a shock layer. Modelsof gas–surface interaction have therefore to be developed under different conditionsthan that of high-speed re-entry.Those models, developed separately, should then be implemented in Computa-tional FluidDynamics (CFD) codes and allow extrapolatation to the actual flight con-ditions. Since they cannot be validated within the flight conditions envelope unlessflight-testing is performed, they need to capture the main physical phenomena andtheir accuracy is of prime importance. This, in turn, allows sizing and designingre-entry probes as well as assessing their performance in flight, with computationalmethods rather than direct ground testing facilities. The general framework of thetesting methodology for high-speed re-entry is summarized in Fig. 9.12.146 O. Chazot9.7 ConclusionThe reproduction of high-speed re-entry conditions in ground-based facilities is areal challenge. Because of the coupling phenomena that characterized this type offlows their accurate experimental simulation on scaled-downmodels is impossible inground-based facilities. It appears as well that only little research has been conductedon dedicated strategies for ground testing of high-speed re-entry flows. Most of thetime, testing is limited to qualification tests that reproduce the heat flux level withouttaking into account all the physical phenomena involved.The two main limitations concern radiation and gas–surface interactions: bothcannot be correctly reproduced at the same time in a single facility. Indeed, the time-scale achieved in impulse facilities is shorter than those relevant for gas–surfaceinteractions, while the correct radiation phenomena are difficult to reproduce inplasma wind tunnels.A solution is to develop models for radiation and gas–surface interaction in therelevant facilities, under controlled environments. Since the conditions in whichthose models can be validated in ground facilities are different from the one encoun-tered in high-speed flows, they specifically require to be physic based in view of theirextrapolation to flight conditions. In such a context, model validation and their incor-poration in CFD codes are crucial for the development of aerospace applications. UQmethodologies are expected to play a major role for this approach as they representthe unique way to give solid basis to the validation process. In order to manifest alltheir benefit and correctly address the problem, UQmethods imperatively need to bearticulated with the experimental procedure. To this end, this brief review exposesthe rationale of experimental testing and how it is linked to the physics of aerospaceflights. It would serve as a basis for the development of Uncertainty Quantificationapply to the validation of high-speed flow modeling.References1. P. Barbante and O. Chazot. Flight extrapolation of plasma wind tunnel of stagnation regionflowfield. J. Thermophys. Heat Transfer, 20(3):493–499, 2006.2. P. L. Chambré and A. Acrivos. On chemical surface reaction in laminar boundary layer flows.J. Appl. Phys., 27(11):1322–1328, 1956.3. O. Chazot, F. Panerai, and V. Van Der Haegen. Aerothermochemistry testing for lifting reentryvehicles. In ESA, editor, 7th European Symposium on Aerothermodynamics, 2011.4. O. Chazot, R. Régnier, and A. Garcia Munõz. Simulation methodology in plasmatron facilityand hypersonic wind tunnels. In 12th International Conference on Method of AerophysicalResearch, 2004.5. G. de Crombrugghe, R. Morgan, and O. Chazot. Theoretical approach and experimental ver-ification of diffusive transport under binary scaling conditions. International Journal of Heatand Mass Transfer, (97):675–682, 2016.6. R. J.-A. Desideri, Glowinski, and J. Periaux. Hypersonic Flows for Re-entry Problems.Springer-Verlag, Berlin, 1991. Viviand, H., Similitude in Hypersonic Aerodynamics, pp.72-97.7. D. Elligton. Binary scaling limits for hypersonic flight. AIAA J., 5(9):1705–1706, 1967.9 Aerospace Flight Modeling and Experimental Testing 1478. J.A. Fay and F.R Riddell. Theory of stagnation point heat transfer in dissociated air. J. Aeron.Sci., 25(2):73–85, 1958.9. N.C. Freeman. Non-equilibrium flow of an ideal dissociating gas. J. Fluid Mech., 4(part 4),1958.10. R. Goulard. The coupling of radiation and convection in detached shock layers. Journal ofQuantitative Spectroscopy and Radiative Transfer, 1.11. W. D. Hayes and R. F Probstein. Hypersonic Flow Theory. Academic Press, 1959.12. W.D.Hayes andR. F. Probstein. Viscous hypersonic similitude. J. Aeron. Sci., 26(12):815–825,1959.13. H. Hornung. Experimental real gas hypersonics. Aeronaut. J., 92:379–389, 1988.14. Kolesnikov A.F. Condition of simulation of stagnation point heat transfer from high enthalpyflow. Fluid Dynamics, 28(1):131–137, 1993.15. W. Kordulla, X. Bouis, and G. Eitelberg. Wind tunnels for space applications at dlr and onera.In Proceeding AAAF, 12th European Aerospace Conference, 1999.16. F. K. Lu and D. E. Marren. Advanced Hypersonic Test Facilities, Prog. in Astronaut. andAeronaut., volume 198. AIAA, Reston, MA, 2002. Principles of Hypersonic Test FacilityDevelopment, pp. 17-25.17. Parck C. Laboratory simulation of aerothermodynamic phenomena: a review. In AIAA, editor,17th Aerospace Ground Testing Conference, 1992.18. Richard Morgan, Tim McIntyre, and al. Impulse facilities for the simulation of hypersonicradiating flows. In AIAA, editor, 38th Fluid Dynamics Conference and Exhibit, 2008.19. Rosner D.E. Chemically frozen boundary layer with surface reaction. Princeton University,Princeton N.J., 1958. March.20. W. Saric, J. Muyalert, and C. Dujaric. Hypersonic experimental and computational capability.improvement and validation. Technical report, AGARD-AR-319, Vol. I, 1996.21. E. Venkatapathy and al. Thermal protection system technologies for enabling future samplereturn missions. Technical report, NASA Ames Research Center, 2009.Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you give appropriatecredit to the original author(s) and the source, provide a link to the Creative Commons license andindicate if changes were made.The images or other third party material in this chapter are included in the chapter’s CreativeCommons license, unless indicated otherwise in a credit line to the material. If material is notincluded in the chapter’s Creative Commons license and your intended use is not permitted bystatutory regulation or exceeds the permitted use, you will need to obtain permission directly fromthe copyright holder.

Uncertainty

Abstract

Extracted data

Uncertainty

Abstract

Extracted data

Related items

Related items