Techniques for generic probabilistic inversion

Abstract

Probabilistic inversion problems are defined, existing algorithms for solving such problems are discussed, and new algorithms based on iterative re-weighting of a sample are introduced. One of these is the well-known iterative proportional fitting whose properties were already studied (Csiszar, Ann. Probab. (3) (1975) 146). A variant on this is shown to have fixed points minimizing an information functional, even if the problem is not feasible, and is shown to have only feasible fixed points if the problem is feasible. The algorithm is not shown to converge, but the relative information of successive iterates is shown to converge to zero. Applications to atmospheric dispersion and environmental transport are discussed.

Introduction

Probabilistic inversion problems can be formulated either from a measure theoretic or from a random variable viewpoint. The former may be more suggestive for theoreticians and the latter is more recognizable for practitioners.

We follow the approach of Kraan (2002) and Kraan and Bedford (2004). Let $(M,\mathcal{B}(M),\lambda )$ and $(N,\mathcal{B}(N),\nu )$ be two Borel probability spaces, where $M,N$ are compact non-empty subsets of ${\mathbb{R}}^m$ and ${\mathbb{R}}^n$, respectively. Let $T:{\mathbb{R}}^m\mapsto {\mathbb{R}}^n$ be a continuous mapping. $\gamma =\lambda \circ T^{-1}$ is called the push forward of $\lambda$ under T, and similarly, $\lambda$ is called the pull back of $\gamma$ under T. In the following $\lambda$ plays the role of a background measure. Measure $\mu$ on $(M,\mathcal{B}(M))$ is called an inverse of T at $\nu$ if $\mu$ is the pull back of $\nu$ under T; that is

$$\forall B\in \mathcal{B}(N),\mu \circ T^{-1}(B)=\nu (B).$$

The problem can then be formulated as

Definition 1 Probabilistic inversion problem

Given $(M,\mathcal{B}(M),\lambda )$ and $(N,\mathcal{B}(N),\nu )$, with $T:{\mathbb{R}}^m\mapsto {\mathbb{R}}^n$ continuous, find a measure $\mu$ absolutely continuous with respect to $\lambda$, $\mu ⪡\lambda$, on $(M,\mathcal{B}(M))$ such that $\mu$ is an inverse of T at $\nu$.

Such problems may be infeasible, and if feasible may have many solutions (Kraan, 2002). Under certain conditions a measure $\mu$ can be found solving (1). If $\nu ⪡\gamma$ the Radon–Nikodym derivative $g:N\mapsto \mathbb{R}$ exists, is non-negative, unique up to sets of $\gamma$-measure zero, and satisfies

$$\forall B\in \mathcal{B}(N),\nu (B)={\int }_{B}g(y)\mathrm{d}\gamma (y).$$

If in addition g is continuous, $g(y)>0$, $y\in N$, and ${\int }_{N}g(y)\log g(y)\mathrm{d}\gamma (y)<\infty$, then define $\tilde{f}≔g\circ T$, and define $\tilde{\mu }$ as a new probability measure:

$$\forall A\in \mathcal{B}(\mathcal{M}),\tilde{\mu }≔{\int }_{T^{-1}(B)}\tilde{f}\mathrm{d}\lambda (x).$$

It is easy to check that $\tilde{\mu }$ is an inverse of T at $\nu$:

$$\tilde{\mu }\circ T^{-1}(B)={\int }_{T^{-1}(B)}\tilde{f}\mathrm{d}\lambda (x)$$

$$={\int }_{T^{-1}(B)}g\circ T(x)\mathrm{d}\lambda (x)={\int }_{B}g(y)\mathrm{d}\gamma (y)=\nu (B).$$

It can be shown that the measure $\tilde{\mu }$ is the unique measure satisfying Eq. (1) and minimizing the relative information with respect to $\lambda$ in the class of measures satisfying Eq. (1) (Kraan, 2002).

In practical problems a reformulation of the problem in terms of random variables may be more convenient. Given a random vector $\mathbf{Y}$ taking values in ${\mathbb{R}}^M$ and a measurable function $G:{\mathbb{R}}^N\to {\mathbb{R}}^M$, find a random vector $\mathbf{X}$ such that $G(\mathbf{X})\sim \mathbf{Y}$, where $\sim$ means that $G(\mathbf{X})$ and $\mathbf{Y}$ share the same distribution. If $G(\mathbf{X})\in \{\mathbf{Y}|\mathbf{Y}\in \mathcal{C}\}$ where $\mathcal{C}$ is a subset of random vectors on ${\mathbb{R}}^M$, then $\mathbf{X}$ is called a probabilistic inverse of G at $\mathcal{C}$. $\mathbf{X}$ is sometimes termed the input to model G, and $\mathbf{Y}$ the output. G corresponds to the mapping T in the above paragraph. Note that T was assumed continuous, whereas G is only assumed to be measurable. If the problem is feasible it may have many solutions and we require a preferred solution; if it is infeasible we seek a random vector $\mathbf{X}$ for which $G(\mathbf{X})$ is ‘as close as possible’ to $\mathbf{Y}$. Usually as a measure of closeness the relative information is used (Kullback, 1959).

Probabilistic inversion problems arise in quantifying uncertainty in physical models with expert judgement (Kraan and Cooke, 2000a). We wish to quantify the uncertainty on parameters $\mathbf{X}$ of some model using expert judgement, but the parameters do not possess a clear physical meaning and are not associated with physical measurements with which experts are familiar. Often the models are derived under assumptions to which the experts do not subscribe. We must then find the observable quantities $\mathbf{Y}$ functionally related with $\mathbf{X}$ that can be assessed by experts. Extracting uncertainties of $\mathbf{X}$ from uncertainties of $\mathbf{Y}$ specified by experts is clearly an inverse problem (see examples in Section 5).

In practical applications the random vector $\mathbf{Y}$ is characterized in terms of some percentiles or quantiles of the marginal distributions $Y_1,\dots ,Y_M$. In this case, we seek a random vector $\mathbf{X}$ such that $G(\mathbf{X})$ satisfies quantile constraints imposed on $Y_1,\dots ,Y_M$. There may be other constraints. These constraints may reflect mathematical desiderata, as when we require independence between variables in Section 5. Physical considerations may also impose constraints on $\mathbf{X}$. In some cases probabilistic inversion problems may have trivial solutions, e.g. if the $\mathbf{X}$ makes the coordinates of $\mathbf{Y}$ completely rank correlated. Such solutions may be rejected on physical grounds; hence physical constraints may stipulate the support of $\mathbf{X}$. Other physical constraints are discussed in the example in Section 6. A few algorithms for solving such problems are available in literature namely: conditional sampling, parameter fitting for uncertain models (PARFUM) (Cooke, 1994), Hora and Young algorithm (Harper et al., 1994) and PREJUDICE (Kraan and Cooke, 2000a, Kraan and Cooke, 2000b). We summarize existing approaches to probabilistic inversion in Section 2.

In this paper we study iterative algorithms for numerically solving probabilistic inversion problems. These methods do not require model inversion they are based on sample re-weighting techniques and are promising because they do not require special knowledge about the problem at hand, or complicated heuristic steering on the part of the user. Moreover, operations on the sample are performed one-at-a-time, so the entire sample need not be held in memory. This means that there is virtually no problem size limitation.

After reviewing existing approaches, we discuss two iterative algorithms. First is known in the literature as iterative proportional fitting (IPF) (Kruithof, 1937). IPF finds a joint distribution satisfying marginal constraints by successively I-projecting an initial distribution on the set of distributions satisfying each marginal constraint. The I-projection of a probability measure p on a set of measures Q is ${\mathrm{argmin}}_{q\in Q}I(q|p)$, where $I(q|p)$ is the relative information of q with respect to p. IPF need not converge, but if it converges, it converges to a solution which is minimally informative with respect to the starting distribution (Csiszar, 1975). A variation on this is an iterative version of the PARFUM algorithm. We show that this algorithm has fixed points minimizing an information functional, even if the problem is not feasible, and that it has only feasible fixed points if the problem is feasible. The algorithm is not shown to converge, but the relative information of successive iterates is shown to converge to zero.

In Section 4 we discuss a sample re-weighting iterative approach to probabilistic inversion. In Section 5 we illustrate the IPF and PARFUM algorithms with an example from atmospheric dispersion modelling (Kraan and Cooke, 2000a) and transfer coefficients in chicken processing line. Iterative algorithms can easily be adopted to satisfy joint as well as marginal constraints on $Y_1,\dots ,Y_M$. We illustrate this by imposing (approximate) convolution constraints in Section 6. The final section gathers conclusions.