Solving Stochastic Inverse Problems using Sigma-Algebras on Contour Maps
Troy Butler, Don Estep, Simon Tavener, Timothy Wildey, Clint Dawson, Lindley Graham
aa r X i v : . [ m a t h . NA ] J u l SIAM/ASA J. U
NCERTAINTY Q UANTIFICATION c (cid:13) xxxx Society for Industrial and Applied MathematicsVol. xx, pp. x x–x Solving Stochastic Inverse Problems using Sigma-Algebras on Contour Maps
T. Butler ∗ , D. Estep † , S. Tavener ‡ , T. Wildey § , C. Dawson ¶ , and L. Graham k Abstract.
We compute approximate solutions to inverse problems for determining parameters in differentialequation models with stochastic data on output quantities. The formulation of the problem andmodeling framework define a solution as a probability measure on the parameter domain for agiven σ − algebra. In the case where the number of output quantities is less than the number ofparameters, the inverse of the map from parameters to data defines a type of generalized contourmap. The approximate contour maps define a geometric structure on events in the σ − algebra forthe parameter domain. We develop and analyze an inherently non-intrusive method of sampling theparameter domain and events in the given σ − algebra to approximate the probability measure. Weuse results from stochastic geometry for point processes to prove convergence of a random samplebased approximation method. We define a numerical σ − algebra on which we compute probabilitiesand derive computable estimates for the error in the probability measure. We present numericalresults to illustrate the various sources of error for a model of fluid flow past a cylinder. Key words.
1. Introduction.
There are a seemingly unlimited number of proposed methods for solv-ing forward and inverse problems of interest in the scientific community. The computationalmodel considered and the framework assumed in formulating the forward and inverse prob-lems play a critical role in defining the solution methodology. In this work, we consider aset-approximation method for approximating solutions to a stochastic inverse problem fordeterministic models formulated within a measure-theoretic framework [7]. The solutions weseek are non-parametrically defined probability measures on the model parameter domaindefined by uncertain model inputs, e.g. viscosity of a fluid, initial or boundary conditions,and/or domain parameters such as locations of holes.The models we consider are deterministic physics-based models for which a limited numberof physical observations of the solution are available. We call the observations quantities ofinterest (QoI) and are particularly interested in the situation where there are a smaller numberof QoI than model parameters. In this case, solution of the deterministic inverse problem iscomplicated by the fact that the map defined by inverting the QoI is a set-valued map. Thisis often referred to as ill-posedness. Common methods such as regularization effectively alterthe QoI map so that it has a well-posed deterministic inverse albeit with an altered geometric ∗ Department of Mathematical and Statistical Sciences, University of Colorado Denver, Denver, CO 80202(
[email protected] ). † Department of Statistics, Colorado State University, Fort Collins, CO 80523 ( [email protected] ). ‡ Department of Mathematics, Colorado State University, Fort Collins, CO 80523( [email protected] ) § Sandia National Labs, Albuquerque, NM 87185 ( [email protected] ). Sandia is a multiprogram laboratoryoperated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin company, for the United StatesDepartment of Energy’s National Nuclear Security Administration under Contract DE-AC04-94-AL85000. ¶ Institute for Computational Engineering and Sciences (ICES), University of Texas at Austin, Austin, Texas 78712( [email protected] ). k Institute for Computational Engineering and Sciences (ICES), University of Texas at Austin, Austin, Texas 78712( [email protected] ). structure. Introducing stochasticity in the form of probability measures on inputs and/oroutputs also fails to address this fundamental geometric issue . The method we use is basedon a computational measure-theoretic approach that fully exploits the set-valued inversesdirectly preserving specific geometric information contained in the map from inputs to outputs[2, 6, 5, 7]. Approximate solutions are defined using simple function approximations to thedensity of the unique probability measure on the parameter domain and given σ − algebra.There are two common steps in computing any simple function approximation to a mea-surable function such as the density (i.e. Radon-Nikodym derivative) of a probability measure:(1) identify measurable sets in a given σ − algebra partitioning the domain of the function;(2) assign a nominal function value on each of these sets.These steps are explicitly used in classical examples applying the Lebesgue Monotone Con-vergence Theorem. First, a specific sequence of partitions on the range of a given measurablereal-valued non-negative function is identified. The preimages of the sets in each partitionare used to identify the measurable sets in the σ -algebra on the domain of the function. Thiscompletes the first step and requires evaluation of the inverse map to identify inverse sets.The second step is trivial in this case where the infimum of each output set is assigned to itspreimage to define the simple function approximations for each partition.A probability measure with density on the parameter domain and original σ − algebrais a solution to the stochastic inverse problem if the QoI maps this parameter density tothe output density. Unique solutions to the stochastic inverse problem exist in a contour σ − algebra [7] embedded in the original σ − algebra. Moreover, combining the DisintegrationTheorem with an Ansatz describing probabilities on the contour events proves unique solutionsto the stochastic inverse problem exist in the original parameter σ − algebra [7]. As shown in[2, 7], the key step in constructing any simple function approximation to the density solvingthe stochastic inverse problem is the approximation of the QoI contour map in the given σ − algebra. The contour map contains all the geometric information available by the QoImap. It is then straightforward to assign probability values to any approximate partitioningsets of the contour events in the original σ − algebra [2, 6, 5, 7] such that a sequence convergesusing either the Monotone Convergence Theorem [2] or the Lebesgue Dominated ConvergenceTheorem [7].The formulations of the forward and inverse problems within this measure-theoretic frame-work, the existence and uniqueness of solutions, and the approximations of solutions by de-terministic approximation techniques have been recently studied [2, 6, 5, 7]. The focus of thiswork is the development and analysis of sampling techniques within a parameter domain toimplicitly define approximating sets of both a contour σ − algebra and the original σ − algebra.In this way we create a new computational algorithm for the probability measure based onthe samples and their approximation properties. We prove that a counting measure approxi-mation converges to the exact probability measure solving the stochastic inverse problem asthe number of samples increases. Finally, for the computed probabilities of events, we derivecomputable error estimates and bounds for the effects of using a finite number of samples and Statistical inverse problems where the QoI map is replaced by a statistical map constitute an entirelydifferent framework for formulating inverse problems, see [7, 2] for more detailed relations to statistical andBayesian inference problems. igma-Algebras on Contour Maps 3 numerical errors in the evaluation of the model.The outline of this paper is as follows. In Section 2, we summarize the measure-theoreticframework and formulation of the stochastic inverse problem. In Section 3, we summarizethe theoretical results involving σ − algebras on contour maps and solutions to the stochasticinverse problem. We then present computational algorithms including a point-sample basedalgorithm. We describe the implicit set-approximation properties defined by the sampling ofpoints in parameter space and define a counting measure approximation to the inverse proba-bility measure. In Section 4, we prove that the counting measure approximation converges tothe correct probability measure for arbitrary events in the original σ − algebra. In Section 5, wesummarize the types and sources of error in the computed probabilities and how we estimateand/or bound these errors. In Section 6, we provide numerical results for a Navier-Stokes flowpast a cylinder with uncertain viscosity and cylinder location.
2. Formulating stochastic inverse problems in a measure-theoretic framework.
LetΛ ⊂ R n denote the parameter domain for the model. We assume Λ is a metric space whosemetric is specified as part of the model. Thus, there is an induced topology on Λ and Borel σ − algebra B Λ . We define the measure space (Λ , B Λ , µ Λ ) using the measure µ Λ induced by themetric [15]. Let Q : Λ → D denote the vector-valued QoI map where D ⊂ R m , m ≤ n , isdefined by functionals of solution to the model. Definition 2.1.
We say that the component maps of m -dimensional locally differentiablevector-valued map Q ( λ ) are geometrically distinct (GD) if the Jacobian of Q has full rankat every point in Λ . In practice, since we deal with measurable events of non-zero measure, we can weakenlocal differentiabillity of Q and/or the full rank of the Jacobian of Q to hold at a.e. point inΛ . Assuming Q is locally differentiable with GD components, µ Λ defines a “push-forward”measure µ D on ( D , B D ), where B D is the Borel σ − algebra on D , and for any A ∈ B D , µ D ( A ) = µ Λ ( Q − ( A )). This yields the measure space ( D , B D , µ D ). The measures µ Λ and µ D are volumemeasures not probability measures. We often assume the probability measures are absolutelycontinuous with respect to the volume measures in which cases we typically refer to theassociated probability density functions (i.e. the Radon-Nikodym derivatives).The inverse sensitivity problem of interest is the direct inversion of the forward stochasticsensitivity analysis problem . In other words, the stochastic inverse problem is to determineprobability measures on measurable space (Λ , B Λ ) such that the push-forward probability mea-sures match given probability measures on ( D , B D ). Fundamentally, this requires a descriptionof the mapping between sets in the various σ − algebras B D and B Λ .When m < n , the Implicit Function Theorem guarantees that for any fixed point in D ,there exists some union of locally continuous ( n − m )-dimensional manifolds defining the set-valued inverse. We use the notation and definitions of [7, 2, 6, 5] and define a generalizedcontour as any such set-valued inverse of the map Q . We assume that D is defined by Q (Λ), Since we use probability densities, which are L functions, to compute probabilities of events, any set ofzero measure has no affect on the final solutions and can simply be ignored or removed from Λ in the formulationof the inverse problem for convenience. The forward problem is standard: given probability measure P Λ on measurable space (Λ , B Λ ) compute theprobability measure P D on measurable space ( D , B D ). T. Butler et. al. i.e. D is the range of the QoI map containing all possible physically observable data that canbe mapped to from Λ. Thus, we can decompose Λ into a union of generalized contours in 1-1correspondence with the points in D . In other words, the map Q defines a type of generalizedcontour map on Λ. There exists a (possibly piecewise-defined) continuous m -dimensionalindexing manifold, called a transverse parameterization , in Λ defining a bijection between D and the generalized contours, see [7] for more details. As way of analogy, a specific transverseparameterization is like a particular path of ascent up a mountain that a hiker plots out usinga contour map of elevations.We note that the contour map obtained from Q − defines an equivalence class relation onΛ. Denote this space of equivalence classes as L where each point in L corresponds to a setof points in Λ. We obtain a measure space ( L , B L , µ L ), where the σ − algebra B L on L can begenerated using inverse images of a collection of Borel sets in B D and the volume measure µ L induced by µ D .We exploit the equivalence relation to define the induced σ − algebra C Λ on Λ that canbe generated from the set of equivalence classes for a set of generating events in B L . For m < n , C Λ is a proper subset of B Λ . We define C Λ as the contour σ − algebra on Λ and callevents in C Λ contour events . The geometric structure of Q on Λ is fully exploited to define thecontour events in the contour σ -algebra C Λ . We emphasize that events in C Λ can be uniquelydetermined by events in B D .The goal is to define solutions to the stochastic inverse problem on (Λ , B Λ ), e.g. in termsof a density defined on points in Λ not points on a contour map . Below, we prove there existsa σ − algebra defined on the generalized contour map on Λ equivalent to the σ − algebra B Λ on Λ. We then use this σ − algebra equivalency in the Disintegration Theorems that follow toprove the existence and uniqueness of solutions to the stochastic inverse problem.
3. Solving the stochastic inverse problem using sigma-algebras on contour maps.
Acommon starting point for constructing a measure on a σ − algebra is to first define an algebraused to generate the σ − algebra of interest. The next steps can vary, but typically we define apremeasure on the algebra which induces an outer-measure in a natural way. Finally, employCarath´eodory’s Theorem to extend the outer-measure uniquely to a complete measure on thegenerated σ − algebra. We can define an algebra with elementary families of elementary setsfor which a clear notion of measure is defined, e.g. by using a partitioning of a space where setsare either disjoint or possibly intersect only at boundaries such as hypercubes or generalizedrectangles in R n . We will exploit such generating sets in much of the theory below.In Section 3.1, we relate σ − algebras on Λ to σ − algebras on the generalized contour maps.In Section 3.2, we summarize the theory of the existence and uniqueness of solutions to thestochastic inverse problem with respect to these σ − algebras. In Section 3.3, we describe theapproximation of solutions. The generalized contours define atype of contour map on Λ. Given the geometry of these generalized contours and a topology,we can define σ − algebras on contour maps (SACOM). In Section 3.2 below, we show that thesolution to the stochastic inverse problem can be solved uniquely using SACOM. The goal isto relate the σ − algebras on Λ (e.g. C Λ or B Λ ) to SACOM.To motivate what follows, consider the case where the map Q is linear. In this case, igma-Algebras on Contour Maps 5 the ( n − m )-dimensional generalized contours are hyperplanes. For simplicity, let Λ = R n ,and assume we use the typical Euclidean distance metric and induced topology giving theBorel σ − algebra on the product space Λ. Let π ℓ ( λ ) = ℓ denote the projection map from Λto the equivalence class ℓ ∈ L containing λ . For arbitrary ℓ ∈ L , let C ℓ := π − ( ℓ ) denotethe associated generalized contour as a subset of Λ. Changing coordinates with respect tothe directions orthogonal and parallel to the hyperplanes gives B Λ = B L ⊗ B C ℓ . Identifyingsuch a product decomposition is useful for applying Fubini’s theorem in order to integratecertain functions using iterated integrals. The key is that the domain of integration, which isa set in the original σ − algebra, can be “cut” into products of lower-dimensional sets in thecomponent measure space σ − algebras. In this case, such sets are “cut” into products of setsboth along and transverse to the generalized contours. Below, we extend this for the generalcase of nonlinear generalized contours indexed by a nonlinear manifold defining a transverseparameterization.Let F C ℓ denote a σ − algebra on a given C ℓ . Below, we describe two natural choices for F C ℓ for each C ℓ defining a family of measurable spaces { C ℓ , F C ℓ } indexed by ℓ ∈ L .Any space contains the so-called trivial σ − algebra consisting of the empty set or the entirespace. Therefore, there exists at least one σ − algebra on each C ℓ , given by T C ℓ = { C ℓ , ∅} .Alternatively, using the induced topology on each C ℓ , we define Borel σ − algebras B C ℓ . Theseare equivalently defined by { A ∩ C ℓ | A ∈ B Λ } . In other words, the Borel σ − algebra on C ℓ canbe defined by the restriction of Borel measurable sets in B Λ to the given generalized contour.These two particular choices for F C ℓ prove useful in the Disintegration Theorems of Section 3.2for showing existence and uniqueness of solutions to the stochastic inverse problem.As in many of the classical proofs of theorems involving product σ − algebras, we makeexplicit use of the generating sets for each σ − algebra. Since each measurable space we consideris metrizable, we assume the generating sets are taken from the implied topology. For each ℓ ∈ L , let F C ℓ denote any family of subsets of C ℓ generating F C ℓ . A standard measure-theoryresult states that F C ℓ is the unique smallest σ − algebra generated by F C ℓ for each ℓ ∈ L .Similarly, we let F L denote a generating set for the Borel σ − algebra B L . Definition 3.1.
We define the transverse product σ − algebra as the smallest σ − algebraon Λ generated by S ℓ ∈ A B C ℓ where A ∈ F L and B C ℓ ∈ F C ℓ for all ℓ . We denote this σ − algebraby N B L {F C ℓ } . Theorem 3.1. C Λ = N B L {T C ℓ } .Proof: This follows immediately from the definition of contour σ − algebra C Λ using theequivalence relation determined by Q − and the generating sets for B L and {T C ℓ } . (cid:3) Theorem 3.2. B Λ = N B L {B C ℓ } . Before we prove this theorem, we note that the inclusion N B L {B C ℓ } ⊂ B Λ is not at allobvious since A ∈ F L may be uncountable and σ − algebras are closed under countable unions.Since we assume Q is locally differentiable, there exists a bijection between L and a transverseparameterization in Λ [2, 7]. Since the transverse parameterization is a piecewise-continuous m -dimensional manifold in separable space Λ, it follows that L is separable. This is used below. See Propositions 1.5 in [14] identifying Λ = R n and from the change of coordinates L = R m and C ℓ = R ( n − m ) . See Lemma 6.2.4 of [1].
T. Butler et. al.
Also, we let E A := π ℓ ( A ) for any A ∈ B Λ . We use this notation elsewhere as convenient. Proof:
Let F Λ denote any Borel generating set of B Λ and consider any A ∈ F Λ . Restrictthe map Q to A , then by the assumption of GD component maps, there exists an ( n − m )-dimensional (piecewise) continuous manifold defining a transverse parameterization on A [7].Moreover, the transverse parameterization is a Borel set in Λ by assumption of the localdifferentiability of Q and is in 1-1 correspondence with E A . Thus, we have that the Borel set A defines the unique Borel set E A ⊂ L . Restricting any generating sets of B Λ to a generalizedcontour C ℓ defines a generating set for B C ℓ . From Definition 3.1, we have that A ∈ N B L {B C ℓ } .It follows that B Λ ⊂ ⊗ B L {B C ℓ } .Suppose we are given Borel generating set F L and family of Borel generating sets { F C ℓ } defined by all the open sets on these spaces. For all ℓ ∈ L , any Borel set on C ℓ is also Borel inΛ. By assumption, Q − ( A ) is Borel in Λ for any A ∈ F L . If A is countable, then S ℓ ∈ A B ℓ ∈ B Λ for any B ℓ ∈ F ℓ . Suppose A is uncountable. Let C L denote a countable dense set in L and E L the collection of open balls in L with rational radius and center in C L . Then every openset in L is a countable union of open balls in E L . In other words, A is countably generatedby E L , so the set S ℓ ∈ A B C ℓ is countably generated by unions of Borel sets and is itself Borel.Thus, S ℓ ∈ A B C ℓ ∈ B Λ , so N B L {B C ℓ } ⊂ B Λ . (cid:3) To solve the stochastic inverse problem, weuse a form of the Disintegration Theorem [7, 9, 11], which is a powerful theoretical tool forrigorous definition of conditional probabilities. We focus on convenient forms of the theoremwritten for probabilities and direct the interested reader to [7] for a more thorough presenta-tion. Using Theorem 3.1 and following the steps of [7], we have,
Theorem 3.3 (Disintegration of Contour Map Probabilities).
Let (Λ , C Λ ) be a measurable space.Assume that P Λ is a probability measure on (Λ , C Λ ) . There exists a family of conditional prob-ability measures { P ℓ } on { ( C ℓ , T Λ ) } giving the disintegration, P Λ ( A ) = Z π L ( A ) (cid:18) Z π − L ( ℓ ) ∩ A dP ℓ ( λ ) (cid:19) dP L ( ℓ ) , ∀ A ∈ C Λ . (3.1)It is clear we can compute the induced probability measure P Λ on (Λ , C Λ ) defined by P Λ ( A ) = P L ( π L ( A )) = P D ( Q ( A )), and P L is defined by a probability density ρ L on ( L , B L )with respect to µ L , P L ( A ) = Z π L ( A ) ρ L dµ L = Z Q ( A ) ρ D dµ D = P D ( Q ( A )) , ∀ A ∈ B L . From Theorem 3.3, the conditional probability measures for P Λ , C Λ are given by P ℓ ( A ) = 1 if π − ( ℓ ) ⊂ A ∈ C Λ and P ℓ ( A ) = 0 otherwise. This proves the following Theorem 3.4.
The stochastic inverse problem has a unique solution on (Λ , C Λ ) . The primary goal is not a probability measure on (Λ , C Λ ), but a probability measure P Λ on (Λ , B Λ ). Since we will often work with densities given as Radon-Nikodym derivatives(i.e. dP Λ /dµ Λ ), we find useful the following corollary of the Disintegration Theorem [7], See Theorem 6.9.7 [1]. See Lemma 1.1 of [14]. igma-Algebras on Contour Maps 7
Corollary 3.1 (Disintegration of the Volume Measure).
There exists a family of volume mea-sures { µ C ℓ } on { ( C ℓ , B C ℓ ) } such that for any A ∈ B Λ , µ Λ ( A ) = Z E A µ C ℓ ( π − L ( ℓ ) ∩ A ) dµ L ( ℓ ) = Z E A Z π − L ( ℓ ) ∩ A dµ C ℓ ( λ ) dµ L ( ℓ ) . Using Theorem 3.2, we arrive at the more common form of the Disintegration Theoremfor probability measures given by
Theorem 3.5.
Let (Λ , B Λ ) be a measurable space. Assume that P Λ is a probability measureon (Λ , B Λ ) . There exists a family of conditional probability measures { P ℓ } on { ( C ℓ , B C ℓ ) } giving the disintegration, P Λ ( A ) = Z π L ( A ) (cid:18) Z π − L ( ℓ ) ∩ A dP ℓ ( λ ) (cid:19) dP L ( ℓ ) , ∀ A ∈ B Λ . (3.2)In the above Disintegration Theorems, the conditional probability measures can be ex-tended to (Λ , C Λ ) or (Λ , B Λ ) by extending the measures P ℓ ( A ) = 0 for all A ⊂ Λ \ C ℓ foreach ℓ ∈ L . Theorem 3.5 guarantees that any probability measure on (Λ , B Λ ) can be decom-posed into a form involving a probability measure on ( L , B L ) uniquely defined by P D andprobability measures on each measurable generalized contour space ( C ℓ , B C ℓ ) defined by theconditional probabilities P ℓ . It is not obvious what these conditional probability measuresare in Theorem 3.5. Clearly, any conditional probability measures on { ( C ℓ , B C ℓ ) } can not bedetermined by observations of Q ( λ ) ∈ D . This motivates the adoption of an Ansatz in whichthe probability measures along { ( C ℓ , B C ℓ ) } are specified. Ansatz:
For all ℓ ∈ L , assume a probability measure P ℓ ( · ) is given on ( C ℓ , B C ℓ ). Theorem 3.6.
Under the Ansatz, the stochastic inverse problem has a unique solution on (Λ , B Λ ) . We prefer a specific Ansatz that can be interpreted as prescribing a “non-probabilistic” or“non-preferential” weighting determined by the disintegration of volume measure to computeprobabilities of events inside of a contour event . However, the approximation method andresulting Algorithm 1 (summarized below) can be easily modified for any Ansatz (see [7] formore details). Standard Choice for Ansatz: P ℓ = µ C ℓ /µ C ℓ ( C ℓ ) , ∀ ℓ ∈ L . If we do not specify an Ansatz, the stochastic inverse problem can only be solved on (Λ , C Λ )where there is only one possible probability measure on each ( C ℓ , T C ℓ ). From Theorems 3.1and 3.2, it is evident that the choice of σ − algebra on Λ dictates the SACOM and whether ornot we must adopt an Ansatz for the solution. This is valid when µ Λ (Λ) < ∞ , which is generally satisfied in practice using compact domains. If this isnot the case, implying certain physical parameters are unbounded, then we may employ techniques similar tothose used in Bayesian analysis for “non-informative” priors. T. Butler et. al.
In [7], we describe several approximation issues that needto be addressed in any practical computation of P Λ . As described below, the fundamental ap-proximation issues of any measure are the approximation of events in the various σ − algebras.The numerical approximation issues involve error in numerical evaluation of the model andcomputation of probability measures on some collection of events. Here, we discuss the eventapproximations. In Section 5, we analyze the effect of errors on the approximate probabilitymeasures computed from the approximating sets. Brief review of approximations.
Repeated application of the Lebesgue Dominated Conver-gence Theorem yields the following,
Theorem 3.7.
Given probability measure P D , absolutely continuous with respect to µ D , on ( D , B D ) with density ρ D and event A ∈ B Λ , there exists a sequence of approximations P Λ ,N ( A ) using simple function approximations to probability densities ρ Λ ,N and ρ D ,M requiring onlycalculations of volumes in Λ that converges to P Λ ( A ) as N, M → ∞ . The proof of Theorem 3.7 details approximating probability densities ρ Λ and ρ D (i.e. theRadon-Nikodym derivatives of the probability measures) in order to apply the LebesgueDominated Convergence Theorem and outlines the computational measure-theoretic Algo-rithm 1 (see [7] for more details). Unsurprisingly, the first step is the generation of partitions { I i } Mi =1 ⊂ D and { b j } Nj =1 ⊂ Λ such that (1) arbitrary events in B D and B Λ are approximatedby unions of sets from these partitions, and (2) simple function approximations to the densi-ties can be computed on these partitions. These partitions can be chosen from an algebra ofelementary sets.Since densities are L functions, we can use convex sets with continuous boundaries todefine a finite partition { I i } Mi =1 ⊂ D and associated simple function approximation ρ D ,M ofsufficient accuracy in the L -norm . The local differentiability of Q implies that Q − ( I i ) is ameasurable event in both C Λ and B Λ with boundary of zero µ Λ -measure. Such events can havetheir measures approximated within any prescribed tolerance by a finite set of hypercubes .This defines an obvious choice for the cells { b j } Nj =1 partitioning Λ. Once the probabilities ofcells { b j } ⊂ Λ have been approximated by Algorithm 1, we may estimate P Λ ( A ) for arbitraryevent A ∈ B Λ in any of the usual measure-theoretic ways such as using inner or outer sums,averages of inner and outer sums, or direct integration of the simple function approximation ρ Λ ,N .Suppose instead that we wish to solve the stochastic inverse problem on (Λ , C Λ ) instead of(Λ , B Λ ). In this case, we can approximate the solution using the same approximate solutionon (Λ , B Λ ) computed using Algorithm 1 since C Λ ⊂ B Λ . In other words, we may use the sameapproximate finite generating set { b j } for B Λ as an approximate finite generating set for C Λ .Unless otherwise stated, we let { I i } Mi =1 denote a partition of D used to approximate eventsin B D and define the simple function approximation ρ D ,M = P Mi =1 p i I i ( q ), where p i = P D ( I i ).In other words, the probability on each set I i ⊂ D is defined by the expected value of I i ( q )on D . This is a natural choice equivalent to using the Integral Mean Value Theorem to definea normalized approximation, i.e. ρ D ,M is a density with this choice of p i . See Theorem 2.41 in [14]. See Theorem 2.40 of [14]. igma-Algebras on Contour Maps 9
Algorithm 1
Approximation of the Inverse DensityGenerate partitions { I i } Mi =1 ⊂ D and { b j } Nj =1 ⊂ ΛFix and normalize the simple function approximation ρ D ,M = P Mi =1 p i I i ( q )Let { A i } Mi =1 ⊂ Λ denote the induced regions of generalized contours partitioning Λ for j = 1 , . . . , N dofor i = 1 , . . . , M do Compute µ Λ ( A i ∩ b j ) and store as ij -component in matrix V . end forend forfor j = 1 , . . . , N do Set P Λ ,N ( b j ) to P Mi =1 p i ( V ij / P Nj =1 V ij ) end for A sample based approximation and counting measure.
Previous implementations of Algo-rithm 1 used regular grids so that { b j } Nj =1 was defined as a set of generalized rectangles orhypercubes [2, 6, 5]. While such sets may be refined and the measure of unions of these setsmade to approximate any Borel set arbitrarily well [14], there are obvious practical difficultieswith using such an approximating set in high dimensions. Here, we consider an alternativewhere the sets { b j } Nj =1 are defined implicitly by a finite collection of samples in Λ satisfyingsome particular properties detailed below. Definition 3.2.
For a fixed number of samples (cid:8) λ ( j ) (cid:9) Nj =1 ⊂ Λ , there is a Voronoi tessella-tion of Λ denoted by {V j } Nj =1 ⊂ Λ defined by V ( λ ( j ) ) := n λ ∈ Λ : d v ( λ ( j ) , λ ) ≤ d v ( λ ( i ) , λ ) , ∀ i = 1 , . . . , N o . Here, d v ( · , · ) denotes a metric on Λ used to define the Voronoi cells . Clearly the goal is to define a set of samples (cid:8) λ ( j ) (cid:9) Nj =1 implicitly defining a tessellation ofΛ that is useful for approximating events, i.e. measurable sets, in B Λ . We often approximateevents in a given σ − algebra by a Voronoi coverage. Definition 3.3.
We say that A N is the Voronoi coverage of A ∈ B Λ and µ k ( A ) its volumedefined by A N := [ λ ( j ) ∈ A, ≤ j ≤ N V ( λ ( j ) ) , and µ k ( A ) = µ Λ ( A k ) = N X j =1 µ Λ ( V ( λ ( j ) )) λ ( j ) ∈ A . Definition 3.4.
A rule for defining any N samples (cid:8) λ ( j ) (cid:9) Nj =1 ⊂ Λ is called B Λ -consistent if µ Λ ( A △ A k ) → as N → ∞ , ∀ A ∈ B Λ , s.t. µ Λ ( ∂A ) = 0 . The metric d v ( · , · ) is possibly different from the metric that induces the volume measure µ Λ and Borel σ -algebra B Λ . Common choices for d v ( · , · ) are the standard Euclidean 2-norm or 1-norm. Clearly any B Λ -consistent rule implies that any generalized rectangle or finite union ofgeneralized rectangles may be approximated arbitrarily well by a Voronoi coverage defined bya finite number of samples. Replacing { b j } Nj =1 with a Voronoi tesselation {V j } Nj =1 generatedfrom a B Λ -consistent rule in Algorithm 1, we define the following, Definition 3.5.
Let ˜ P Λ ,N denote the counting probability measure (or simply countingmeasure on (Λ , B Λ ) defined by ˜ P Λ ,N ( A ) := N X j =1 P Λ ,N ( V j ) λ ( j ) ( A ) . Note that ˜ P Λ ,N is simply a measure of point masses at each sample from (cid:8) λ ( j ) (cid:9) Nj =1 . Atechnical point is that we may use a counting measure to estimate the probability of any A ∈ B Λ . However, we only prove convergence to the exact probability when µ Λ ( ∂A ) = 0.It is possible to have compact Borel sets with boundaries of non-zero measure. We describein Section 4.1 how to approximate the probability of such sets. In practical computations,we are not concerned with such sets with fractal boundaries. Furthermore, such sets donot enter into any of the computational algorithms. This is due to the assumption of localdifferentiability of the map Q , which implies that all generalized contours, and subsequently,all induced regions of generalized contours defined by Q − ( I ) for any Borel set I ∈ B D aresets with boundaries of zero µ Λ -measure. A Monte Carlo approach.
There are many possible B Λ -consistent rules that we may choose.For example, we may define a sequence of uniformly refined grids of which the grid points arenumbered and sequentially sampled in a serpentine manner. This choice produces Voronoi cellsthat are generalized rectangles (for certain choices of N ) and has been used previously, e.g. see[2, 6, 5]. While this produces small Voronoi cells everywhere in the domain, the approximationproperties of the cells have negative dependence on the dimension of the parameter space interms of requiring large N to achieve reasonable µ Λ -approximations by the Voronoi coverage.We may choose to normalize the volume measure (if possible) and generate N i.i.d. samplesfrom this distribution corresponding to a “uniform” sampling density (see Section 4.2 below).This amounts to a standard Monte Carlo (MC) approach for sampling on (Λ , B Λ , µ Λ ). Fur-thermore, we may use the standard MC approximation that all Voronoi cells have the samevolume which greatly reduces the computational demands of the algorithm, see Algorithm 2.The counting measure definition is consistent with the MC approach of Algorithm 2 toestimate the probabilities of individual Voronoi cells V j inside of a particular induced regionof generalized contours by counting the number of samples within this event. Specifically, inthe algorithm, c ( i o ( j )) counts the number of Voronoi cells we associate within the inducedgeneralized region of contours defined by Q − ( I i o ( j ) ).Standard MC sampling algorithms tend to produce random samples that are clusteredresulting in large variations in the sizes of the corresponding Voronoi cells. This may notprovide good set-approximation properties for small sample sizes. Sampling points equidistantalong a space filling curve or using Latin hypercube sampling are other options. Or, we maysimply choose the samples deterministically based on other information we know or we wantto impose a certain resolution on the problem, e.g. we may wish to obtain a fine resolution of igma-Algebras on Contour Maps 11 Algorithm 2
A Monte Carlo Approximation of the Inverse DensityLet (cid:8) λ ( j ) (cid:9) Nj =1 ⊂ Λ denote uniform i.i.d. random samples from B Λ -consistent rule, and {V j } Nj =1 ⊂ Λ denote the associated Voronoi tessellation of Λ. for j = 1 , . . . , N do Assign a nominal value of Q j to V j , e.g. in the continuous case use Q j = Q ( λ ( j ) ). end for Generate partition { I i } Mi =1 ⊂ D .Fix and normalize the simple function approximation ρ D ,M = P Mi =1 p i I i ( q ).Initialize M × c and N × i o to zeros. for j = 1 , . . . , N do Set i = 1 and flag = 0 while i ≤ M and flag = 0 doif Q j ∈ I i then c ( i ) = c ( i ) + 1, i o ( j ) = i , flag = 1. else i = i + 1. end ifend whileend forfor j = 1 , . . . , N do Set P Λ ,N ( V j ) to p i o ( j ) / c ( i o ( j )) end for some particular subset of Λ for some design purposes. Some pseudo-MCMC sampling might bebest. The underlying theory only requires that the samples be generated from a B Λ -consistentrule. Exploring all possibilities is beyond the scope of this paper and we limit presentationto the basic theory of convergence of counting measures for B Λ -consistent rules with specialattention paid to random sampling rules satisfying criterion as discussed in Section 4.2 below.
4. Convergence of counting measures.4.1. General theory.
The proof of convergence is based upon classical measure-theoreticarguments using simple function approximations defined using some approximate generatingset to a Borel σ − algebra. We use a triangle inequality to separate the error in probabilityinto two parts involving the approximation properties of the Voronoi coverages to generalizedrectangles and the approximation properties of the generalized rectangles to other Borel-measurable sets. Lemma 4.1.
For any A ∈ B Λ and ǫ > , there exists simple function approximation ρ D ,M and K < ∞ generalized rectangles in R n partitioning compact Λ with | P Λ ,K ( A ) − P Λ ( A ) | < ǫ .Proof: Let ǫ > D into M generalized rectangles, { I i } Mi =1 , such that the L error of ρ D ,M and ρ D is less than ǫ/
3. For each A i := Q − ( I i ) ∈ B Λ ,which have boundaries of zero µ Λ -measure, there is a finite number K i of generalized rectanglespartitioning Λ such that | P Λ ,K i ( A i ) − P Λ ( A ) | < ǫ/ (3 M ). Moreover, for any A ∈ B Λ there isa finite number K A of generalized rectangles, { R j } K A j =1 in Λ such that µ Λ ( A △ S j R j ) < ǫ/ The conclusion follows from Theorem 3.7 where P Λ ,K is constructed as in Algorithm 1 with { b j } Kj =1 given by K sufficiently fine generalized rectangles partitioning Λ from which the finitenumber of generalized rectangles used above may be constructed by finite unions. (cid:3) Theorem 4.1.
Given probability measure P D , absolutely continuous with respect to µ D , on ( D , B Λ ) with density ρ D and a B Λ -consistent rule for generating samples, then for any event A ∈ B Λ with µ Λ ( ∂A ) = 0 , there exists a sequence of approximations of simple functions ρ D ,M and counting measures ˜ P Λ ,N ( A ) such that ˜ P Λ ,N ( A ) → P Λ ( A ) a.s. as N, M → ∞ .Proof:
Let ǫ > λ (1) , λ (2) , . . . denote a sequence of samples in Λ computedfrom the B Λ -consistent rule. We have by the assumptions of Q that for any fixed I i ∈ B D withcontinuous boundary that the induced region of generalized contours A i := Q − ( I i ) ∈ B Λ hasboundary with zero µ Λ -measure. By assumption of ρ D , there exists a sequence of partitions { I } Mi =1 ⊂ D such that for any fixed but arbitrary M and I i ∈ B D we have that µ Λ ( A i,N △ A i ) → N → ∞ where A i,N denotes the Voronoi coverage of A i . By construction of P Λ ,N ( V i ), thedefinition of ˜ P Λ ,N , and the fact that the Borel sets on Λ can be generated by the collection ofall generalized rectangles, there is a sufficiently fine partition of Λ by K generalized rectanglesand an N such that (cid:12)(cid:12)(cid:12) ˜ P Λ ,N ( A i ) − P Λ ,K ( A i ) (cid:12)(cid:12)(cid:12) < ǫ/ N > N .By Lemma 4.1, there is a sufficiently fine partition of Λ by K generalized rectanglessuch that | P Λ ,K ( A i ) − P Λ ( A i ) | < ǫ/
2. Let K denote the number of generalized rectanglesneeded to generate all of the above generalized rectangles by finite unions and intersections.Then by a standard triangle inequality and Lemma 4.1, we have that for partition { I i } Mi =1 asabove defining simple function approximation ρ D ,M to ρ D that (cid:12)(cid:12)(cid:12) ˜ P Λ ,N ( A i ) − P Λ ( A i ) (cid:12)(cid:12)(cid:12) < ǫ for all N > N for each 1 ≤ i ≤ M . The result extends to any A ∈ B Λ with µ Λ ( ∂A ) = 0 by applyingthe Lebesgue Dominated Convergence Theorem to a sequence of ρ D ,M constructed from asequence of partitions { I i } Mi =1 ∈ B D of generalized rectangles in D and using Theorem 3.6. (cid:3) We can approximate the probability of general A ∈ B Λ using the counting measure witherror less than any ǫ >
0. For such an ǫ >
0, we use Lemma 4.1 with ǫ/ A . Let B denote this finite union ofgeneralized rectangles. Then B has the property that µ Λ ( ∂B ) = 0 and Theorem 4.1 appliesto the set B . Let B N denote the Voronoi coverage of B with N samples in Λ satisfying (cid:12)(cid:12)(cid:12) ˜ P N, Λ ( B ) − P Λ ( B ) (cid:12)(cid:12)(cid:12) < ǫ/ ǫ >
0. Then using ( A △ B N ) = ( A △ B ) △ ( B △ B N ) and atriangle inequality, the result follows. This proves the following Theorem 4.2.
Given probability measure P D , absolutely continuous with respect to µ D , on ( D , B Λ ) with density ρ D and a B Λ -consistent rule for generating samples. For any event A ∈ B Λ and ǫ > , there exists simple function ρ D ,M , finite N samples defining countingmeasure ˜ P Λ ,N , and a set B ∈ B Λ with µ Λ ( ∂B ) = 0 such that ˜ P Λ ,N ( B △ A ) < ǫ . We present a result that can beviewed as a strong law of large numbers for Voronoi tessellations. We direct the interestedreader to [19, 17] for a more thorough exposition on the subject within the context of stochasticgeometry and point processes.
Definition 4.1. A sampling distribution , denoted by F , is any distribution absolutely con-tinuous with respect to µ Λ such that the corresponding sampling density f ( λ ) > for almostevery λ ∈ Λ . igma-Algebras on Contour Maps 13 The following Lemma from [17] has been modified to conform to our notation. To ourknowledge, [17] was the first to summarize and explicitly prove the fundamental elements ofLemma 4.1. Below, we modify and expand the proof presented in [17] to highlight detailsof the convergence and uniqueness that are of practical use, e.g., when considering designof adaptive sampling procedures to determine more accurate Voronoi coverages of specificevents [8] or in understanding issues related to a posteriori error estimates described below inSection 5.
Lemma 4.1.
Given sampling density f ( λ ) > almost everywhere on Λ , if A ∈ B Λ suchthat µ Λ ( ∂A ) = 0 , then almost surely µ Λ ( A N △ A ) → , N → ∞ . Proof.
Suppose A ∈ B Λ , µ Λ ( A ) > µ Λ ( ∂A ) = 0. Let A δ := ( A + δB (0 , ∩ Λ denotethe Minkowski sum of A and δB (0 ,
1) restricted to domain Λ, where B (0 ,
1) is the unit ballin R n . Let A δ = (( A c ) δ ) c . Since µ Λ ( ∂A ) = 0, for all ǫ > δ ( ǫ ) > ∂A ⊂ A δ ( ǫ ) \ A δ ( ǫ ) and µ Λ ( A δ ( ǫ ) \ A δ ( ǫ ) ) < ǫ .Let r ( λ ( j ) ) := max λ ∈V ( λ ( j ) ) d v ( λ, λ ( j ) ) and r N := max ≤ j ≤ N r ( λ ( j ) ). We now prove that r N → h ∈ N , let H h denote the finite set of hypercubes in R n withedge-length 1 /h partitioning compact Λ. By Kninchine’s strong law of large numbers,max H ∈H h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P Nj =1 { λ ( j ) ∈ H } N − Z H f ( λ ) dµ Λ ( λ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) → N → ∞ and min H ∈H h Z H f ( λ ) dµ Λ ( λ ) > h. Suppose r N does not converge to zero almost surely, then there exists ǫ > r N ≥ ǫ infinitely often. By construction, r N ≥ r K for any N ≥ K , and there exists afixed j such that r ( λ ( j ) ) ≥ ǫ for all k ≥ j . Choosing a set of hypercubes intersecting Λwith edge-length sufficiently small (and positive) such that the maximum distance in the d v -metric between any two points in the hypercubes is less than ǫ/ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) P Nj =1 { λ ( j ) ∈ H } N − R H f ( λ ) dµ Λ ( λ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) → λ ( i ) with i = j such that its distance is less than ǫ/ λ ( j ) . Thus r N → ǫ >
0, there exists a Voronoi coverage A N of A such that A δ ( ǫ ) ⊂ A N ⊂ A δ ( ǫ ) for all N > N ( ǫ ) almost surely, which implies µ Λ ( A N △ A ) → (cid:3) This implies that any random sampling scheme with sampling density f ( λ ) > B Λ -consistent rule with probability 1. Theorem 4.1 then applies to counting measurescomputed using random sampling as long as the sampling density is positive. If f ( λ ) = 0 on any set A ∈ B such that µ Λ ( A ) > A will result in a sequence of zeros. If A shares any volume with a region of induced generalized contours of non-zero probability, thenaccording to Theorem 3.5 , P Λ ( A ) > f ( λ ) = 1, i.e. we sample from Λ according toits volume measure. However, we may choose a non-uniform intensity to improve resolutionor accuracy of the approximate probability measure for certain induced regions of generalizedcontours, e.g. in locations where k∇ Q k is largest. Such ideas of computational accuracy witha finite number of samples and adaptive sampling are the subjects of future work.We note that in Algorithm 2, we use a standard MC approximation that can be inter-preted as approximating the volumes of Voronoi cells to be equal. We may instead opt tomore accurately estimate the volumes for each V j , in which case the algorithm more closelyresembles Algorithm 1 where we explicitly use volumes of the partition of Λ to obtain P Λ ,N .Similarly, if there exists A ∈ B Λ such that µ Λ ( A ) > f ( λ ) = 1 for λ ∈ A , then the stan-dard MC approximation that all µ Λ ( V j ) are equal for all j no longer applies. The necessarymodification is that approximation to the ratio of volumes of µ Λ ( V j ) / P { k : i o ( k )= i o ( j ) } µ Λ ( V k )replace division by c ( i o ( j )) in the last for-loop . Similar modifications are required for anon-standard choice of the Ansatz. Below, we assume such modifications to Algorithm 2 aremade as necessary to correctly compute the counting measure ˜ P Λ ,N .
5. Sources of error in the counting measure.
There are two types of error in approxi-mating probability measure P Λ using Algorithm 2: stochastic and deterministic. Stochasticerror arises from using N finite samples (cid:8) λ ( j ) (cid:9) Nj =1 to implicitly define a Voronoi tessellation {V j } of Λ. Deterministic error arises from the numerical evaluation of the map Q for each ofthe N samples. Each of these errors affects ˜ P Λ ,N in different ways. The deterministic errormay lead to misidentification of induced regions of generalized contours (possibly leading toincorrect component values of the vectors c and i o in Algorithm 2). We let ˜ P Λ ,N,h denotethe computed counting measure using numerical computations Q h approximating the map Q ,where h denotes some numerical discretization parameter. The choice of samples defines theset-approximation properties of events in B Λ in terms of unions and complements of Voronoicells (the first step in Algorithm 2). In fact, the ability to estimate induced regions of gener-alized contours is limited by the choice of samples, so it is the stochastic source of error thatwe examine first.Below, P Λ denotes the exact probability measure on (Λ , B Λ ) for the standard Ansatz given a fixed simple-function approximation to P D . In practice, when repeated experimentsand the subsequent measurements are used to empirically determine P D , such an approxima-tion would be determined by the binning of the relevant QoI. In general, we decompose theerror in the computed probability measure ˜ P Λ ,N,h as P Λ ( A ) − ˜ P Λ ,N,h ( A ) = (cid:16) P Λ ( A ) − ˜ P Λ ,N ( A ) (cid:17)| {z } I − (cid:16) ˜ P Λ ,N,h ( A ) − ˜ P Λ ,N ( A ) (cid:17)| {z } II , ∀ A ∈ ˜ B Λ ⊂ B Λ . (5.1) With any Ansatz such that the generalized contours in A have non-zero probability. In this case, the counting vector c can be entirely removed from Algorithm 2. The error analysis can be altered for a non-standard Ansatz where volume measures and counting arereplaced by the probabilistic weighting of the contour events. igma-Algebras on Contour Maps 15
Here, term I is the error in approximating the probability of event A by the counting measure˜ P Λ ,N . Term II is the error in using the numerical map Q h to identify induced regions ofgeneralized contours in Algorithm 2. The events A belong to a σ − algebra ˜ B Λ that is a subsetof the original σ − algebra B Λ . σ − algebras. We find convenient thefollowing,
Definition 5.1.
For a given set of samples on Λ , (cid:8) λ ( j ) (cid:9) Nj =1 , let B Λ ,N denote the numerical σ − algebra generated by the implicitly defined Voronoi tessellation {V j } Nj =115 . By the definition or Borel sets and Voronoi cells, B Λ ,N is a proper subset of B Λ for anyfinite N , i.e. any event in B Λ ,N is µ Λ -measurable.Definition 5.1 is illustrative of the fundamental types of events we may compute prob-abilities for using Algorithm 2 with no set-approximation errors, i.e. with no term I error.Specifically, for a fixed number of samples, a fixed approximation to ρ D , and given exact map Q in Algorithm 2, we may compute ˜ P Λ ,N ( A ) exactly for all A ∈ B Λ ,N . From Lemma 4.1, wehave that any A ∈ B Λ can be sufficiently approximated in µ Λ -measure by a Voronoi coveragealmost surely, i.e. we are almost surely guaranteed to find an element of B Λ ,N for some N that approximates A to the desired accuracy. This was exploited in the proof of Theorem 4.1that the counting measures converge. Here, we focus on the events A ∈ B Λ of non-zero˜ P Λ ,N -measure and the errors in the computed probabilities of such events.Below, we show that for certain events in B Λ that there are, in a probabilistic sense,optimal and unique approximations to these events in B Λ ,N . Lemma 5.1. If B, C ∈ B Λ ,N and µ Λ ( B △ C ) > , then there exists at least one λ ( k ) , k ∈{ , , . . . , N } ), such that λ ( k ) is in one and only one of the measurable sets B or C .Proof: By definition, the σ − algebra B Λ ,N is constructed by the closure of complements andcountable unions of the finite set of individual Voronoi cells {V j } Nj =1 . Thus, by construction,any sets of non-zero µ Λ -measure in B Λ ,N must contain at least the interior of a single Voronoicell. Since B △ C ∈ B Λ ,N , the conclusion follows. (cid:3) Theorem 5.1.
Suppose (cid:8) λ ( j ) (cid:9) Nj =1 is fixed and choose any A ∈ B Λ such that there exists atleast one λ ( j ) ∈ A . There exists unique (up to a set of zero µ Λ -measure) B ∈ B Λ ,N such thatfor any P D (absolutely continuous with respect to µ D ), for any simple function approximationto ρ D , and for any choice of Ansatz, we have ˜ P Λ ,N ( A ) = ˜ P Λ ,N ( B ) .Proof: Fix an arbitrary A ∈ B Λ containing at least one sample and let J denote the setof indices such that λ ( j ) ∈ A for j ∈ J and B := ∪ j ∈J V j . By construction, B is a Voronoicoverage of A and by definition of counting measure ˜ P Λ ,N , for any probability measure P D absolutely continuous with respect to µ D and any simple function approximation to ρ D , wehave that ˜ P Λ ,N ( A ) = ˜ P Λ ,N ( B ). Lemma 5.1 gives uniqueness (see the Appendix for details). (cid:3) Theorem 5.1 guarantees that for any event in the original σ -algebra B Λ of which we wantto compute its probability that there is a “best” approximation in the numerical σ -algebra Borel σ − algebras on R are commonly generated using generating sets defined, for example, by all half-openintervals. Analogously, B Λ ,N is equivalently generated by the set of all Voronoi coverages for all A ∈ B Λ . Thisparticular definition uses the minimal generating set. B Λ ,N as long as the original event in B Λ contains at least one sample in Λ. While we useprobabilities to prove this, the probabilities on D are independent of the choices of sampleson Λ (the probability measure on outputs exists independently of our choice of approximatingevents in the parameter domain). In other words, Theorem 5.1 is a statement about set-approximation not the approximation of probabilities since ˜ P Λ ,N is exact.The set-approximation error can be determined a priori to any computation of the model(and bounded in µ Λ -measure independently of any probability measure) either globally orfor a certain collection of events of physical importance or interest in B Λ . For the sake ofsimplicity and also for practical reasons, we consider any event of interest to belong to B Λ ,N ,i.e. in Eq. (5.1) we set ˜ B Λ = B Λ ,N . In other words, we assume that all the events we wishto compute the probabilities of can be represented exactly (or with negligible error) in B Λ ,N .This prevents the problem of trying to compute probabilities of complicated events that existeven in standard Borel σ − algebras, which are usually not of interest in a realistic settingwhere events are often defined by hypercubes or balls. Moreover, all such events A ∈ B Λ ,N have the property that µ Λ ( ∂A ) = 0.Recall that { I i } Mi =1 denotes the partition on D used to define the simple-function approx-imation to the density of P D as ρ D ,M = M X i =1 P D ( I i ) I i ( q ) . Suppose that every A i := Q − ( I i ) ∈ B Λ ,N , then it follows that the error given by term I is zero for all A ∈ B Λ ,N . However, it is unlikely that the induced regions of generalized contours A i are exactly represented in B Λ ,N and this is the source of error in term I of Eq. (5.1). Theexact density of P Λ given the above density ρ D ,M can be written as ρ Λ = M X i =1 P D ( I i ) A i ( λ ) , (5.2)and the exact density of P Λ ,N given ρ D ,M can be written as˜ ρ Λ ,N = M X i =1 P D ( I i ) A i,N ( λ ) , (5.3)where A i,N is the unique Voronoi coverage of A i guaranteed to exist by Theorem 5.1 for each i = 1 , , . . . , M . Thus, for any A ∈ B Λ ,N , we have that P Λ ( A ) − P Λ ,N ( A ) = M X i =1 P D ( I i ) (cid:18) µ Λ ( A ∩ A i ) µ Λ ( A i ) − E ( A, A i ) (cid:19) . (5.4)Here E ( A, A i ) is a computable constant determined by the Voronoi cells defining the uniquecovers of A and A i,N . If we use the standard MC implementation of Algorithm 2, then E ( A, A i ) is equal to the number of samples in A ∩ A i over the number of samples in A i . If wedo not use the standard MC approximation or use non-uniform sampling so that the variant of igma-Algebras on Contour Maps 17 Algorithm 2 discussed in Section 4.1 is used, then E ( A, A i ) is equal to µ Λ ( A ∩ A i,N ) /µ Λ ( A i,N ).We also consider these options in deriving computable estimates of term II in Section 5.2below. The error representation of Eq. 5.4 is uncomputable as it requires knowledge of theexact volume of A i .To derive the computable bounds for term I described below, we require the following Assumption 5.1.
Assume that for a fixed approximation to P D that a sufficient number ofsamples in Λ are taken so that there is at least one sample λ ( j ) in A i,N for each i such thatthe associated Voronoi cell V j shares no boundary with a Voronoi cell from any other regionof approximate induced generalized contours A k,N , k = i . Theorem 5.2.
If Assumption 5.1 holds, then there exists signed computable lower and uppera posteriori bounds on term I of Eq. 5.1 for every A ∈ B Λ ,N .Proof: Let B i,N denote the union of Voronoi cells belonging to A i,N that share no boundarywith A k,N for all i , k = i . By Assumption 5.1, B i,N = ∅ and µ Λ ( B i,N ) > i . Let C i,N denote the union of A k,N and all Voronoi cells sharing boundary with A k,N . The sets B i,N and C i,N are identifiable (e.g. using nearest neighbor searches) and computable since they belongto B Λ ,N for all i . Furthermore, we have that µ Λ ( A ∩ B i,N ) µ Λ ( C i,N ) ≤ µ Λ ( A ∩ A i ) µ Λ ( A i ) ≤ µ Λ ( A ∩ C i,N ) µ Λ ( B i,N ) . (5.5)Substitution of Eq. (5.5) into Eq. (5.4) gives the following signed computable lower andupper a posteriori bounds on I , I ≤ M X i =1 P D ( I i ) max (cid:26)(cid:18) µ Λ ( A ∩ B i,N ) µ Λ ( C i,N ) − E ( A, A i ) (cid:19) , (cid:18) µ Λ ( A ∩ C i,N ) µ Λ ( B i,N ) − E ( A, A i ) (cid:19)(cid:27) , (5.6) I ≥ M X i =1 P D ( I i ) min (cid:26)(cid:18) µ Λ ( A ∩ B i,N ) µ Λ ( C i,N ) − E ( A, A i ) (cid:19) , (cid:18) µ Λ ( A ∩ C i,N ) µ Λ ( B i,N ) − E ( A, A i ) (cid:19)(cid:27) . (5.7)This completes the proof (cid:3) . As described above inSection 5.1, we assume ˜ B Λ = B Λ ,N in Eq. (5.1). By assumption of local differentiability ofmap Q , we assume the numerical method used to compute Q h produces a piecewise continuousapproximation to Q and there exists a piecewise continuous error function e such that Q ( λ ) = Q h ( λ ) + e ( λ ) , ∀ λ ∈ Λ . We assume there exists computable a posteriori error estimates for e denoted by e h such that Q ( λ ) ≈ Q h ( λ ) + e h ( λ ) , ∀ λ ∈ Λ . It is not uncommon for a sufficient number of samples to be assumed in standard inequalities for empiricaldistribution functions, e.g. forms of the multivariate Dvoretzky-Kiefer-Wolfowitz inequality require the squareroot of the number of samples to be greater than the dimension of the random variable divided by the desired α -level [12, 18]. Analogous to Algorithm 1, let A i,N and A i,N,h denote the Voronoi coverages of Q − ( I i ) and Q − h ( I i ), respectively. Following Algorithm 2, we have that, for any A ∈ B Λ ,N , term II inEq. (5.1) can be written as II = M X i =1 P D ( I i ) (cid:8) j : λ ( j ) ∈ A i,N,h ∩ A (cid:9)(cid:8) j : λ ( j ) ∈ A i,N,h (cid:9) − (cid:8) j : λ ( j ) ∈ A i,N ∩ A (cid:9)(cid:8) j : λ ( j ) ∈ A i,N (cid:9) ! . (5.8)Let J i,A , and J i , denote (cid:8) j : λ ( j ) ∈ A i,N,h ∩ A (cid:9) and (cid:8) j : λ ( j ) ∈ A i,N,h (cid:9) , respectively. Notethat the numbers J i,A and J i are computable. Terms involving A i,N are not computable sincethey require the exact map Q . We may invert Q h + e h to more accurately estimate A i,N , anddefine the computable approximations n j : λ ( j ) ∈ A i,N ∩ A o ≈ n j : λ ( j ) ∈ A, and Q h ( λ ( j ) ) + e h ( λ ( j ) ) ∈ I i o , (5.9) n j : λ ( j ) ∈ A i,N o ≈ n j : Q h ( λ ( j ) ) + e h ( λ ( j ) ) ∈ I i o . (5.10)Let J i,A,e and J i,e denote the computable approximations of Eqs. (5.9) and (5.10), respectively.Substituting these approximations into Eq. (5.8) and re-arranging terms gives the computablea posteriori estimate II ≈ M X i =1 P D ( I i ) (cid:18) J i,A J i,e − J i,A,e J i J i J i,e (cid:19) (5.11)As discussed in Section 4.1, if we replace the counting vector in Algorithm 2 with estimatesof volumes for the Voronoi cells, then the above summations of number of cells are replacedby the summations of the volumes of the associated Voronoi cells. For example, in this casewe set J i,A to be the number determined by the following X { j : λ ( j ) ∈ A i,N,h ∩ A } µ Λ ( V j ) . Computable error estimates and/orerror bounds are useful in adaptive algorithms. This is beyond the scope of this paper, but isthe subject of future work on adaptive sampling strategies for computing accurate probabilitiesof specified events, e.g. rare events.In the numerics below, we consider an alternative use of the a posteriori error estimatesof the map Q h . Specifically, we use the a posteriori error estimates to improve the functionalevaluations. This is motivated by previous work on improved linear functionals where aposteriori error estimates are used to improve the pointwise accuracy of the QoI map [3] andimprove the pointwise accuracy of probability distributions propagated through likelihoodsinvolving a QoI map [4].To show the effect of using a posteriori error estimates to improve the accuracy of computedprobabilities of a counting measure we first compute a so-called “reference solution” that isan accurate approximation to P Λ on a fine grid of uniformly spaced rectangles partitioningΛ. When computing the reference solution, we use a posteriori error estimates on the finegrid to correct for the deterministic error of the map. We follow this by the random sampling igma-Algebras on Contour Maps 19 algorithm for computing counting measures both with and without the a posteriori errorestimate correcting for the deterministic error. Finally, we demonstrate that the errors of thecounting measures with respect to the reference solution evaluated on generating events fora given ˜ B Λ are reduced globally when the a posteriori error estimate is used to improve thepointwise accuracy of Q h .
6. Numerical Results.6.1. Improving probabilities with a posteriori error estimates.
We take as the modelthe Navier-Stokes equation ∂ u ∂t − ν ∇ · ∇ u + u · ∇ u + ∇ p = g u , x ∈ Ω , t ∈ ( t , t f ) , (6.1)where Ω is defined on a rectangular domain with a cylindrical hole as shown in the left plotof Figure 6.1. We describe the boundary conditions as specified for a similar problem in[10] with fixed parameter values. For the rectangular boundary, the inflow velocity is set to(1 , u · n = 0) are set along the top and bottomboundaries, and a natural outflow condition is set on the right boundary. We use no-slipboundary conditions ( u = 0) on the cylindrical boundary. The diameter of the cylinder is 1 / ν ) − . We focus on steady-statesolutions, so we consider the situation where the fluid viscosity is uncertain but known to bewithin [0 . , .
1] resulting in Reynolds numbers bounded by 50. Thus, the parameters forthis problem are the viscosity and the vertical displacement of the center of the cylinder. Let λ ∈ Λ = [0 . , . × [ − . , . ⊂ R denote the parameter domain with viscosity the firstcomponent. −2 −1 0 1 2 3 4 5 6−2−1.5−1−0.500.511.52 Figure 6.1.
Left: Physical domain Ω in Eq. 6.1. Right: The QoI map Q ( λ ) is the first component of velocityat the point (1 , . in Ω . The viscosity ν is taken as λ in the parameter domain. The vertical displacementof the center of the cylindrical hole is taken as λ in the parameter domain. For the numerical solution of the model, we use SUPG-PSPG-LSIC stabilization withpiecewise linear finite elements for the forward problem. The QoI we consider is the firstcomponent of the velocity at the point (1 , .
2) which is just behind the cylinder and off thecenter-line. This QoI defines a differentiable map shown on the right plot of Figure 6.1. Toestimate the error in the numerically evaluated QoI, we use an a posteriori error estimate using solution to an adjoint problem. We solve the adjoint problem using SUPG-PSPG-LSICstabilization with piecewise quadratic finite elements. This approach has been demonstratedto produce reliably accurate error estimates for a variety of QoI in [10]. −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.700.511.522.5 Q( λ ) Figure 6.2.
Left Plot: Plots of the Beta (4 , distributions centered at Q (0 . , . (solid line withleftmost peak), Q (0 . , . (dashed line with middle peak), and Q (0 . , − . (dotted line withrightmost peak). Remaining plots show approximations of P Λ ( R j ) on × rectangles { R j } partitioning Λ sothat moving left-to-right the associated output densities have means at Q (0 . , . , Q (0 . , . , and Q (0 . , − . , respectively. The parameter values (0 . , . , (0 . , . , and (0 . , − . are denoted by white boxes in these plots. We create an approximate reference solution to P Λ in the following way. First, we use62,500 samples taken from a regularly spaced 250 ×
250 set of points in Λ. For each sampledparameter, we compute both the QoI and the associated a posteriori error estimate. Since D iscompact, we assume the output density is approximated parametrically by a Beta distributionBeta( α, β ) within the range of the QoI map. The parameters α and β are chosen so the mean ofthe distribution is at a specified “exact” observed value and the distribution “appears Normal”in shape . For the cases shown here, we set α = 4 and β = 5 and took the mean QoI valuesfrom a more accurate (i.e. a posteriori error corrected) mapping of known parameter values,see the leftmost plot of Figure 6.2. We used three separate parameters within Λ that mapto nearby QoI values to demonstrate the sensitivity of inverse probability measures to smallchanges in the location of high probability intervals in D . We approximate each ρ D witha simple function ρ D,M defined on a grid of 200 uniformly spaced points defining 199 binspartitioning D and binning 1 E + 6 i.i.d. samples from the associated Beta distributions. Wethen follow the steps of Algorithm 1 using the 62,500 samples of Λ evaluated on the errorcorrected map. In Figures 6.2 and 6.3, we show the resulting approximate plots of P Λ ( R j )computed on 50 ×
50 uniformly sized rectangles { R j } partitioning Λ. In other words, we usethe 250 ×
250 partition to define P Λ , ≈ P Λ and create plots of probabilities on a 50 ×
50 gridof rectangles, which gives a good visual approximation to the shape of the density function.Below, we generally refer to any computation involving P Λ , as a reference solution.We observe in the plots of the reference solution in Figures 6.2 and 6.3 that in two cases theassociated densities will exhibit a type of “bi-modality” across two distinct regions of inducedgeneralized contours that appear to fan out as viscosity (i.e. the λ parameter) is increased. We could also simply truncate Normal distributions and renormalize. We simply avoid this nuisance here. igma-Algebras on Contour Maps 21Figure 6.3.
Left: to Right: A 3-D view of approximations of P Λ ( R j ) on × rectan-gles { R j } partitioning Λ so that moving left-to-right the associated output Beta densities have meansat Q (0 . , . , Q (0 . , . , and Q (0 . , − . , respectively. The parameter values (0 . , . , (0 . , . , and (0 . , − . are denoted by white boxes in these plots. This is due in part to the fact that the peaks of the associated output Beta distributions withmeans at Q (0 . , . Q (0 . , − . ,
5) distribution withmean at Q (0 . , − . B Λ as the σ − algebra generated from a 10 ×
10 uniform rectangular grid of Λ,i.e. we assume any event we want to measure in probability can be defined by unions of asubset of 100 uniformly sized rectangles in Λ. We plot the reference solution computed on thegenerating events of ˜ B Λ in the left plots of Figures 6.4 and 6.5. λ λ λ λ λ λ Figure 6.4.
Reference solution (left), ˜ P Λ , k,h (middle), and ˜ P Λ , k,h (right) evaluated on the 100 gener-ating events of ˜ B Λ We use 15,000 i.i.d. uniform random samples in Λ to compute two separate countingmeasures denoted by ˜ P Λ , k,h and ˜ P Λ , k,h . We compute ˜ P Λ , k,h using Q h to evaluate thesamples in Λ and compute ˜ P Λ , k,h using Q h + e h (i.e. using computed a posteriori error estimates to correct for error in the numerically evaluated QoI) to evaluate the same set ofsamples in Λ. We plot ˜ P Λ , k,h and ˜ P Λ , k,h given the generating events of ˜ B Λ in the middleand right plots, respectively, of Figures 6.4 and 6.5. λ λ λ λ λ λ Figure 6.5. ˜ P Λ , k,h (middle), and ˜ P Λ , k,h (right) evaluated onthe 100 generating events of ˜ B Λ In Figures 6.6 and 6.7, we show the signed errors ˜ P Λ , k,h − P Λ , and ˜ P Λ , k,h − P Λ , ,respectively, on the generating sets of ˜ B Λ . We observe a general reduction in the magnitudeof the errors across large areas of the domain when we compute the counting measure usingthe a posteriori error estimates to correct for the numerical error. λ λ −4−20246x 10 −3 −3 λ λ Figure 6.6.
Plots of ˜ P Λ , k,h − P Λ , on the generating sets of ˜ B Λ .
7. A higher dimensional example.
We apply Algorithm 2 to compute a counting measurefor a 15-dimensional parameter space defined by uncertain coefficients and initial conditions ofan MSEIRS model. The MSEIRS model is a generalization of the well studied SIR epidemicmodel that takes into account groups of infants protected by maternal antibodies (the “M”class) and groups of exposed and latent infected but not infectious (the “E” class) [16]. The igma-Algebras on Contour Maps 23 λ λ −4−20246x 10 −3 −3 λ λ Figure 6.7.
Plots of ˜ P Λ , k,h − P Λ , on the generating sets of ˜ B Λ . coupled nonlinear system of differential equations defining this model is given by dMdt = B ( S + E + I + R ) − ( δ + µ M ) M dSdt = δM − βSI − ( µ G + ι ) S + f R dEdt = βSI − ( ǫ + µ G ) E dIdt = ǫE − ( γ + µ I + µ G ) I dRdt = γI − ( µ G + f ) R + ιS (7.1)We define the uncertain parameters of this model in Table 7.1 including the bounds (shownwithout dimension and a single unit of time is 1 week, a unit of population is 1 E
6, andbirth/death rates are normalized to a population of size 300 million) that we use to define Λ.We sample the 15-dimensional Λ uniformly 1 E Q and Q , respectively). These quantities are uncertain (e.g., they may be estimated by surveys ofhospital data during the disease outbreak), and it is generally impossible to know these valuesexactly at any particular time during the outbreak. Here, we consider the problem of param-eter identification under uncertainty. We seek to determine the most probable configurations(defined by events) of parameters and initial conditions associated with the uncertain outputdata.We use a tensor product of marginal probability densities defined by shifted and scaled(Beta) B (4 ,
5) distributions to model the uncertainty in the joint QoI (denoted by Q =( Q , Q )). These densities are defined by first solving the model with the reference parameterand initial conditions shown in Table 7.1, which produced reference QoI values of approxi-mately Q = 1 .
54 and Q = 2 .
43 (in millions). We determined the rescaling and shifting ofthe Beta distributions so that the range of any sampled QoI value stays within ± .
15 millionpeople of the reference QoI value. We first invert only the marginal densities and then invert
B avg. birth rate [2 . E − , . E −
4] 3 . E − δ avg. temporary immunity [1 / , /
4] 0.16 µ M avg. infant death rate [4 E − , E − . E − β contact/infectivity rate [1 . E − , . E −
3] 3 . E − µ G avg. general death rate [2 . E − , . E −
4] 2 . E − /ǫ avg. infection time [0 . ,
1] 0.7 µ I avg. infected death rate [4 . E − , . E −
5] 1 . E − /γ avg. recovery time [0 . , .
33] 0 . f avg. loss of immunity rate [0 . , .
25] 0.18 ι avg. immunization rate [0 . , . M initial infants [2 . , .
5] 3.25 S initial susceptibles [260 , E initial exposed [0 . , .
5] 0.425 I initial infected [0 . ,
4] 3.8 R initial recovered/immunized [10 ,
20] 13
Table 7.1
Uncertain parameters and initial conditions and their interval bounds in the MSEIRS model given byEq. 7.1. The last column of values are the reference parameter values used to define the mean of the outputdensities. the joint density. To invert the marginal densities using Algorithm 2, we approximate eachmarginal using simple function approximations computed on a uniform grid of 200 subinter-vals on the computed ranges of Q and Q . To invert the marginal density, we use 50 × Q and Q tocompute the simple function approximation.In high dimensional parameter spaces, visualization and analysis of probability measuresand/or densities can be difficult. Here, to highlight some of the underlying geometric structure,we show plots of some of the more interesting marginals for which Q and Q exhibit differentsensitivities to the underlying parameters in Figure 7.1. Most of the marginals over pairsof parameters not including those shown in Figure 7.1 were approximately uniform due tothe lack of sensitivity of the QoI maps with respect to the parameters (i.e., the generalizedcontours were nearly flat in the directions of these parameters). We see that inverting thedensity for either Q , Q , or Q all determine different marginal distributions on ( γ, δ ). Ingeneral, each QoI can only produce “good estimates” of two or three parameters in terms ofa particular reference parameter or initial condition being in a relatively small event of highprobability. Using both QoI provides good estimates of five of the 15 parameters and initialconditions. By utilizing additional geometrically distinct QoI, we can further improve resultsin terms of localizing the probability into regions of decreasing volume as described in [7]. Appendix A. Proof of Uniqueness in Theorem 5.1.
Suppose there is also a C ∈ B Λ ,N satisfying the equality for arbitrary output probabilitymeasures and µ Λ ( B △ C ) >
0. By Lemma 5.1, there exists a λ ( k ) for some k ∈ { , . . . , N } ineither B or C but not both. Without loss of generality, let K C \ B denote the set of indicessuch that λ ( k ) ∈ C \ B for all k ∈ K C \ B . Let I K C \ B ⊂ D denote the measurable event defined igma-Algebras on Contour Maps 25 γ δ γ δ −3 γ β β I γ δ Figure 7.1.
Thin plate spline approximations to some of the counting measure densities for various inverseproblems. The top left and middle plots are marginals over ( M , δ ) and ( γ, δ ) , respectively, computed frominverting the Beta density on Q . The top right plot is a marginal over ( γ, δ ) computed from inverting thedensity on Q . The bottom left and middle plots are marginals over ( γ, β ) and ( γ, I ) , respectively, computedfrom inverting the Beta density on Q . The bottom right plot is a marginal over ( γ, δ ) computed from invertingthe joint output density on Q . by I K C \ B = Q ( ∪ k ∈K C \ B V ( λ ( k ) )) ∈ B D and P D be defined by a uniform probability measure on I K C \ B . This implies ρ D is a simple function and we choose ρ D ,M as this exact density. Thereare two cases to consider: (1) µ Λ ( B \ C ) = 0, or (2) µ Λ ( B \ C ) = 0. For case (1), we immediatelyhave that ˜ P Λ ,N ( C ) = ˜ P Λ ,N ( B ), which is a contradiction. For case (2), let K B \ C denote theset of indices such that λ ( k ) ∈ B \ C for any k ∈ K B \ C . If Q ( λ ( k ) ) / ∈ I K C \ B for all k ∈ K B \ C ,then we immediately arrive at the contradiction ˜ P Λ ,N ( C ) = ˜ P Λ ,N ( B ) by construction of ˜ P Λ ,N .If there exists any k ∈ K B \ C such that Q ( λ ( k ) ) ∈ I K C \ B , then the corresponding Voronoicells are identified as approximating induced regions of generalized contours with non-zeroprobabilities whose global approximation also includes Voronoi cells indexed by K C \ B . Weare free to choose any Ansatz with the required modifications to Algorithm 2 and choosingthe Ansatz such that probabilities of parts of generalized contours through B not included in C are set to zero leads to a contradiction. (cid:3) REFERENCES [1]
V.I. Bogachev , Measure Theory (Volume 2) , Springer, 2007.[2]
J. Breidt, T. Butler, and D. Estep , A Measure-Theoretic Computational Method for Inverse Sensitiv-ity Problems I: Method and Analysis , SIAM Journal on Numerical Analysis, 49 (2011), pp. 1836–1859.[3]
T. Butler, P. Constantine, and T. Wildey , A posteriori error analysis of parameterized linearsystems using spectral methods , SIAM. J. Matrix Anal. Appl., 33 (2012), pp. 195–209.[4]
T. Butler, C. Dawson, and T. Wildey , Propagation of uncertainties using improved surrogate models ,SIAM/ASA Journal on Uncertainty Quantification, 1 (2013), pp. 164–191. [5]
T. Butler and D. Estep , A numerical method for solving a stochastic inverse problem for parameters ,Annals of Nuclear Energy, 52 (2013), pp. 86–94.[6]
T. Butler, D. Estep, and J. Sandelin , A Computational Measure Theoretic Approach to InverseSensitivity Problems II: A Posteriori Error Analysis , SIAM Journal on Numerical Analysis, 50 (2012),pp. 22–45.[7]
T. Butler, D. Estep, S. Tavener, C. Dawson, and J.J. Westerink , A Measure-Theoretic Compu-tational Method For Inverse Sensitivity Problems III: Multiple Quantities of Interest , SIAM Journalon Uncertainty Quantification, (2014), pp. 1–27. In press.[8]
T. Butler, L. Graham, C. Dawson, D. Estep, and J.J. Westerink , Quantifying uncertainty ofland classification within the advanced circulation (adcirc) model I: A computational framework . inpreparation.[9]
J.T. Change and D. Pollard , Conditioning as disintegration , Statistica Neerlandica, 51 (1997),pp. 287–317.[10]
E. Cyr, J. Shadid, and T. Wildey , Approaches for adjoint-based a posteriori analysis of stabilizedfinite element methods , SIM Journal on Scientific Computing, 36 (2014), pp. A766–A791.[11]
C. Dellacherie and P.A. Meyer , Probabilities and Potential , North-Holland Publishing Co., Amster-dam, 1978.[12]
L.P. Devroye , A uniform bound for the deviation of empirical distribution functions , Journal of Multi-variate Analysis, 7 (1977), pp. 594 – 597.[13]
J.R. Dormand and P.J. Prince , A family of embedded runge-kutta formulae , Journal of Computationaland Applied Mathematics, 6 (1980), pp. 19–26.[14]
Gerald B. Folland , Real Analysis: Modern Techniques and Their Applications , Wiley, 1999.[15]
A. Friedman , Foundations of Modern Analysis , Dover Publiciations, 1982.[16]
H.W. Hethcote , The mathematics of infectious diseases , SIAM Review, 42 (2000), pp. 599–653.[17]
E. Khmaladze and N. Toronjadze , On the almost sure coverage property of voronoi tessellation: Ther1 case , Advances in Applied Probability, 33 (2001), pp. 756–764.[18]
P. Massart , The tight constant in the dvoretzky-kiefer-wolfowitz inequality , The Annals of Probability,18 (1990), pp. 1269–1283.[19]