# “Stochastic Inverse Problems” and Changes-of-Variables ††thanks: This work
was supported in part by the NNSA Office of Defense Nuclear Nonproliferation
Research and Development, NA-22 F2019 Nuclear Forensics Venture
(LA19-V-FY2019-NDD3Ac).

Peter W. Marcy Los Alamos National Laboratory, Los Alamos, NM
<EMAIL_ADDRESS>Rebecca E. Morrison University of Colorado Boulder, Boulder,
CO<EMAIL_ADDRESS>
###### Abstract

Over the last decade, a series of applied mathematics papers have explored a
type of inverse problem—called by a variety of names including “inverse
sensitivity”, “pushforward based inference”, “consistent Bayesian inference”,
or “data-consistent inversion”—wherein a solution is a probability density
whose pushforward takes a given form. The formulation of such a _stochastic_
inverse problem can be unexpected or confusing to those familiar with
traditional Bayesian or otherwise statistical inference. To date, two classes
of solutions have been proposed, and these have only been justified through
applications of measure theory and its disintegration theorem. In this work we
show that, under mild assumptions, the formulation of and solution to _all_
stochastic inverse problems can be more clearly understood using basic
probability theory: a stochastic inverse problem is simply a change-of-
variables or approximation thereof. For the two existing classes of solutions,
we derive the relationship to change(s)-of-variables and illustrate using
analytic examples where none had previously existed. Our derivations use
neither Bayes’ theorem nor the disintegration theorem explicitly. Our final
contribution is a careful comparison of changes-of-variables to more
traditional statistical inference. While taking stochastic inverse problems at
face value for the majority of the paper, our final comparative discussion
gives a critique of the framework.

Keywords: Jacobian, reparameterization, uncertainty quantification,
statistical inference, Bayesian analysis

##  1 Introduction

Scientific models often relate unknown parameters or functions to observable
quantities. Any choice of unknown input to the model will produce an
observable that can be compared with real measurements. In the absence of
noise or measurement error, this _forward_ problem can be thought of as a
well-defined operator from parameter/function-space to observable-space. To
solve the _inverse_ problem (IP) is to estimate the unknown inputs from a
finite collection of observations. In the more formal language of Tikhonov et
al. (1995), if $\mathcal{F}$ is an operator between metric (or Banach, or
Hilbert) spaces representing parameters and data,
$\mathcal{F}:\mathcal{P}\mapsto\mathcal{Q}$, the solution to the IP is defined
by the operator equation

$\displaystyle\boldsymbol{\theta}$ $\displaystyle\text{\quad such that
\quad}\mathcal{F}(\boldsymbol{\theta})=\boldsymbol{y}\quad\quad\big{(}\boldsymbol{\theta}\in\mathcal{P},\
\boldsymbol{y}\in\mathcal{Q}\big{)}\ .$ (1)

This operator formulation is typically associated with function spaces, where
solving the IP means solving equations (often integral equations) given a
finite number of potentially noisy observations (Groetsch, 1993; Stuart, 2010;
Kirsch, 2011; Aster et al., 2019; Lesnic, 2021). Equality within (1) may not
be possible, and a solution might instead be required to minimize some
functional involving $\mathcal{F}(\boldsymbol{\theta})$ and $\boldsymbol{y}$,
as in the case of weighted or regularized least-squares estimation.

When $\mathcal{P}$ and $\mathcal{Q}$ are subsets of $\mathbb{R}^{p}$ and
$\mathbb{R}^{q}$ (respectively), the forward operator is a Euclidean map
$\boldsymbol{g}$, and the noise is typically written into the IP explicitly.
The solution is

$\displaystyle\boldsymbol{\theta}$ $\displaystyle\text{\quad such that
\quad}\boldsymbol{g}(\boldsymbol{\theta})+\boldsymbol{\epsilon}=\boldsymbol{y}\quad\quad\big{(}\boldsymbol{\theta}\in\mathcal{P},\
\boldsymbol{y}\in\mathcal{Q},\ \boldsymbol{\epsilon}\sim
F_{\boldsymbol{E}}\big{)}\ .$ (2)

The term $\boldsymbol{\epsilon}$ is a realization of a random variable
$\boldsymbol{E}$ corresponding to the error process with cumulative
distribution function $F_{\boldsymbol{E}}$. (The use of the cumulative
function allows for continuous as well as discrete and mixed random variables.
Also, the statement within (2) assumes additive error, but is trivially
modified for multiplicative or other error structures.) Together,
$\boldsymbol{g}(\boldsymbol{\theta})$ and
$F_{\boldsymbol{E}}(\boldsymbol{\epsilon})$ determine a model whose likelihood
function forms the basis for _statistical_ IPs (Pawitan, 2001; Davison, 2003;
Kaipio and Somersalo, 2005; Tarantola, 2005; Stuart, 2010; Reid, 2010;
Chiachio-Ruano et al., 2022). As stated in (2), the solution to the
statistical IP is a point estimate derived in some fashion from the
likelihood. In certain cases, the maximum likelihood estimate will be a
weighted or regularized least-squares estimate. In the Bayesian (Bernardo and
Smith, 1994; Robert, 2007; Gelman et al., 2014) and Fiducial (Hannig, 2009;
Hannig et al., 2016) inferential contexts, the solution to a statistical IP is
an entire probability distribution from which point estimates may be obtained.

Over the past century, the field of inverse problems has extended into nearly
every domain of science, engineering, and technology, cementing its status as
fundamental to applied mathematics and statistics. One particular, more
recently developed class of so-called _stochastic_ inverse problems (SIPs)
identifies a solution as a probability density whose pushforward takes a given
form (Butler et al., 2014, 2015a, 2015b; Mattis et al., 2015; Butler et al.,
2018a; Mattis and Butler, 2019; Uy and Grigoriu, 2019; Butler et al., 2020a;
Bruder et al., 2020; Butler et al., 2020b; Tran and Wildey, 2021). These SIPs
and their solutions have also gone under the names of “(stochastic) inverse
sensitivity problems” (Breidt et al., 2011; Butler et al., 2012, 2014; Butler
and Estep, 2013; Graham et al., 2017); “measure-theoretic inverse problems”
(Butler et al., 2017; Presho et al., 2017); “consistent Bayesian” or
“pushforward-based inference” (Butler et al., 2018a, b; Walsh et al., 2018);
“data/observation-consistent inversion” (Butler et al., 2018b, 2020c; Butler
and Hakula, 2020; Mattis et al., 2022); or “random parameter models” (Swigon
et al., 2019).

The formulation and solution of these SIPs can look peculiar to those familiar
with more traditional inverse problems. For example, each observable quantity
is an entire probability density, not of a fixed realization of a random
variable modelled conditionally. Also, there must be at least as many unknown
parameters as observables. Moreover, derivations of the two different classes
of solutions that have been proposed—those of Breidt et al. (2011) and Butler
et al. (2018a) which we call “BBE” (short for “Breidt, Butler, Estep”) and
“BJW” (short for “Butler, Jakeman, Wildey”), respectively—rely heavily upon
measure theory, specifically the disintegration theorem.

The goal of this paper is to explore the formulation and solution to SIPs
using only introductory probability and mathematical statistics while also
giving a careful examination of the existing literature. While measure theory
is often the language of choice, we feel that in this context it can obscure
concepts that are quite uncomplicated. To begin to explore some of the SIP
peculiarities, we first address the case when $p=q$. For $p\geq q$, we then
offer a class of “intuitive” solutions and show how these relate to a change-
of-variables (CoV) from observable- to parameter-space (Section 3). We give a
simple algorithm to obtain samples from intuitive solutions, and the
underlying reasoning proves useful for later results. The BBE and BJW
solutions are investigated through theoretical derivation and analytic
examples (Sections 4, 5). Analytic results have heretofore never been given.
We show that under mild assumptions, the BBE and BJW solutions can be derived
from CoVs. Some discussion of related work is provided in Section 6. We then
show that _any_ solution to an SIP must be related to a CoV (Section 7). For
$p>q$, the results rely on _auxiliary variables_ to augment the
underdetermined system. A comparative critique of SIPs/CoVs versus inference
is given in Section 8 along with two illustrative examples. Our main
conclusion is that SIPs are significantly different from statistical inverse
problems and should be handled with greater care.

In the remainder of this section, we give a formal definition of the SIP and a
statement of the traditional CoV theorem.

###  1.1 Stochastic Inverse Problem Formulation and Assumptions

In an SIP, a
$q$-dimensional vector of observable quantities or data is taken to be a
random variable $\boldsymbol{Y}$ with given probability density function (pdf)
$f_{\boldsymbol{Y}}(\boldsymbol{y})$. Further, there is also a function or
“forward-map” from parameter- to data-space

$\displaystyle\boldsymbol{g}:\big{(}\mathcal{P}\subseteq\mathbb{R}^{p}\big{)}\
$ $\displaystyle\longmapsto\ \big{(}\mathcal{Q}\subseteq\mathbb{R}^{q}\big{)}$
$\displaystyle\left[\begin{array}[]{l}\Theta_{1}\\\ \vdots\\\ \vdots\\\
\Theta_{p}\end{array}\right]$
$\displaystyle\longmapsto\left[\begin{array}[]{lcl}Y_{1}&=&g_{1}(\boldsymbol{\Theta})\\\
\vdots&&\\\ Y_{q}&=&g_{q}(\boldsymbol{\Theta})\end{array}\right],$

with $p\geq q$, which is either known analytically or which can be evaluated
as a “black-box”, such as a computer model that solves a set of differential
equations. The goal is to obtain a density for random variables in the pre-
image, $\boldsymbol{\Theta}\in\mathcal{P}$, that will transform (exactly or in
some approximate sense) to $f_{\boldsymbol{Y}}(\boldsymbol{y})$ under the
forward-map $\boldsymbol{g}$. In other words, given $f_{\boldsymbol{Y}}$ and
$\boldsymbol{g}$, the solution is

$\displaystyle f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle\text{\quad such that
\quad}\boldsymbol{\Gamma}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{g}(\boldsymbol{\Theta})\sim
f_{\boldsymbol{Y}}$ (3)

i.e., a density that pushes forward or propagates “correctly”. _Throughout
this paper we assume that a solution exists_ , though in applications one may
need to be careful and check that the range $\boldsymbol{g}(\mathcal{P})$
contains $\mathcal{Q}$, the support of the given
$f_{\boldsymbol{Y}}(\boldsymbol{y})$.

Within the classical IP definition (1), taking $\mathcal{P}$ and $\mathcal{Q}$
to be subsets of Euclidean space, and replacing $\boldsymbol{y}$ and
$\boldsymbol{\theta}$ with continuous random variables $\boldsymbol{Y}$ and
$\boldsymbol{\Theta}$, then an SIP defined by (3) appears to be a natural
variant of the more traditional (1). The simplest example is when the operator
$\mathcal{F}$ is a linear map $\boldsymbol{g}$ between Euclidean spaces and
represented by the invertible matrix $\boldsymbol{A}$. The solution to the
classical linear IP is of course
$\boldsymbol{\theta}=\boldsymbol{A}^{-1}\boldsymbol{y}$; the solution to the
linear SIP is the density of the random variable
$\boldsymbol{A}^{-1}\boldsymbol{Y}$.

In practice, a random sample from
$f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ constitutes a solution as well,
albeit an _approximate_ solution whose quality increases with sample size. Let
us explicitly state this as a non-controversial assumption.

_ Assumption (A0): The SIP defined by (3) is approximately solved when a
random sample $\boldsymbol{\theta}^{(1)},\ldots,\boldsymbol{\theta}^{(M)}$ is
obtained such that
$\boldsymbol{g}\big{(}\boldsymbol{\theta}^{(1)}\big{)},\ldots\boldsymbol{g}\big{(}\boldsymbol{\theta}^{(M)}\big{)}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}f_{\boldsymbol{Y}}$.
_

The main assumption that we use in this paper is the following.

_ Assumption (A1): The maps $g_{1},\ldots,g_{q}$ are in $C^{1}(\mathcal{P})$
(i.e. they are continuously differentiable on the domain) and the Jacobian
$\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\big{|}$ has
full row-rank except for on a set $\mathcal{P}_{0}$ of measure zero. Without
loss of generality, the left $q\times q$ block of the Jacobian
$\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}_{1:q}}\big{|}$
is invertible on all of $\mathcal{P}\setminus\mathcal{P}_{0}$. _

First of all, it is not unreasonable to assume the Jacobian has linearly
independent rows almost everywhere (“a.e.”, meaning for all but a set of
measure zero) in $\mathcal{P}$. If the rows were dependent a.e., then one
could consider only the independent outputs and instead solve the
corresponding reduced SIP. Second, (A1) allows the domain to be partitioned
into subdomains $\mathcal{P}_{0}$ where the Jacobian determinant vanishes and
disjoint sets $\mathcal{P}_{1},\ldots,\mathcal{P}_{m}$ where the functions
$g_{1},\ldots,g_{q}$ produce full row-rank Jacobian.

In the first paragraph of this section it was taken as given that:

_ Assumption (A2): The random variables $\boldsymbol{Y}$ and
$\boldsymbol{\Theta}$ are absolutely continuous and thus admit probability
densities with respect to Lebesgue measure. _

This assumption is not strictly necessary but is, nonetheless, the usual
context for the SIPs. The results of this paper carry over to discrete random
variables as well, but we shall from here on only speak in terms of
probability densities.

###  1.2 Changes-of-Variables

For completeness we state the change-of-variables (CoV) theorem, which can be
found in standard textbooks on probability and mathematical statistics. The
following is taken from Shao (2003), pg. 23 and is also similar to Casella and
Berger (2002), pg. 185.

######  Theorem (Change-of-Variables).

Let $\boldsymbol{U}$ be a random $p$-vector with a Lebesgue pdf
$f_{\boldsymbol{U}}(\boldsymbol{u})$ and let
$\boldsymbol{V}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{t}(\boldsymbol{U})$,
where $\boldsymbol{t}$ is a Borel function from
$(\mathbb{R}^{p},\mathcal{B}^{p})$ to $(\mathbb{R}^{p},\mathcal{B}^{p})$. Let
$A_{1},\ldots,A_{d}$ be disjoint sets in $\mathcal{B}^{p}$ such that
$\mathbb{R}^{p}-(A_{1}\cup\cdots\cup A_{d})$ has Lebesgue measure 0 and
$\boldsymbol{t}$ on $A_{i}$ is one-to-one with a nonvanishing Jacobian, i.e.,
$|\frac{\partial\boldsymbol{t}}{\partial\boldsymbol{u}}|\neq 0$ on $A_{i}$ for
$i=1,\ldots,d$. Then $\boldsymbol{V}$ has the following Lebesgue pdf

$\displaystyle f_{\boldsymbol{V}}(\boldsymbol{v})$
$\displaystyle=\sum_{i=1}^{d}f_{\boldsymbol{U}}\big{(}\boldsymbol{u}=\boldsymbol{t}^{-1}_{i}(\boldsymbol{v})\big{)}\Big{|}\frac{\partial\boldsymbol{t}^{-1}_{i}}{\partial\boldsymbol{v}}\Big{|}\
,$ (4)

where $\boldsymbol{t}^{-1}_{i}$ is the inverse function of $\boldsymbol{t}$ on
$A_{i}$ (i=1, …, d).

Being that one of the goals of the paper is to keep exposition at the level of
basic probability and away from measure theory, we note that the theorem above
will be the only place in which Borel sigma algebras are mentioned, and the
last time we need to refer to Lebesgue measure.

##  2 CoV Solutions to the SIP When $p=q$: Immediate Issues and Insights

Taking $\boldsymbol{t}=\boldsymbol{g}$ and
$\boldsymbol{U}=\boldsymbol{\Theta}$ in the CoV Theorem, the pushforward of
$f_{\boldsymbol{\Theta}}$ through $\boldsymbol{g}$ has density
$f_{\boldsymbol{\Gamma}}(\boldsymbol{\gamma})=\sum_{i=1}^{d}f_{\boldsymbol{\Theta}}\big{(}\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})\big{)}\big{|}\partial\boldsymbol{g}_{i}^{-1}/\partial\boldsymbol{\gamma}\big{|}$.
If $\boldsymbol{g}$ is a one-to-one function of $\mathcal{P}$ onto
$\mathcal{Q}$, then $d=1$, and no $i$ subscript is necessary. In this case one
can go the other direction and pull $f_{\boldsymbol{Y}}$ back through
$\boldsymbol{h}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{g}^{-1}$ to
get
$f_{\boldsymbol{H}}(\boldsymbol{\eta})=f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\eta})\big{)}\big{|}\partial\boldsymbol{g}/\partial\boldsymbol{\eta}\big{|}$.
Thus, by taking
$\boldsymbol{\Theta}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{H}$,
the SIP is solved uniquely (a.e.) by the density

$\displaystyle f^{\text{CoV}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}\
.$ (5)

If $\boldsymbol{g}$ is not a uniquely invertible function, its many-to-one
nature allows for multiple solutions to the SIP. We show in the next
proposition that there are in fact _infinitely_ many solutions. This property
does not appear to have been observed within the existing literature. For
simplicity and ease of exposition we assume that $\boldsymbol{g}$ is exactly
$m$-to-1 onto all of $\mathcal{Q}$. The general case can be proved similarly,
but involves keeping track of all the distinct sets where the function is two-
to-one, three-to-one, etc. within the $d$-partition of the domain; as such the
indexing is more laborious and no further insight is added.

######  Proposition 1.

Assume (A1) and let $\boldsymbol{g}$ be an $m$-to-1 ($m>1$) function of
$\mathcal{P}$ onto $\mathcal{Q}$. Then there exists a continuously-indexed,
infinite family of solutions to the SIP.

###### Proof.

Let $\mathbbm{1}\\{\boldsymbol{\theta}\in S\\}$ denote an indicator function
that is 1 when $\boldsymbol{\theta}$ is in the set $S$, and 0 otherwise.

Consider the density

$\displaystyle
f^{\text{CoV}}_{\boldsymbol{\Theta},\boldsymbol{w}}(\boldsymbol{\theta})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}\sum_{j=1}^{m}w_{j}\
\mathbbm{1}\\{\boldsymbol{\theta}\in\mathcal{P}_{j}\\}$ (6)

continuously indexed by mixture weights summing to unity:
$\sum_{j=1}^{m}w_{j}=1$. By the CoV Theorem, the pushforward density for
$\boldsymbol{\Gamma}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{g}(\boldsymbol{\Theta})$
is

$\displaystyle f_{\boldsymbol{\Gamma}}(\boldsymbol{\gamma})$
$\displaystyle=\sum_{i=1}^{m}f^{\text{CoV}}_{\boldsymbol{\Theta},\boldsymbol{w}}\big{(}\boldsymbol{\theta}=\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})\big{)}\
\Big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{\gamma}}\Big{|}$
$\displaystyle=\sum_{i=1}^{m}\left(f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})\big{)}\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}_{\boldsymbol{\theta}=\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})}\sum_{j=1}^{m}w_{j}\
\mathbbm{1}\\{\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})\in\mathcal{P}_{j}\\}\right)\Big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{\gamma}}\Big{|}$
$\displaystyle=f_{\boldsymbol{Y}}(\boldsymbol{\gamma})\sum_{i=1}^{m}\left(\sum_{j=1}^{m}w_{j}\
\mathbbm{1}\\{\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})\in\mathcal{P}_{j}\\}\right)$
$\displaystyle=f_{\boldsymbol{Y}}(\boldsymbol{\gamma})\sum_{i=1}^{m}w_{i}\
\mathbbm{1}\\{\boldsymbol{g}_{i}^{-1}(\boldsymbol{\gamma})\in\mathcal{P}_{i}\\}\quad=f_{\boldsymbol{Y}}(\boldsymbol{\gamma})$

as desired. In the penultimate line, all of the $i\neq j$ terms vanish
tautologically. ∎

Section A of the appendix gives a concrete demonstration of the method used in
the proof above. In the general case, any sets in the domain that get mapped
to the same range subset can be given their own set of weights summing to
unity. In summary, when $p=q$, the existence of a neighborhood within
$\mathcal{P}$ where $\boldsymbol{g}$ is not one-to-one implies an infinite
number of solutions to the SIP. The following related result states that there
are no other solutions to the SIP beyond those found by CoVs.

######  Proposition 2.

Suppose that $p=q$ and assumption (A1) holds for the forward map
$\boldsymbol{g}$. If $f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ solves the
SIP (3), then it is (a.e.) derivable from changes-of-variables.

###### Proof.

By (A1), the domain $\mathcal{P}\setminus\mathcal{P}_{0}$ can be partitioned
into $\mathcal{P}_{1},\ldots,\mathcal{P}_{d}$ such that $\boldsymbol{g}$ has
$C^{1}(\mathcal{Q})$ inverses
$\boldsymbol{g}_{1}^{-1},\ldots,\boldsymbol{g}_{d}^{-1}$ defined on these
sets. The density $f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ solves the
SIP, so then by the CoV Theorem, its pushforward under $\boldsymbol{g}$ is
$f_{\boldsymbol{Y}}(\boldsymbol{y})=\sum_{i=1}^{d}f_{\boldsymbol{\Theta}}\big{(}\boldsymbol{g}_{i}^{-1}(\boldsymbol{y})\big{)}\
\big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{y}}\big{|}$.
Let the range sets be denoted
$\mathcal{Q}_{i}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{g}(\mathcal{P}_{i})$
(while not a proper partition, it still holds that
$\mathcal{Q}=\bigcup_{i=1}^{d}\mathcal{Q}_{i}$). Taking

$\displaystyle f_{\boldsymbol{Y}_{i}}(\boldsymbol{y})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}w_{i}^{-1}f_{\boldsymbol{\Theta}}\big{(}\boldsymbol{g}_{i}^{-1}(\boldsymbol{y})\big{)}\
\Big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{y}}\Big{|}\
\mathbbm{1}\\{\boldsymbol{y}\in\mathcal{Q}_{i}\\}$ $\displaystyle w_{i}$
$\displaystyle=\int_{\mathcal{Q}_{i}}f_{\boldsymbol{\Theta}}\big{(}\boldsymbol{g}_{i}^{-1}(\boldsymbol{y})\big{)}\
\Big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{y}}\Big{|}d\boldsymbol{y}\quad=\int_{\mathcal{P}_{i}}f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})d\boldsymbol{\theta}\
,$

a CoV from each $\mathcal{Q}_{i}$ to $\mathcal{P}_{i}$ under
$\boldsymbol{g}_{i}^{-1}$ yields densities
$f^{\text{CoV}}_{\boldsymbol{\Theta}_{i}}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}w_{i}^{-1}f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})\mathbbm{1}\\{\boldsymbol{\theta}\in\mathcal{P}_{i}\\}$.
The weighted combination $\sum_{i=1}^{d}w_{i}\
f^{\text{CoV}}_{\boldsymbol{\Theta}_{i}}(\boldsymbol{\theta})$ of the $d$
range subset CoVs is then the original given solution to the SIP, up to a set
of exceptional points having measure zero. ∎

Regardless of whether $\boldsymbol{g}$ is invertible everywhere, if the
pushforward of $f_{\boldsymbol{\Theta}}$ is $f_{\boldsymbol{Y}}$, then the
original density can be written concisely as a CoV from $\mathcal{Q}$ to
$\mathcal{P}\setminus\mathcal{P}_{0}$. To see this note that by the Inverse
Function Theorem,
$\big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{y}}\big{|}_{\boldsymbol{y}=\boldsymbol{g}(\boldsymbol{\theta})}=\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\big{|}^{-1}$
for all $i$ so that

$\displaystyle f_{\boldsymbol{Y}}(\boldsymbol{y})$
$\displaystyle=\sum_{i=1}^{d}f_{\boldsymbol{\Theta}}\big{(}\boldsymbol{g}_{i}^{-1}(\boldsymbol{y})\big{)}\
\Big{|}\frac{\partial\boldsymbol{g}_{i}^{-1}}{\partial\boldsymbol{y}}\Big{|}$
$\displaystyle\Rightarrow\quad
f_{\boldsymbol{Y}}(\boldsymbol{g}(\boldsymbol{\theta}))$
$\displaystyle=\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}^{-1}\
\sum_{i=1}^{d}f_{\boldsymbol{\Theta}}\big{(}\boldsymbol{g}_{i}^{-1}(\boldsymbol{g}(\boldsymbol{\theta}))\big{)}$
$\displaystyle=\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}^{-1}\
\sum_{i=1}^{d}f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})\
\mathbbm{1}\\{\boldsymbol{\theta}\in\mathcal{P}_{i}\\}$
$\displaystyle\Rightarrow\quad f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=f_{\boldsymbol{Y}}(\boldsymbol{g}(\boldsymbol{\theta}))\
\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}\ .$
(7)

This is a general statement that a density can be written in terms of its
pushforward, and will be used in Proposition 4.

##  3 Intuitive Solutions to SIPs

It was demonstrated in the last section that when $p=q$, the SIP (3) is solved
by using at least one density given by the CoV theorem. However, for the more
general case of $p>q$, there is also a family of densities that solve the SIP
(3) which we call “intuitive solutions”.

Let the parameter space be partitioned as
$\boldsymbol{\Theta}^{\top}=\big{(}\boldsymbol{\Theta}_{1:q}^{\top},\boldsymbol{\Theta}_{(q+1):p}^{\top}\big{)}$.
Temporarily fix
$\boldsymbol{\Theta}_{(q+1):p}=\boldsymbol{\theta}^{*}_{(q+1):p}$ and consider
the distribution that comes from solving the “square” SIP:

$\displaystyle
f_{\boldsymbol{\Theta}_{1:q}|(\boldsymbol{\Theta}_{(q+1):p}=\boldsymbol{\theta}^{*}_{(q+1):p})}^{\text{CoV}}$
$\displaystyle\text{\quad such that
\quad}\boldsymbol{g}\big{(}\boldsymbol{\Theta}_{1:q},\boldsymbol{\theta}^{*}_{(q+1):p}\big{)}\sim
f_{\boldsymbol{Y}}$ (8)

via a $q$-dimensional change-of-variables. This implies that the intuitive
solution $f_{\boldsymbol{\Theta}}^{\text{Int}}$ given by

$\displaystyle f_{\boldsymbol{\Theta}}^{\text{Int}}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}f_{\boldsymbol{\Theta}_{1:q}|\boldsymbol{\Theta}_{(q+1):p}}^{\text{CoV}}\cdot
f_{\boldsymbol{\Theta}_{(q+1):p}}$ (9)

solves the SIP for any choice of $f_{\boldsymbol{\Theta}_{(q+1):p}}$.

Sampling from this density constitutes a solution by (A0), and this can be
performed using a simple and intuitive (hence the name) Monte Carlo routine
within the algorithm below.

Algorithm 1 Intuitive solutions to the stochastic inverse problem (SIP)

Given densities $f_{\boldsymbol{Y}}$ and $f_{\boldsymbol{\Theta}_{(q+1):p}}$,

0: Sample $\boldsymbol{y}^{*}\sim f_{\boldsymbol{Y}}$ and
$\boldsymbol{\theta}_{(q+1):p}^{*}\sim f_{\boldsymbol{\Theta}_{(q+1):p}}$.

0: Solve the (likely nonlinear) equation
$\displaystyle\boldsymbol{y}^{*}-\boldsymbol{g}\big{(}\boldsymbol{\theta}_{1:q}^{*},\boldsymbol{\theta}_{(q+1):p}^{*}\big{)}=\boldsymbol{0}$
(10) for $\boldsymbol{\theta}_{1:q}^{*}$ to form the random sample
$\big{(}\boldsymbol{\theta}_{1:q}^{*\top},\boldsymbol{\theta}_{(q+1):p}^{*\top}\big{)}^{\top}$,
a realization of $\boldsymbol{\Theta}\sim
f^{\text{Int}}_{\boldsymbol{\Theta}}$.

0: Repeat $\langle 1\rangle-\langle 2\rangle$.

Algorithm 1 requires no sophisticated sampling techniques, but it does require
a nonlinear solver. Typical solvers benefit from the ability to evaluate and
use the Jacobian, but this is not always required. This algorithm was
mentioned in Swigon et al. (2019) for the simplest case of $p=q$, though the
authors ultimately avoided the use of solvers and instead employed a
Metropolis-Hastings algorithm based upon an approximated Jacobian term.
Intuitive solutions as given above were not considered as an alternative
within the work stemming from the two major approaches of BBE (Breidt et al.,
2011) and BJW (Butler et al., 2018a).

When the Jacobian of the function $\boldsymbol{g}$ vanishes within the domain
$\mathcal{P}$, the input space must be partitioned into $d>1$ disjoint sets
according to the CoV Theorem. This will lead to the same infinite class of
solutions as in Proposition 1 of the previous section. When using
deterministic equation solvers, the distribution of starting values for
$\boldsymbol{\theta}^{*}_{1:q}$ will determine the particular density out of
the infinite family.

The main characteristic of intuitive solutions is made explicitly clear in
Algorithm 1. In the first step, $\boldsymbol{y}^{*}\sim f_{\boldsymbol{Y}}$
and $\boldsymbol{\theta}_{(q+1):p}^{*}\sim f_{\boldsymbol{\Theta}_{(q+1):p}}$
are drawn _independently_ , not from some joint density. This implies that the
solution has the feature that the response and $(p-q)$ of the parameters do
not covary—the output is independent of these inputs! Thus, when forming an
intuitive solution to the SIP, one might specify a density for $(p-q)$ of the
least important input parameters, as determined by an initial sensitivity
analysis.

As a final note about intuitive solutions, the relationship between the
density (9) and the sampling steps of Algorithm 1 will prove useful in later
sections (specifically, 4.1, 5.4, and 7).

##  4  Breidt et al. 2011 (BBE) and Related Work

The IP under current consideration was first (to the best of our knowledge)
described in Breidt et al. (2011) under the name “inverse sensitivity problem”
in the first of three papers (Parts I, II, III). The authors give an algorithm
to compute an approximate pdf of the solution, which we denote as
$\widehat{f}^{\text{BBE}}_{\boldsymbol{\Theta}}$ for a single output ($q=1$).
However, the authors do not specify the exact solution
$f^{\text{BBE}}_{\boldsymbol{\Theta}}$ that should be approximated by
$\widehat{f}^{\text{BBE}}_{\boldsymbol{\Theta}}$. This approximate solution
relies upon (A1), (A2), and the presence of derivative information, obtained
either analytically or estimated via adjoint techniques. The algorithm and
hence the approximate solution
$\widehat{f}^{\text{BBE}}_{\boldsymbol{\Theta}}$ depends heavily upon
discretization and as such is restricted to very low dimension $p$.

Part II (Butler et al., 2012) gives a rigorous error analysis of Breidt et al.
(2011), treating both statistical error due to sampling of the data density
$f_{\boldsymbol{Y}}$, and numerical error due to solving differential
equations during the likelihood calculation. Part III (Butler et al., 2014)
considers the multiresponse scenario, i.e., $q>1$. The Part III approach
deviates from Parts I and II in that the authors discretize events in the
parameter space as opposed to manifolds in the observable space. This
framework leads to the first mention of the disintegration theorem of measure
theory.

The last methodological developments within the BBE framework are Butler et
al. (2017) and Mattis and Butler (2019). Butler et al. (2017) discusses the
necessary event approximations and considers an adaptive sampling algorithm to
solve the SIP. Mattis and Butler (2019) focuses on computational estimates of
events from samples, using adjoints to enhance low-order piecewise-defined
surrogate response surfaces.

Several subsequent works have explored the use of
$\widehat{f}^{\text{BBE}}_{\boldsymbol{\Theta}}$ in various contexts. Butler
and Estep (2013) illustrates the BBE approach using the Brusselator
reaction–diffusion model. Butler et al. (2015b) compares the multiresponse
methodology of Butler et al. (2014) to other UQ methods in the context of
material damage from vibrational data. Mattis et al. (2015) applies the
multiresponse BBE to a groundwater contamination problem. Eight parameters
(including source location coordinates, the contaminant source flux, etc.) are
inferred from seven data points (the concentrations of the contaminant in
seven wells). Similarly, Butler et al. (2015a) describes the SIP for an
application in hydrodynamic models. The authors infer two parameters—the two-
dimensional vector of the Manning’s $n$ parameter—given maximum water
elevations at two of twelve possible observation stations. Note that, although
data is available from all twelve stations, only one or two may be used in any
given SIP, due to the constraint that $p\geq q$. To address this peculiarity
of the method, the authors introduce the notion of “geometrically distinct”
data (or quantities of interest, “QoI”): since $q\leq p$, then each data point
should contain as much information as possible (analogous to linear
independence). At the same time (and as pointed out by the authors), two
different sets of geometrically distinct data may lead to very different
solutions for the otherwise same SIP. Graham et al. (2017) also estimates
probable Manning’s $n$ fields using a storm surge model and the BBE SIP
framework. Finally, Presho et al. (2017) uses the BBE framework together with
the generalized multiscale finite element method for uncertainty
quantification within two-phase flow problems.

More recently, Uy and Grigoriu (2019) provides useful examples and commentary
on SIPs, with an emphasis on Breidt et al. (2011); we comment further on this
paper in Section 6.

###  4.1 The BBE Solution to the SIP is a Change-of-Variables

In their solutions to the univariate and multivariate inverse sensitivity
problem, Breidt et al. (2011) and Butler et al. (2014) provide algorithms to
discretize a probability density that is never actually given in closed form.
We begin our reanalysis by deriving what this density must be. The derivation
will clarify and extend arguments of Uy and Grigoriu (2019) (Sec. 2) without
the explicit use of the disintegration theorem.

Some notation and concepts used within the BBE approach are needed. In the BBE
setup there is a transverse parameterization
$\boldsymbol{t}(\boldsymbol{\theta})$ of dimension $q$, and a
$(p-q)$-dimensional contour parameterization
$\boldsymbol{c}(\boldsymbol{\theta})$ of the forward map contours
$\boldsymbol{g}$ in $\mathcal{Q}$. Furthermore, there exists a one-to-one
function $\boldsymbol{g}^{t}$ between $\boldsymbol{t}(\boldsymbol{\theta})$
and $\mathcal{Q}$ induced by $\boldsymbol{g}$, although in practice this
function is typically unknown. The function $\boldsymbol{g}^{t}$ simply
transforms the transverse coordinates uniquely to the observable space
$\mathcal{Q}$ and thus operates as a $q$-dimensional forward map.

It is easiest to derive the density of the BBE solution by first thinking of
how to sample from it, and this can be done akin to Algorithm 1 for the
intuitive solutions. Given a random sample $\boldsymbol{y}^{*}\sim
f_{\boldsymbol{Y}}$, the solution to
$\boldsymbol{y}^{*}-\boldsymbol{g}^{t}(\boldsymbol{t}^{*})=\boldsymbol{0}$ is
a random draw from the density
$f_{\boldsymbol{T}}(\boldsymbol{t})=f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}^{t}(\boldsymbol{t})\big{)}\big{|}\frac{\partial\boldsymbol{g}^{t}}{\partial\boldsymbol{t}}\big{|}$.
Next, a uniformly sampled point along the contour indexed by
$\boldsymbol{t}^{*}$ is a realization $\boldsymbol{c}^{*}\sim
f_{\boldsymbol{C}|\boldsymbol{T}}(\boldsymbol{c}\ |\
\boldsymbol{T}=\boldsymbol{t}^{*})$ for the special case of a uniform
distribution. Thus the pair $(\boldsymbol{t}^{*},\boldsymbol{c}^{*})$ is a
random draw from the joint distribution
$f_{\boldsymbol{T},\boldsymbol{C}}=f_{\boldsymbol{T}}\cdot
f_{\boldsymbol{C}|\boldsymbol{T}}$. Putting this back in
$\boldsymbol{\theta}$-coordinates requires one more change-of-variables and
assumes that the map
$(\boldsymbol{t},\boldsymbol{c})\mapsto\boldsymbol{\theta}$ has nonvanishing
Jacobian a.e. The resulting density is thus given by the following iterated
change-of-variables:

$\displaystyle f^{\text{BBE}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=f_{\boldsymbol{T},\boldsymbol{C}}\Big{(}\boldsymbol{c}(\boldsymbol{\theta}),\boldsymbol{t}(\boldsymbol{\theta})\Big{)}\Big{|}\frac{\partial(\boldsymbol{t},\boldsymbol{c})}{\partial\boldsymbol{\theta}}\Big{|}$
$\displaystyle=f_{\boldsymbol{T}}(\boldsymbol{t}(\boldsymbol{\theta}))f_{\boldsymbol{C}|\boldsymbol{T}}\Big{(}\boldsymbol{c}(\boldsymbol{\theta})\
|\
\boldsymbol{t}(\boldsymbol{\theta})\Big{)}\Big{|}\frac{\partial(\boldsymbol{t},\boldsymbol{c})}{\partial\boldsymbol{\theta}}\Big{|}$
$\displaystyle=\left(f_{\boldsymbol{Y}}\Big{(}\boldsymbol{g}^{t}\big{(}\boldsymbol{t}(\boldsymbol{\theta})\big{)}\Big{)}\Big{|}\frac{\partial\boldsymbol{g}^{t}}{\partial\boldsymbol{t}}\Big{|}_{\boldsymbol{t}(\boldsymbol{\theta})}\
f_{\boldsymbol{C}|\boldsymbol{T}}\Big{(}\boldsymbol{c}(\boldsymbol{\theta})\
|\
\boldsymbol{t}(\boldsymbol{\theta})\Big{)}\right)\left|\frac{\partial(\boldsymbol{t},\boldsymbol{c})}{\partial\boldsymbol{\theta}}\right|$
(11)

for $f_{\boldsymbol{C}|\boldsymbol{T}}$ uniform. We have just sketched a proof
of the following result.

######  Proposition 3.

Suppose that $\boldsymbol{g}^{t}(\boldsymbol{t})$ (induced by the full forward
map $\boldsymbol{g}$) is a continuously differentiable, invertible map from
the $q$-dimensional image of $\mathcal{P}$ under
$\boldsymbol{t}(\boldsymbol{\theta})$ to $\mathcal{Q}$. Also suppose that
together the transverse and contour parameterizations
$\boldsymbol{t}(\boldsymbol{\theta})$ and
$\boldsymbol{c}(\boldsymbol{\theta})$ jointly map to the parameter space
$\mathcal{P}$ with nonvanishing Jacobian according to (A1). Then the exact BBE
solution is equivalent to an iterated change-of-variables given by (11).

The main difficulty with directly using the CoV form of the BBE solution above
is the lack of closed-form transverse and contour parameterizations. Hence,
examples in which the exact density $f^{\text{BBE}}(\boldsymbol{\theta})$ is
available are hard to come by. The following two examples do however admit
closed-form solutions and illustrate the concepts above.

###  4.2 Example: Linear Map

Here we find the exact solution density corresponding to Section 3.1 of Butler
et al. (2014). Let
$\boldsymbol{g}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}\boldsymbol{\theta}$
so that $\boldsymbol{Y}=\boldsymbol{A}\boldsymbol{\Theta}$. Suppose that
$\boldsymbol{A}^{\perp}$ is a $(p-q)\times p$ representation of the null space
of row($\boldsymbol{A}$).

Consider the transverse parameterization that follows the Jacobian of the
forward map, namely,
$\boldsymbol{t}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}\boldsymbol{\theta}$.
The connection between the data-space and transverse coordinates is the
trivial one:
$\boldsymbol{y}=\boldsymbol{g}^{t}(\boldsymbol{t})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{I}_{q}\boldsymbol{t}$.
Thus we have
$f_{\boldsymbol{T}}(\boldsymbol{t})=f_{\boldsymbol{Y}}(\boldsymbol{t})$.
Contours of $\boldsymbol{g}$ are described by
$\boldsymbol{c}(\boldsymbol{\theta};\boldsymbol{t})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}^{\perp}\boldsymbol{\theta}$.
The choice of uniformly distributed contours implies that
$f_{\boldsymbol{C}|\boldsymbol{T}}(\boldsymbol{c}|\boldsymbol{t})\propto\mathbbm{1}\big{\\{}\boldsymbol{l}(\boldsymbol{t})\leq\boldsymbol{c}\leq\boldsymbol{u}(\boldsymbol{t})\big{\\}}$,
where $\boldsymbol{l}$ and $\boldsymbol{u}$ are given lower and upper bounds
(in order to define a valid probability distribution). If the pre-specified
bounds are aligned with $\boldsymbol{A},\boldsymbol{A}^{\perp}$, then
$\boldsymbol{l}$ and $\boldsymbol{u}$ will not depend on $\boldsymbol{t}$.

Transforming $(\boldsymbol{T},\boldsymbol{C})$ to $\boldsymbol{\Theta}$ comes
from the invertible augmented system

$\displaystyle\begin{bmatrix}\boldsymbol{t}\\\ \boldsymbol{c}\end{bmatrix}$
$\displaystyle=\left(\boldsymbol{A}_{\text{aug}}\stackrel{{\scriptstyle\text{def}}}{{=}}\begin{bmatrix}\boldsymbol{A}\\\
\boldsymbol{A}^{\perp}\end{bmatrix}\right)\boldsymbol{\theta}$

having constant Jacobian $\left|\boldsymbol{A}_{\text{aug}}\right|$, and this
implies the result

$\displaystyle f^{\text{BBE}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle\propto f_{\boldsymbol{Y}}(\boldsymbol{A}\boldsymbol{\theta})\
\mathbbm{1}\big{\\{}\boldsymbol{l}(\boldsymbol{A}\boldsymbol{\theta})\leq\boldsymbol{A}^{\perp}\boldsymbol{\theta}\leq\boldsymbol{u}(\boldsymbol{A}\boldsymbol{\theta})\big{\\}}\
.$

The top row of Figure 1 shows samples of the BBE solution overlaid upon
contours of the forward map
$g(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}[-1/3,4/3]\boldsymbol{\theta}$
(left). The observable density $f_{Y}(y)$ is a $N\big{(}1/2,(1/4)^{2}\big{)}$
truncated to (0,1); this is plotted over the pushforward histogram of the BBE
solution samples (right), confirming the solution.

Figure 1:  Left: Samples of $\boldsymbol{\Theta}\sim
f^{\text{BBE}}_{\boldsymbol{\Theta}}$ and heatmap of
$f^{\text{BBE}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ with contours of
the function $g(\boldsymbol{\theta})$ overlaid; Right: The pushforward of the
BBE solution samples ($g(\boldsymbol{\Theta})$, histogram) compared to the
given density $f_{Y}(y)$. Top row: Linear map in Section 4.2; Bottom row:
Nonlinear map in Section 4.3.

###  4.3 Example: Nonlinear Map

Now we work through Example 2 of Butler et al. (2014). Let
$Y=g(\boldsymbol{\Theta})$ with
$g(\theta_{1},\theta_{2})\stackrel{{\scriptstyle\text{def}}}{{=}}\tfrac{1}{2}\big{(}\theta_{1}^{2}+\theta_{2}^{2}\big{)}$
with $0<\theta_{1},\theta_{2}\leq 1$. (Note we have slightly modified the
forward map by introducing the factor of $1/2$, simply to clean up the
results.) Suppose $f_{Y}(y)$ is a pdf defined on $(0,1)$, such as a
$Beta(a,b)$ distribution.

The symmetry of $g(\boldsymbol{\theta})$ allows for tidy transverse and
contour parameterizations via polar coordinates. Consider the transverse
parameter defined by the distance to the origin
$r(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}$
(i.e., the radius), and the contour parameter defined by the polar angle to
the positive $\theta_{1}$-axis,
$\phi(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\text{atan2}(y=\theta_{2},x=\theta_{1})$.
Observe that the polar angle is increasingly restricted as $r$ goes from 1 to
$\sqrt{2}$.

The connection between the data-space and transverse coordinates is
$y=g^{t}(r)\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{2}r^{2}$. Thus we
have $f_{R}(r)=f_{Y}(\frac{1}{2}r^{2})\cdot r$. Here, the choice of uniformly
distributed contours implies that

$\displaystyle f_{\Phi|R}(\phi|r)$ $\displaystyle=\left\\{\begin{array}[]{r
l}\frac{2}{\pi}\quad\mathbbm{1}\big{\\{}\ 0\ <\phi<\
\tfrac{\pi}{2}\big{\\}}&\quad 0<r\leq 1\\\ \frac{1}{\phi_{2}-\phi_{1}}\
\mathbbm{1}\big{\\{}\phi_{1}<\phi<\phi_{2}\big{\\}}&\quad
1<r\leq\sqrt{2}\end{array}\right.$ $\displaystyle\phi_{1}$
$\displaystyle=\text{atan2}\left(y=\sqrt{r^{2}-1},x=1\right)$
$\displaystyle\phi_{2}$
$\displaystyle=\text{atan2}\left(y=1,x=\sqrt{r^{2}-1}\right).$

Transforming $(R,\Phi)$ to $\boldsymbol{\Theta}$ comes with Jacobian
$\frac{1}{r}$ evaluated at $r=\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}$, and
therefore

$\displaystyle f^{\text{BBE}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=f_{Y}\Big{(}\tfrac{1}{2}(\theta_{1}^{2}+\theta_{2}^{2})\Big{)}\
f_{\Phi|R}\Big{(}\text{atan2}(\theta_{2},\theta_{1})\ |\
\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}\Big{)}\ .$

The bottom row of Figure 1 shows the solution with contours of the nonlinear
forward map (left). The observable density $f_{Y}(y)$ is a $Beta(8,12)$, and
this is plotted over the pushforward histogram of the BBE solution samples
(right).

##  5 Butler et al. 2018 (BJW) and Related Work

Butler, Jakeman, and Wildey (“BJW”) (Butler et al., 2018a) consider the SIP
(3) as above and call the use of their solution “consistent Bayesian
inference” or “pushforward based inference”. Here the use of the term
“consistent” does not refer to the statistical limiting sense, but rather to
the fact that the solution pushes forward to a given density. (We prefer
neither of these terms since we will show that this BJW solution is neither
Bayesian [Sections 5.2, 5.4] nor inference [Section 8]). The exact solution to
the SIP proposed by Butler et al. (2018a) is

$\displaystyle
f_{\boldsymbol{\Theta}}^{\text{BJW}}\big{(}\boldsymbol{\theta}\big{)}$
$\displaystyle=\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}\big{(}\boldsymbol{\theta}\big{)}\frac{f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}\
,$ (12)

where $\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}$ is a given density and
$\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}$ is its pushforward through
$\boldsymbol{g}$. The approximate solution
$\widehat{f}_{\boldsymbol{\Theta}}^{\text{BJW}}\big{(}\boldsymbol{\theta}\big{)}$
is what gets used in practice, and this has the same form except the density
in the denominator is replaced by an approximate density
$\widehat{\underaccent{\wtilde}{f}}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}$.

The BJW solution to the SIP was at least initially called Bayesian because
Bayes’ Rule was used in the derivation of the solution, and because the form
of the solution explicitly features a given initial density
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
(potentially playing a role similar to a prior distribution) times a weighting
function (akin to a likelihood function). Indeed, one must specify this
initial $p$-dimensional density
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$, similar
to the choice of a prior distribution during Bayesian inference. Unlike
standard Bayesian methods, however, the derivation also explicitly invokes the
disintegration theorem of measure theory, and the applied solution relies on
kernel density estimation for the denominator term. The practical reliance
upon density estimation restricts the application of BJW to very few
observable QoI (small $q$).

Subsequent works explore numerical aspects, applications, and extensions of
Butler et al. (2018a). Butler et al. (2018b) studies the convergence of kernel
density approximate solutions to the analytic BJW-derived densities. Walsh et
al. (2018) and Butler et al. (2020a) propose algorithms for optimal
experimental design which maximize the expected information gain between
initial and updated densities. Butler and Hakula (2020) applies the SIP
framework to a drum manufacturing process. The QoI are two (of twenty
possible, observable) eigenmodes of the drum vibration, and the parameters are
two diffusion parameters. Butler et al. (2020b) generalizes the BJW solution
to “stochastic” forward maps; we will discuss this generalization further in
Section 5.3 after an example in which (12) has closed form. Bruder et al.
(2020) uses multi-fidelity methods and Gaussian process regression models to
efficiently solve the SIP. Tran and Wildey (2021) focuses on a materials
science application, and also augments the BJW approach with a regression
model based on Gaussian processes to decrease the computational expense.
Finally, focusing on problems yielding large amounts of time-series data,
Mattis et al. (2022) learn the parameter-to-observable (QoI) map from data,
thus allowing for specification of the observable probability distribution.
Learning this map may rely on clustering the data and creating a partition on
the parameter-space. Once this is completed through the “Learning Uncertain
Quantities” framework, then the SIP can be solved as in previous works above.

###  5.1 Example: Linear Transformation of a Multivariate Gaussian Vector

Here we provide an analytic solution to a problem that was partly solved by
Butler et al. (2020a). Let
$\boldsymbol{g}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}\boldsymbol{\theta}$
so that $\boldsymbol{Y}=\boldsymbol{A}\boldsymbol{\Theta}$. Suppose that
$\boldsymbol{Y}\sim N_{q}(\boldsymbol{\mu}_{y},\boldsymbol{\Sigma}_{y})$ and
an initial density for $\boldsymbol{\Theta}$ is also multivariate Gaussian:
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})\sim
N_{p}(\boldsymbol{\mu}_{\theta},\boldsymbol{\Sigma}_{\theta})$. The
pushforward of $\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}$ under
$\boldsymbol{g}$ has distribution $\boldsymbol{\Gamma}\sim
N_{p}(\boldsymbol{A}\boldsymbol{\mu}_{\theta},\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})$.
Using these three given densities, the BJW solution (12) is

$\displaystyle f^{\text{BJW}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle\propto\exp\left\\{-(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta})^{\top}\boldsymbol{\Sigma}_{\theta}^{-1}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta})/2\right\\}$
$\displaystyle\ \ \cdot\
\frac{\exp\left\\{-(\boldsymbol{A}\boldsymbol{\theta}-\boldsymbol{\mu}_{y})^{\top}\boldsymbol{\Sigma}_{y}^{-1}(\boldsymbol{A}\boldsymbol{\theta}-\boldsymbol{\mu}_{y})/2\right\\}}{\exp\left\\{-(\boldsymbol{A}\boldsymbol{\theta}-\boldsymbol{A}\boldsymbol{\mu}_{\theta})^{\top}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})^{-1}(\boldsymbol{A}\boldsymbol{\theta}-\boldsymbol{A}\boldsymbol{\mu}_{\theta})/2\right\\}}$

After expanding, combining terms, and completing the quadratic form, it is
seen that the answer must be

$\displaystyle f^{\text{BJW}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=N_{p}(\wtilde{\boldsymbol{\mu}},\wtilde{\boldsymbol{\Sigma}})$
$\displaystyle\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\wtilde{\boldsymbol{\Sigma}}\left(\boldsymbol{A}^{\top}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{\mu}_{y}-\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})^{-1}\boldsymbol{A}\boldsymbol{\mu}_{\theta}+\boldsymbol{\Sigma}_{\theta}^{-1}\boldsymbol{\mu}_{\theta}\right)$
(13)
$\displaystyle=\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}(\boldsymbol{\mu}_{y}-\boldsymbol{A}\boldsymbol{\mu}_{\theta})+\boldsymbol{\mu}_{\theta}$
(14) $\displaystyle\wtilde{\boldsymbol{\Sigma}}$
$\displaystyle=\left(\boldsymbol{A}^{\top}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}-\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})^{-1}\boldsymbol{A}+\boldsymbol{\Sigma}_{\theta}^{-1}\right)^{-1}$
(15)
$\displaystyle=\boldsymbol{\Sigma}_{\theta}-\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T}\left\\{(\boldsymbol{\Sigma}_{y}^{-1}-(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T})^{-1})^{-1}+\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T}\right\\}^{-1}\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}$
(16)

Note that when $\boldsymbol{A}$ is a square invertible matrix,
$\wtilde{\boldsymbol{\Sigma}}=\boldsymbol{A}^{-1}\boldsymbol{\Sigma}_{y}\boldsymbol{A}^{-\top}$,
and $\wtilde{\boldsymbol{\mu}}=\boldsymbol{A}^{-1}\boldsymbol{\mu}_{y}$.
Otherwise, as shown in the appendix (Lemma 1),
$\boldsymbol{A}\wtilde{\boldsymbol{\mu}}=\boldsymbol{\mu}_{y}$ and
$\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{\top}=\boldsymbol{\Sigma}_{y}$,
confirming that the pushforward indeed has the distribution
$\boldsymbol{A}\boldsymbol{\Theta}\sim
N_{q}(\boldsymbol{\mu}_{y},\boldsymbol{\Sigma}_{y})$.

###  5.2 Sequential Updating(?)

With its apparent ability to update an initial distribution through (12), the
BJW solution has been compared to (and presented as an alternative to)
Bayesian inference. We show however that such comparisons are misleading.

Suppose that an analyst wants to compare two $q$-vectors
$\boldsymbol{Y}_{1},\boldsymbol{Y}_{2}$ of data to a model $\boldsymbol{g}$
that produces $q$ comparable outputs. The model takes $p$ inputs through the
vector $\boldsymbol{\theta}$, and the goal is to estimate these unknown
quantities. Because $\boldsymbol{g}$ is a well-defined function, it cannot
simultaneously map $\boldsymbol{\theta}$ to both $\boldsymbol{Y}_{1}$ and
$\boldsymbol{Y}_{2}$. The analyst therefore decides to use the BJW method to
update an initial distribution
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}$ before updating it once more
so that both pieces of information can be used. This procedure yields a final
answer

$\displaystyle\underaccent{\wtilde}{\ut
f}_{\boldsymbol{\Theta}}\big{(}\boldsymbol{\theta}\big{)}\frac{f_{\boldsymbol{Y}_{2}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{\ut
f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}$
$\displaystyle=\left(\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}\big{(}\boldsymbol{\theta}\big{)}\frac{f_{\boldsymbol{Y}_{1}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}\right)\frac{f_{\boldsymbol{Y}_{2}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{\ut
f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}$
$\displaystyle=\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}\big{(}\boldsymbol{\theta}\big{)}\frac{f_{\boldsymbol{Y}_{2}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}$

since $\underaccent{\wtilde}{\ut f}_{\boldsymbol{\Gamma}}\equiv
f_{\boldsymbol{Y}_{1}}$ (as $\underaccent{\wtilde}{\ut
f}_{\boldsymbol{\Theta}}$ solved the initial SIP). This is the BJW solution
using $\boldsymbol{Y}_{2}$, and as such, the final answer does not depend on
$\boldsymbol{Y}_{1}$! If done in reverse order, the final answer would not
depend upon $\boldsymbol{Y}_{2}$. In a Bayesian analysis, both observations
would be used. (Thankfully, after 2018, the terminology in the literature
related to BJW moved away from the name “consistent Bayes” and the analogy of
Bayesian inference.)

Thus sequential updating is not a defensible way to deal with _replicate_
measurements within an SIP. The next section explores a recent method that
allows for replicates.

###  5.3 The Extension of BJW to “Stochastic Maps” is Actually an Application
of BJW

In Butler et al. (2020c), the authors extend the BJW framework to so-called
“stochastic maps”—forward maps that include either embedded or additive
parameters meant to capture irreducible aleatoric uncertainty. In their terms,
if $\boldsymbol{g}(\boldsymbol{\theta})$ is a deterministic map, then the
corresponding stochastic map
$\hat{\boldsymbol{g}}(\boldsymbol{\theta},\boldsymbol{\epsilon})$ incoporates
either embedded or additive noise. In the additive case, we then have the
following system defining the SIP:

$\displaystyle\boldsymbol{Y}$
$\displaystyle=\hat{\boldsymbol{g}}(\boldsymbol{\Theta},\boldsymbol{E})\quad\quad\hat{\boldsymbol{g}}(\boldsymbol{\theta},\boldsymbol{\epsilon})=\boldsymbol{\theta}+\boldsymbol{\epsilon}$
(17)

But this is just a special case of BJW. That is, there is nothing inherently
stochastic in how the new $\boldsymbol{\epsilon}$ term interacts with
$\boldsymbol{g}$; if $\boldsymbol{g}(\cdot)$ is deterministic, then so is
$\hat{\boldsymbol{g}}(\boldsymbol{\theta},\boldsymbol{\epsilon})$. It is when
the parameters are turned into random variables
($\boldsymbol{\theta}\rightarrow\boldsymbol{\Theta}$ and
$\boldsymbol{\epsilon}\rightarrow\boldsymbol{E}$) that the response becomes
stochastic, and this has nothing to do with $\hat{\boldsymbol{g}}$.

The addition of these $\boldsymbol{\epsilon}$ ($\boldsymbol{E}$) parameters
simply accommodates the case when there are more observables than original
parameters of interest $\boldsymbol{\theta}$. As such, this redefined forward
map allows the BJW solution to be used when there are replicate measurements
on the same observable. The BJW solution for the system above gives the
“discrepancy” tradeoff between the $p+q$ random parameters
$\boldsymbol{\Theta}$ and $\boldsymbol{E}$ such that the marginal propagated
density $f_{\boldsymbol{Y}}$ is unaffected.

####  5.3.1 Gaussian Mean Estimation as “Stochastic Map” Inversion

To see the difference between the “stochastic map” BJW solution and a
statistical solution, consider the simplest possible scenario where the goal
is to estimate a single parameter of interest $\mu$ ($q=1$) from $n$
independent replicate samples. The statistical model is
$(Y_{i}|\mu)=\mu+E_{i}$ for
$E_{i}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}N(0,\sigma^{2}_{\epsilon})$
($i=1,\ldots,n$), or equivalently, $\boldsymbol{Y}\sim
N_{n}(\mu\boldsymbol{1}_{n},\sigma^{2}_{\epsilon}I_{n})$ for
$\sigma^{2}_{\epsilon}$ given. Within the BJW framework, the map is
$\hat{\boldsymbol{g}}(\mu,\boldsymbol{\epsilon})=\mu\boldsymbol{1}_{n}+\boldsymbol{\epsilon}$
so that $\boldsymbol{Y}=M\boldsymbol{1}_{n}+\boldsymbol{E}$ (where $M$ is the
random variable form of $\mu$).

The stochastic BJW solution comes from a system of $n$ equations in $n+1$
unknowns. Note that if any $E_{i}$ is fixed at a constant value or given a
known distribution, the system becomes square and $M$ is then determined by a
single observable. For example, if $E_{1}$ is set, then $Y_{1}$ completely
determines the density for $M=Y_{1}-E_{1}$; adding more observables
$Y_{n+1},\ldots$ only provides information on $\epsilon_{n+1},\ldots$. In
other words, the information on $M$ is contained within a single, arbitrary
$Y_{1}$ and does not increase as more information is gathered. However, if no
$E_{i}$ is set, then the BJW solution for the $n+1$ parameters instead
produces a distribution for $M$ that converges to a point mass.

The forward map is
$\boldsymbol{g}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}\boldsymbol{\theta}$
where
$\boldsymbol{A}\stackrel{{\scriptstyle\text{def}}}{{=}}[\boldsymbol{1}_{n}\vdots\boldsymbol{I}_{n}]$
and
$\boldsymbol{\theta}_{n}\stackrel{{\scriptstyle\text{def}}}{{=}}(\mu,\epsilon_{1},\ldots,\epsilon_{n})^{\top}$.
The observable distribution is $\boldsymbol{Y}\sim
N_{n}(\boldsymbol{\mu}_{y},\boldsymbol{\Sigma}_{y})$ with
$\boldsymbol{\Sigma}_{y}=\sigma^{2}_{y}\boldsymbol{I}_{n}$. When the initial
distribution is $\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}\equiv
N_{n+1}(\boldsymbol{\mu}_{\theta},\boldsymbol{\Sigma}_{\theta})$ with
$\boldsymbol{\mu}_{\theta}^{\top}=(\mu_{0},\boldsymbol{0}_{n}^{\top})$ and
$\boldsymbol{\Sigma}_{\theta}=\text{diag}(\sigma_{0}^{2},\sigma^{2}_{\epsilon},\ldots,\sigma^{2}_{\epsilon})$,
the BJW solution is
$N_{n+1}(\wtilde{\boldsymbol{\mu}},\wtilde{\boldsymbol{\Sigma}})$ with mean
and covariance given by (13, 15). We derive the closed-form moments in Section
B.2 of the appendix. The matrix algebra is tedious but reveals some
interesting features of the BJW solution, namely that the mean and variance
take the form

$\displaystyle\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\begin{bmatrix}\overline{\mu}_{y}+\mathcal{O}(n^{-1})\\\
(\boldsymbol{\mu}_{y}-\overline{\mu}_{y}\boldsymbol{1}_{n})+\mathcal{O}(n^{-1})\boldsymbol{1}_{n}\end{bmatrix}\quad\quad\wtilde{\boldsymbol{\Sigma}}=\begin{bmatrix}\mathcal{O}(n^{-1})&-\mathcal{O}(n^{-1})\boldsymbol{1}_{n}^{\top}\\\
-\mathcal{O}(n^{-1})\boldsymbol{1}_{n}&\sigma^{2}_{y}\boldsymbol{I}_{n}-\mathcal{O}(n^{-1})\boldsymbol{J}_{n}\end{bmatrix}\
.$ (18)

The marginal variance of $M$ decreases like $\mathcal{O}(n^{-1})$, implying
that the distribution of $M$ concentrates on $\overline{\mu}_{y}$, a simple
average of the means of $\boldsymbol{Y}$. Marginally, $\boldsymbol{E}$
converges to a distribution with covariance $\sigma^{2}_{y}\boldsymbol{I}_{n}$
that is centered on a residual vector. Thus all the irreducible uncertainty is
contained within $\boldsymbol{E}$, but not within the parameter of interest,
contradicting the intended purpose of the BJW approach.

For any $n$, $\wtilde{\boldsymbol{\Sigma}}$ is a dense matrix, meaning that
$M$ and $\boldsymbol{E}$ covary and all elements of $\boldsymbol{E}$ covary
with one another. Furthermore, the mean of $M$ is a linear combination of
$\overline{\mu}_{y}$ and $\mu_{0}$. While this is reminiscent of Bayesian
inference, it should be noted that the terms in the BJW solution are
considerably more complicated and thus harder to interpret. In addition, a
statistical solution involves simple scalar arithmetic while the BJW solution
comes from inverting and multiplying potentially large matrices.

###  5.4 The BJW Solution to the SIP is a Change-of-Variables

Due to its generality, the use of measure theory in describing and solving
problems is often preferable. However, if restricting one’s attention to
absolutely continuous random variables (A2) in practice, then the language of
measure theory can obscure concepts which are quite simple in nature. In the
next proposition we show that the BJW density can be derived in a
straightforward manner under Assumption (A1). The second half of its proof
relies upon the notion of _auxiliary variables_. For some choice (but
typically, infinitely many choices) of $p-q$ auxiliary variables
$\boldsymbol{Y}^{c}$ via maps
$\boldsymbol{g}^{c}(\cdot)\stackrel{{\scriptstyle\text{def}}}{{=}}\left(g_{q+1}(\cdot),\ldots,g_{p}(\cdot)\right)$,
the following augmented system is locally invertible:

$\displaystyle\boldsymbol{\Theta}=\left[\begin{array}[]{l}\Theta_{1}\\\
\vdots\\\ \Theta_{q}\\\\[3.87498pt] \Theta_{q+1}\\\\[3.87498pt] \vdots\\\
\Theta_{p}\end{array}\right]\quad\begin{array}[]{c}\boldsymbol{g}_{\text{aug}}\\\
\longmapsto\end{array}\quad\left[\begin{array}[]{lcl}Y_{1}&=&g_{1}(\boldsymbol{\Theta})\\\
\vdots&&\\\ Y_{q}&=&g_{q}(\boldsymbol{\Theta})\\\\[3.87498pt] \hline\cr
Y_{q+1}&\stackrel{{\scriptstyle\text{def}}}{{=}}&g_{q+1}(\boldsymbol{\Theta})\\\
\vdots&&\\\
Y_{p}&\stackrel{{\scriptstyle\text{def}}}{{=}}&g_{p}(\boldsymbol{\Theta})\end{array}\right]=\begin{bmatrix}\boldsymbol{Y}\\\
\boldsymbol{Y}^{c}\end{bmatrix}\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{Y}_{\text{aug}}\
.$ (33)

For the time being we shall take the existence of
$\boldsymbol{g}_{\text{aug}}$ defining the locally invertible system as given,
but more details will be given in Section 7.

######  Proposition 4.

Under assumption (A1) the exact BJW solution can be derived from a CoV
solution for $p\geq q$.

###### Proof.

Fix some initial $p$-dimensional density
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$.

First consider the case that $p=q$. Using the same reasoning leading to (7),
the initial density can be written in terms of its pushforward on
$\mathcal{P}\setminus\mathcal{P}_{0}$, regardless of the global invertibility
of $\boldsymbol{g}$. Equivalently, the Jacobian determinant is

$\displaystyle\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}$
$\displaystyle=\frac{\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}(\boldsymbol{g}(\boldsymbol{\theta}))}$
(34)

wherever $\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}>0$ (and hence
$\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}>0$) on
$\mathcal{P}\setminus\mathcal{P}_{0}$. Now we have

$\displaystyle
f_{\boldsymbol{\Theta}}^{\text{BJW}}\big{(}\boldsymbol{\theta}\big{)}$
$\displaystyle=f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}\frac{\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}\big{(}\boldsymbol{\theta}\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}=f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}\
\Big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\Big{|}$

which is the density corresponding to a CoV, as desired. The explicit mixture
of CoVs from $\mathcal{Q}$ to $\mathcal{P}$ can be constructed exactly as in
the proof of Proposition 2. Note that in this $p=q$ case,
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}$ did not affect the solution
to the SIP since
$\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\big{|}$ is
invariant to this choice.

The case of $p>q$ follows exactly the same reasoning as above but relies on
auxiliary variables, as per (33). Additionally, it is slightly more
transparent to start with $f^{\text{CoV}}$, and show that for a certain
choice, it becomes the exact BJW solution. Replacing “$\boldsymbol{g}$”,
“$\boldsymbol{Y}$”, and “$\boldsymbol{\Gamma}$” above with
“$\boldsymbol{g}_{\text{aug}}$”, “$\boldsymbol{Y}_{\text{aug}}$”, and
“$\boldsymbol{\Gamma}_{\text{aug}}$”, we have by previous reasoning that

$\displaystyle\Big{|}\frac{\partial\boldsymbol{g}_{\text{aug}}}{\partial\boldsymbol{\theta}}\Big{|}$
$\displaystyle=\frac{\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}_{\text{aug}}}(\boldsymbol{g}_{\text{aug}}(\boldsymbol{\theta}))}$
$\displaystyle f^{\text{CoV}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=f_{\boldsymbol{Y}_{\text{aug}}}\big{(}\boldsymbol{g}_{\text{aug}}(\boldsymbol{\theta})\big{)}\
\Big{|}\frac{\partial\boldsymbol{g}_{\text{aug}}}{\partial\boldsymbol{\theta}}\Big{|}$
$\displaystyle=\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})\
\frac{f_{\boldsymbol{Y}_{\text{aug}}}\big{(}\boldsymbol{g}_{\text{aug}}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}_{\text{aug}}}\big{(}\boldsymbol{g}_{\text{aug}}(\boldsymbol{\theta})\big{)}}$
$\displaystyle=\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})\frac{f_{\boldsymbol{Y}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}(\boldsymbol{\theta})\big{)}}\
\frac{f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}\big{(}\boldsymbol{g}^{c}(\boldsymbol{\theta})\
|\
\boldsymbol{y}=\boldsymbol{g}(\boldsymbol{\theta})\big{)}}{\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}}\big{(}\boldsymbol{g}^{c}(\boldsymbol{\theta})\
|\ \boldsymbol{\gamma}=\boldsymbol{g}(\boldsymbol{\theta})\big{)}}\ ,$

where
$\boldsymbol{\Gamma}_{\text{aug}}=\\{\boldsymbol{\Gamma},\boldsymbol{\Gamma}^{c}\\}$
and $\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}_{\text{aug}}}$ is the
pushforward of $\underaccent{\wtilde}{f}$ through
$\boldsymbol{g}_{\text{aug}}$. Setting
$f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}\equiv\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}}$
leads to cancellation that establishes that
$f^{\text{BJW}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ is indeed
derivable from a CoV solution. ∎

In the proof of Proposition above, Bayes’ Theorem was not invoked and the
disintegration theorem did not need to be called explicitly. All that was
required was that the joint density of unobserved and observed variables
$\boldsymbol{Y}_{\text{aug}}=\\{\boldsymbol{Y},\boldsymbol{Y}^{c}\\}$ could be
factored as a conditional times a marginal density. Of course this fact is
directly related to the disintegration theorem, but will nonetheless be easier
for most audiences to digest.

The statement at the end of the proof,
$f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}\stackrel{{\scriptstyle\text{set}}}{{=}}\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}}$,
may seem mysterious, but actually provides the density to assign to the set of
auxiliary variables within a CoV exercise to match the BJW solution. First
choose a complementary map $\boldsymbol{g}^{c}(\cdot)$ that yields a full-rank
augmented system (a.e.). Under the augmented map
$\boldsymbol{g}_{\text{aug}}$, the initial density
$\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}$ pushes forward to the
$p$-dimensional joint density
$\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma},\boldsymbol{\Gamma}^{c}}=\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}}\cdot\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}}$.
The BJW solution to the SIP is thus a density obtained through a linear
combination of CoVs from $\mathcal{Q}$ to $\mathcal{P}$ with
$\boldsymbol{Y}_{\text{aug}}\sim
f_{\boldsymbol{Y}}(\boldsymbol{y})\cdot\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}}(\boldsymbol{y}^{c}|\boldsymbol{y})$.

In theory, one could sample from
$f^{\text{BJW}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ using Monte Carlo
reasoning akin to Algorithm 1. One would first generate
$\boldsymbol{y}^{*}\sim f_{\boldsymbol{Y}}$, then use this value to obtain a
conditional draw
$\boldsymbol{y}^{c*}\sim\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}=\boldsymbol{y}^{*}}$
to form
$\boldsymbol{y}_{\text{aug}}^{*}\stackrel{{\scriptstyle\text{def}}}{{=}}(\boldsymbol{y}^{*},\boldsymbol{y}^{c*})$.
Finally, the solution to
$\boldsymbol{y}_{\text{aug}}^{*}-\boldsymbol{g}_{\text{aug}}\big{(}\boldsymbol{\theta}^{*}\big{)}=\boldsymbol{0}$
would give a realization $\boldsymbol{\theta}^{*}$ of the BJW solution.
However, this is hypothetical since in practice one does not have the ability
to simulate from the unknown conditional density
$\underaccent{\wtilde}{f}_{\boldsymbol{\Gamma}^{c}|\boldsymbol{\Gamma}}$.

##  6 Other Related Work

In order to complete our literature review of SIPs, we need to discuss two
more recent works: Swigon et al. (2019) and Uy and Grigoriu (2019).

Swigon et al. (2019) demonstrates the critical role of the Jacobian
determinant in two types of parameter estimation problems. In the first, the
data is given by a nonlinear transformation of the parameters, i.e.,
$\boldsymbol{Y}=\boldsymbol{g}(\boldsymbol{\Theta})$, and the authors call
this situation (i.e., our SIP) a “random parameter model”. Here the answer is
simply given by the CoV, similar to our conclusion for $p=q$, and as such
depends on the Jacobian determinant. In the second type of parameter
estimation problem, the data is corrupted by error, i.e.,
$\boldsymbol{Y}=\boldsymbol{g}(\boldsymbol{\theta})+\boldsymbol{E}$, for some
fixed $\boldsymbol{\theta}$, which the authors called a “random measurement
error model” (and we call a statistical IP). In this context the Jacobian
determinant appears in the construction of a default prior for
$\boldsymbol{\Theta}$. The authors take both these modeling approaches as
worthy of consideration, and while offering some discussion akin to our
Section 8, do not tie their results into the work stemming from BBE and BJW.

Uy and Grigoriu (2019) restricts attention to the case that $p>q$ and approach
SIPs from a different angle: given that the problem is underdetermined, what
additional knowledge about $\boldsymbol{\Theta}$ is required to pose and solve
SIPs? (By contrast, traditional IPs can use “extra” information to regularize
solutions in overdetermined systems, but here represents necessary missing
information for the underdetermined system.) In short, one can either specify
moments of $\boldsymbol{\Theta}$ or its distributional family, similar to
methods based on the principle of maximum entropy. Interestingly, if
additional information is provided on the family of distributions to which
$\boldsymbol{\Theta}$ belongs, then the new inverse problem can be solved
using standard Bayesian methods (see Section 3.2.2). For each option, the
authors explore two possibilities about the data: (1) pdf information about
$\boldsymbol{Y}=\boldsymbol{g}(\boldsymbol{\Theta})$ is known, or (2) samples
from $\boldsymbol{Y}$ are given. While mainly focusing on the BBE solutions,
Uy and Grigoriu (2019) stresses the fundamental underdetermined nature of SIPs
(previously noted in Breidt et al. (2011)). We will take this a little further
in Section 8 (IV).

Uy and Grigoriu (2019) shows that the uniform contour ansatz can lead to
arbitrarily bad approximations of the true density, and in some cases cannot
recover the true density at all. Moreover, the authors specify this additional
information in order to find the true distribution of $\boldsymbol{\Theta}$,
so that it can be queried in other contexts, and not just to match
$\boldsymbol{g}(\boldsymbol{\Theta})$. This goal causes a philosophical breach
from Breidt et al. (2011) and Butler et al. (2018a), in which different data
points may lead to different solutions, even within the same problem.

##  7 All Solutions to SIPs are Changes-of-Variables

Our main proposition follows from the CoV theorem and contains elements of
Proposition 4. We give the proof here because it is both constructive and
instructive.

######  Proposition 5.

Let $f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ be any exact solution to
the SIP of (3), and suppose the forward map $\boldsymbol{g}$ is as described
in Assumption (A1). Then $f$ is derivable from a CoV (a.e.) for at least one
choice of auxiliary variables.

###### Proof.

The $p=q$ case is covered in Proposition 2, so we now assume $p>q$. After
constructing auxiliary variables, the proof follows in much the same fashion.

Under (A1) we have that the component functions of $\boldsymbol{g}$ (i.e.,
$g_{1},\ldots,g_{q}$) produce a Jacobian with left $q\times q$ block
$\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}_{1:q}}\big{|}$
that is invertible on $\mathcal{P}_{0}$. Now observe that there exist maps
$g_{q+1},\ldots,g_{p}$ from $\left(\mathcal{P}\cap\mathbb{R}^{p-q}\right)$ to
$\mathbb{R}^{p-q}$ such that the corresponding lower-right $(p-q)\times(p-q)$
sub-Jacobian is invertible; Again, let $\boldsymbol{g}_{\text{aug}}$ denote
the full $p$-dimensional function. Any choice of these $(p-q)$ functions
results in the definition of auxiliary variables
$\boldsymbol{Y}^{c}\stackrel{{\scriptstyle\text{def}}}{{=}}(Y_{q+1},\ldots,Y_{p})^{\top}$.
The simplest augmented system is one that uses the identity maps to define the
components of $\boldsymbol{Y}^{c}$

$\displaystyle\left[\begin{array}[]{l}\Theta_{1}\\\ \vdots\\\
\Theta_{q}\\\\[3.87498pt] \Theta_{q+1}\\\\[3.87498pt] \vdots\\\
\Theta_{p}\end{array}\right]\begin{array}[]{c}\boldsymbol{g}_{\text{aug}}\\\
\longmapsto\end{array}\left[\begin{array}[]{lcl}Y_{1}&=&g_{1}(\boldsymbol{\Theta})\\\
\vdots&&\\\ Y_{q}&=&g_{q}(\boldsymbol{\Theta})\\\\[3.87498pt] \hline\cr
Y_{q+1}&\stackrel{{\scriptstyle\text{def}}}{{=}}&g_{q+1}(\boldsymbol{\Theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\Theta_{q+1}\\\
\vdots&&\\\
Y_{p}&\stackrel{{\scriptstyle\text{def}}}{{=}}&g_{p}(\boldsymbol{\Theta})\hskip
11.0pt\stackrel{{\scriptstyle\text{def}}}{{=}}\Theta_{p}\end{array}\right]$
(49)

which produces Jacobian determinant

$\displaystyle\left|\frac{\partial\boldsymbol{g}_{\text{aug}}}{\partial\boldsymbol{\theta}}\right|\quad$
$\displaystyle=\quad\text{det}\begin{bmatrix}\frac{\partial
g_{1}}{\partial\theta_{1}}&\cdots&\frac{\partial
g_{1}}{\partial\theta_{q}}&\frac{\partial
g_{1}}{\partial\theta_{q+1}}&\cdots&\frac{\partial
g_{1}}{\partial\theta_{p}}\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\frac{\partial g_{q}}{\partial\theta_{1}}&\cdots&\frac{\partial
g_{q}}{\partial\theta_{q}}&\frac{\partial
g_{q}}{\partial\theta_{q+1}}&\cdots&\frac{\partial
g_{q}}{\partial\theta_{p}}\\\\[3.87498pt] \hline\cr
0&\cdots&0&1&\boldsymbol{0}^{\top}&0\\\
\vdots&\vdots&\vdots&\boldsymbol{0}&\ddots&\boldsymbol{0}\\\
0&\cdots&0&0&\boldsymbol{0}^{\top}&1\\\
\end{bmatrix}\quad=\quad\left|\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}_{1:q}}\right|\
.$ (50)

This augmented system is invertible due to the fact that the upper-left
$(q\times q)$ block
$\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}_{1:q}}$ is itself
invertible. The pushforward of $f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
is thus possible by the CoV theorem. Furthermore, because
$f_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$ is assumed to solve the SIP,
its pushforward is a $p$-dimensional distribution
$f_{\boldsymbol{Y}_{\text{aug}}}$ whose marginal is
$f_{\boldsymbol{Y}}(\boldsymbol{y})$ for the first $q$ entries. In fact,
$f_{\boldsymbol{Y}_{\text{aug}}}=f_{\boldsymbol{Y}}\cdot
f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}$ with

$\displaystyle
f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}(\boldsymbol{y}^{c}|\boldsymbol{y})\stackrel{{\scriptstyle\text{def}}}{{=}}f_{\boldsymbol{\Theta}_{(q+1):p}|\boldsymbol{\Theta}_{1:q}}(\boldsymbol{y}^{c}|\boldsymbol{y})$

under the identity maps.

Going the other direction, the reasoning in the proof of Proposition 2 can
again be used. Specifically, Assumption (A1) guarantees local inverses
$\boldsymbol{h}_{1},\ldots,\boldsymbol{h}_{d}$ of $\boldsymbol{g}$ defined on
the appropriate sets, and we have that the posited solution can be written as
a particular $d$-combination of CoV solutions from the augmented range to
$\mathcal{P}\setminus\mathcal{P}_{0}$ except for on a set of measure zero. The
presumed solution was then actually a CoV using a particular choice of
auxiliary variables $\boldsymbol{Y}^{c}$.

Thus we have shown that no matter how one has obtained an exact solution to
the SIP, it could have actually been derived from CoV solutions for any choice
of auxiliary variables to make an augmented forward map
$\boldsymbol{g}_{\text{aug}}$ locally invertible. ∎

Some readers will recognize the augmented system within the proof as that used
within standard constructive proof of the Implicit Function Theorem (given the
Inverse Function Theorem); see e.g. Lee (2013) pg. 661. This auxiliary
variable strategy is exactly what is taught in a first course on mathematical
statistics. For example, to find the distribution of the product $Y$ of two
random variables $\Theta_{1}$ and $\Theta_{2}$, one may augment the one-
dimensional system
$Y\stackrel{{\scriptstyle\text{def}}}{{=}}\Theta_{1}\Theta_{2}$ with the
auxiliary variable $Y^{c}\stackrel{{\scriptstyle\text{def}}}{{=}}\Theta_{2}$
(for example), obtain $f_{Y,Y^{c}}$ from the two-dimensional CoV, and
integrate the joint density to get $f_{Y}$.

The proposition above shows that the answer will be completely determined
(a.e.) from the choice of $p-q$ dimensional density
$f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}$. Intuitive solutions (9) can be viewed
as the simplest cases that use $(p-q)$ identity maps (as in the proof above),
together with $f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}}\equiv
f_{\boldsymbol{\Theta}_{(q+1):p}}$ chosen by the analyst.

In order to sample from a CoV density
$f^{\text{CoV}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$, it may appear
that the Jacobian determinant
$\big{|}\frac{\partial\boldsymbol{g}_{\text{aug}}}{\partial\boldsymbol{\theta}}\big{|}$
(or at the very least
$\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}_{1:q}}\big{|}$)
is readily available (Swigon et al., 2019). However, this is not true. Instead
of using the explicit density to obtain random draws, one could instead use a
simple Monte Carlo routine similar to Algorithm 1. One would first generate
$\boldsymbol{y}^{*}\sim f_{\boldsymbol{Y}}$, then use this value to obtain a
conditional draw $\boldsymbol{y}^{c*}\sim
f_{\boldsymbol{Y}^{c}|\boldsymbol{Y}=\boldsymbol{y}^{*}}$. The vector
$\boldsymbol{y}_{\text{aug}}^{*}\stackrel{{\scriptstyle\text{def}}}{{=}}(\boldsymbol{y}^{*},\boldsymbol{y}^{c*})$
is a random draw from $f_{\boldsymbol{Y}_{\text{aug}}}$. Finally, the solution
to
$\boldsymbol{y}_{\text{aug}}^{*}-\boldsymbol{g}_{\text{aug}}\big{(}\boldsymbol{\theta}^{*}\big{)}=\boldsymbol{0}$
would give a realization $\boldsymbol{\theta}^{*}$ of the CoV solution.

###  7.1 Example: Linear Transformation of a Multivariate Gaussian Vector

Let us return to the example of Section 5.1. Again, let
$\boldsymbol{g}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}\boldsymbol{\theta}$
so that $\boldsymbol{Y}=\boldsymbol{A}\boldsymbol{\Theta}$; the dimension of
$\boldsymbol{A}$ is $q\times p$, and this matrix is assumed to have full row-
rank. Suppose that the observation distribution is $\boldsymbol{Y}\sim
N_{q}(\boldsymbol{\mu}_{y},\boldsymbol{\Sigma}_{y})$.

When $p=q$, $\boldsymbol{A}$ is invertible so the SIP is solved by
$\boldsymbol{\Theta}=\boldsymbol{A}^{-1}\boldsymbol{Y}$, and this has density

$\displaystyle
f^{\text{CoV}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})=N_{p}\Big{(}\boldsymbol{A}^{-1}\boldsymbol{\mu}\
,\ \boldsymbol{A}^{-1}\boldsymbol{\Sigma}\boldsymbol{A}^{-\top}\Big{)}\ .$
(51)

When $p>q$, $\boldsymbol{A}$ can be augmented into an invertible matrix in an
infinite number of ways. Suppose that $\boldsymbol{A}^{\perp}$ is a
$(p-q)\times p$ representation of the null space of row($\boldsymbol{A}$). The
augmented matrix
$\boldsymbol{A}_{\text{aug}}\stackrel{{\scriptstyle\text{def}}}{{=}}\begin{bmatrix}\boldsymbol{A}\\\
\boldsymbol{A}^{\perp}\end{bmatrix}$ is then invertible. Furthermore
$\boldsymbol{A}^{\perp}$ defines $p-q$ auxiliary variables
$\boldsymbol{Y}^{c}$. If the joint distribution of
$(\boldsymbol{Y},\boldsymbol{Y}^{c})$ is
$N_{p}(\boldsymbol{\mu}_{+},\boldsymbol{\Sigma}_{+})$ (making
$\boldsymbol{\mu}_{y}$ the first $q$ entries of $\boldsymbol{\mu}_{+}$ and
$\boldsymbol{\Sigma}_{y}$ the upper left $q\times q$ block of
$\boldsymbol{\Sigma}_{+}$), then the CoV solution to the SIP is

$\displaystyle f^{\text{CoV}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})$
$\displaystyle=N_{p}\Big{(}\boldsymbol{A}_{\text{aug}}^{-1}\boldsymbol{\mu}_{+}\
;\
\boldsymbol{A}_{\text{aug}}^{-1}\boldsymbol{\Sigma}_{+}\boldsymbol{A}_{\text{aug}}^{-1}\Big{)}\
.$

The density above begs the question: how does one choose a density for
_unobservable_ quantities in order to augment and close the system? We return
to this fundamental question in the next section. The issue of
underdeterminedness was observed back in Breidt et al. (2011) (pg. 1839) for
the bivariate Gaussian case ($p=2,q=1$) by considering moment conditions
(instead of auxiliary variables), but this did not curtail the investigation
or application of SIPs.

##  8 Changes-of-Variables and Inference

In the previous section we showed that any solution to an SIP (3) can be
viewed as a CoV. The solution to the IP given in (2) is typically a matter of
statistical _inference_. What is the connection between these two solutions?
In this section we explore this question and more broadly examine the
interface of CoVs and inference. While multiple works provide some discussion
to distinguish SIP solutions from more traditional statistical
inference—Section 2 in Breidt et al. (2011); Section 7 in Butler et al.
(2018a); Section 4.3 in Butler et al. (2014); Remark 2.3 in Butler et al.
(2020c); and throughout Uy and Grigoriu (2019), especially Sections 1 and
4.2—there is still much to address, and we do so here.

In any inferential context, a change of the likelihood function $\propto
f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{\theta})$ with respect to
$\boldsymbol{y}$ is obviously a CoV. In the Bayesian or Fiducial paradigms, a
change of the target distribution $f(\boldsymbol{\theta}|\boldsymbol{y})$ with
respect to $\boldsymbol{\theta}$ is a reparameterization that is also quite
conspicuously a CoV. A less obvious connection between inference and CoVs is
the very concept of Bayesian inference itself. If
$f_{\boldsymbol{\Theta}}^{\text{pri}}(\boldsymbol{\theta})$ is a proper prior
distribution, and
$f_{\boldsymbol{\Theta}}^{\text{post}}(\boldsymbol{\theta})\propto
f(\boldsymbol{y}|\boldsymbol{\theta})\cdot
f_{\boldsymbol{\Theta}}^{\text{pri}}(\boldsymbol{\theta})$ denotes the
posterior distribution, there is some “forward map” $\boldsymbol{g}_{T}$ that
pushes prior forward to posterior. In fact this is the motivation behind
_transport maps_ where, in an inference setting, the goal is to learn exactly
that function $\boldsymbol{g}_{T}$. This CoV is however entirely distinct from
the SIP solutions discussed here: with transport maps, the function
$\boldsymbol{g}_{T}$ maps _from parameter to parameter_ space, not from
parameter to observable space (Marzouk et al., 2016; Baptista et al., 2021).

The solutions to (2) and (3) can sometimes agree, but as we will argue, this
does not mean that statistical IPs and SIPs are analogous notions. We will
first give examples where Bayesian posterior and Generalized Fiducial (GF)
distributions are identical to a CoV before giving five (interrelated) ways in
which SIPs are fundamentally different from inference. We close by providing
two illustrative examples.

Bayesian and GF solutions can appear to be CoVs when $p=q$. Consider again the
linear map to a Gaussian observable (Section 7.1) having known covariance
$\boldsymbol{\Sigma}_{y}$. When $p=q$ and an improper prior is chosen
$f_{\boldsymbol{\Theta}}^{\text{pri}}(\boldsymbol{\theta})\propto 1$ for
$\boldsymbol{\Theta}$, the posterior is
$f_{\boldsymbol{\Theta}}^{\text{post}}(\boldsymbol{\theta})\propto
f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{A}\boldsymbol{\theta})\cdot 1$
which agrees with the CoV solution (51) after a simple rearrangement. A
similar result is also obtained in the appendix of Swigon et al. (2019). The
density above is also the uncertainty distribution in the GF paradigm (Hannig
et al., 2016, Ex. 3). In fact, when $p=q$ and $\boldsymbol{g}$ is an
invertible nonlinear map, the GF solution is
$f_{\boldsymbol{\Theta}}^{\text{GF}}(\boldsymbol{\theta})\propto
f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{g}(\boldsymbol{\theta}))\big{|}\frac{\partial\boldsymbol{g}}{\partial\boldsymbol{\theta}}\big{|}$
(Hannig et al., 2016, Eqns 3,4) which is the CoV solution (5). The agreement
in these special cases is superficial, as a closer look in the first point
below will show a sharp contrast in the interpretation of the terms involved.

(I) In SIPs, observables are _populations_ with given parameters. It is in the
application of the CoV methodology that one can clearly see the difference
between SIPs and any inferential framework based upon a likelihood function.
In an inferential context, after the data is observed, $\boldsymbol{Y}$
becomes $\boldsymbol{y}_{\text{obs}}$ and the data pdf
$f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{\theta})$ will yield the
likelihood $L(\boldsymbol{\theta};\boldsymbol{y}_{\text{obs}})$. The
likelihood is a function of $\boldsymbol{\theta}$ and indexed by the data
$\boldsymbol{y}_{\text{obs}}$. In the Gaussian examples above, the _unknown_
mean of $\boldsymbol{Y}$ features the parameters:
$\boldsymbol{A}\boldsymbol{\theta}$ or $\boldsymbol{g}(\boldsymbol{\theta})$.
On the other hand, the CoV solution is based upon the density for the
observables $f_{\boldsymbol{Y}}(\boldsymbol{y})$ which is actually
$f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{y}_{\text{obs}})$ where the
observed “data” actually operate as _known_ population parameters! This
density has both $\boldsymbol{y}$ and $\boldsymbol{y}_{\text{obs}}$ as
arguments, and the solution will have the $\boldsymbol{y}$ argument replaced
with the forward map:
$f_{\boldsymbol{Y}}(\boldsymbol{y}=\boldsymbol{g}(\boldsymbol{\theta});\
\boldsymbol{y}_{\text{obs}})$. The agreement of the Bayesian/GF and CoV
solutions for the $p=q$ cases above is due to the symmetry of the Gaussian
density with respect to its $\boldsymbol{y}$ and mean arguments. Moreover,
SIPs require all such observable population parameters to be given: as seen
above, SIPs involving Gaussian observables require the covariance
$\boldsymbol{\Sigma}_{y}$ to be given. This is quite different from
statistical IPs where quite often covariances are parameters to estimate.

SIP solutions are not inference because they essentially treat everything on a
“population” level to begin with: the observables _are not a realization_ from
an unknown population density, they _are the population_ with a known density.
In the case of inference, the data model
$f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{g}(\boldsymbol{\theta}))$
features the forward map and population parameters $\boldsymbol{\theta}$ to
infer. The likelihood function gives desirable large-sample properties such as
the consistency and asymptotic normality of the maximum likelihood estimator
(Reid, 2010). In the case of SIP solutions, the density for the observables
$f_{\boldsymbol{Y}}(\boldsymbol{y};\boldsymbol{y}_{\text{obs}})$ is a
population-level description with known parameters and depends upon neither
$\boldsymbol{g}$ nor $\boldsymbol{\theta}$. As such,
$f_{\boldsymbol{Y}}(\boldsymbol{y})$ does not feature any
$\boldsymbol{g}(\boldsymbol{\theta})$ in it, so can neither be interpreted as
a conditional density nor a likelihood given the map
$\boldsymbol{g}(\boldsymbol{\theta})$. Hence there is no likelihood function
or the corresponding theory.

In fact, the phrases “increasing sample size” and “collecting more data” are
not really even coherent phrases in the world of SIPs; adding observables is
simply changing the population to invert. Despite the fact that the BJW
solution was first presented as “consistent Bayesian inference”, it too
suffers from from this changing-populations issue because none of the
$f_{\boldsymbol{Y}}$’s for the new data explicitly contain any
$\boldsymbol{g}(\boldsymbol{\theta})$ terms. In contrast, within any
statistical inference procedure, each new piece of data will contribute
knowledge of $\boldsymbol{\theta}$ via $\boldsymbol{g}(\boldsymbol{\theta})$
through a joint likelihood function.

(II) SIPs require $p\geq q$. Perhaps the most obvious sign that an SIP is
different from a statistical IP is the hard stipulation that $p\geq q$, i.e.,
more parameters than observables or data. On the other hand, in any
likelihood-based inference setting (Bayesian, GF, frequentist, etc.), an
analyst takes great care to avoid unregularized or saturated models and the
inherent risk of “over-fitting”. This is tied to the problem of _prediction_.

Along with obtaining an uncertainty distribution for $\boldsymbol{\Theta}$,
another common goal is prediction: inferring a new value given those values
already observed. In the Bayesian/GF frameworks, the (posterior) predictive
density is $f_{\boldsymbol{Y}^{*}|\boldsymbol{Y}}=\int
f_{\boldsymbol{Y}^{*}|\boldsymbol{\Theta},\boldsymbol{Y}}\
f_{\boldsymbol{\Theta}|\boldsymbol{Y}}d\boldsymbol{\theta}$. This density is
trustworthy when the models $f_{\boldsymbol{Y}|\boldsymbol{\Theta}}$ and
$f_{\boldsymbol{Y}^{*}|\boldsymbol{\Theta},\boldsymbol{Y}}$ are adequate.
Adequacy can be assessed by analyzing residuals and out-of-sample predictions.
When $p\geq q$, a non-Bayesian introduces regularization or smoothness
penalties (such as in LASSO) to better constrain the problem, and a Bayesian
codifies complexity penalization into the prior specification. In SIPs there
is no way to accommodate such regularization: all solutions embody over-
fitting. Prediction, however, is still possible within the CoV framework,
through some new forward map $\boldsymbol{g}^{*}$. After obtaining
$\boldsymbol{\Theta}\sim f_{\boldsymbol{\Theta}}^{\text{CoV}}$, a predictive
density is simply $\boldsymbol{g}^{*}(\boldsymbol{\Theta})$. However, as
$\boldsymbol{Y}\rightarrow\boldsymbol{\Theta}\rightarrow\boldsymbol{g}^{*}(\boldsymbol{\Theta})$
is a sequence of population-level reparameterizations (I), there are no
falsifiable models that can be checked, and hence no way to defensibly justify
prediction.

The SIP framework has no extension into the realm of $p<q$; instead, one has
to force a problem into becoming an SIP and out of the realm of inference.
Approaches to ensure that $p\geq q$ include the following. First, the analyst
can pick a new subset of $p$ or fewer observables, even though each subset
will lead to a different answer. This is the approach followed in Butler et
al. (2015a), in which 2 of 12 possible water levels are selected, and in
Butler and Hakula (2020), in which 2 of 20 possible eigenmodes of a drum are
selected. A related approach is to combine the $q$ original observables into
$p$ or fewer, such as by taking averages and pooling variances. Per the LUQ
framework developed by Mattis et al. (2022), large quantities of time series
data are transformed into samples of a small number of QoIs (via PCA-based
feature extraction). Alternatively, the analyst can add parameters; the
original forward map is effectively augmented to make $p\geq q$. One way to do
this is to add a new parameter for each new observable, as suggested in the
stochastic map framework Butler et al. (2020c). In this case, the number of
unknowns is inflated from $p\mapsto p+q>q$.

(III) SIPs accommodate replicates in an awkward fashion. This is very related
to (II) above but can be an issue regardless of $p$’s relation to $q$. Suppose
for example that a computer model has $p$ parameters ($\boldsymbol{\theta}$)
as inputs and $q$ outputs. If there were $n$ vectors of measurements that can
be related to the computer model output, one could not immediately treat the
estimation of $\boldsymbol{\theta}$ as an SIP. Either the data vectors would
have to be combined into a single $q$-vector (to be treated like a population,
I), or following Butler et al. (2020c), $n$ additional $q$-vectors of unknowns
$\boldsymbol{\epsilon}_{i}$ would be added ($p\mapsto p+nq$).

This awkward, forced adjustment contrasts the world of inference where
replicate measurements are a virtual cornerstone of statistical design and
analysis of experiments.

Figure 2:  A demonstration of the fundamental underdetermined nature of the
SIP for $p=2$ and $q=1$. The only way to recover the true
$f_{\boldsymbol{\Theta}}$ distributions (left) from given only $f_{Y}$ (right)
is to use unknown (and unknowable) auxiliary variables whose joint
distribution could be arbitrarily complex.

(IV) The underdeterminedness of an SIPs may preclude the recovery of any
“true” solution. Any solution to (3) solves the SIP. If one is faced with a
problem where there is supposed to be a single, true distribution to be
recovered, then one will almost certainly not be able to recover it for the
following reasons. As we showed in this work, an infinite number of solutions
will exist when $p>q$, as there are infinitely many choices for auxiliary
variables. Any true distribution comes from joint density of hypothetical
observables with the given observables, and this distribution could be very
complicated. Figure 2 shows three potential “true” solutions that come from
complicated joint distributions for augmented observables
$\boldsymbol{Y}_{\text{aug}}$. Is one choice of auxiliary variables and
densities more valid or preferable than others? In general, no.

Even without direct mention of auxiliary variables, there is no _a priori_
reason to prefer BBE, BJW, or any of the possible intuitive solutions. This is
quite different from Bayesian inference where the choice of prior distribution
becomes irrelevant as more data is collected. Moreover, even in the case that
$p=q$, there is the risk of infinitely many solutions if the forward map is
not uniquely invertible on a given domain, as shown in Section A of the
appendix.

(V) In practice, SIP solutions may not exist. Throughout this work we have
taken it as given that the SIP solution exists, i.e., that the range
$\boldsymbol{g}(\mathcal{P})$ contains $\mathcal{Q}$, and that the system of
equations has a solution. Underdetermined systems can, in reality, be
“inconsistent” in that they have none. Therefore, the existence of a solution
needs to be checked and not assumed. Moore (1977), for example, gives a test
to check for solutions to nonlinear equations within given bounds.

###  8.1 Example: CoV Versus Simple Linear Regression

Suppose that two measurements of a response variable are given
$\boldsymbol{y}_{\text{obs}}=(-1,1)^{\top}$ together with an uncertainty
matrix of $\Sigma_{y}=\sigma^{2}\boldsymbol{I}_{2}$. An analyst wants to
relate these to a predictor variable $x$, having values of $(-1,1)$, through a
forward map which is linear in its parameters:
$g_{x}(\boldsymbol{\theta})=\theta_{1}+\theta_{2}x=[1\ x]\boldsymbol{\theta}$.

To treat this as an SIP, the analyst might take $\boldsymbol{Y}\sim
N(\boldsymbol{y}_{\text{obs}},\boldsymbol{\Sigma}_{y})$, or that the given
measurement values were actually population parameters of a Gaussian
distribution (I). The linear forward map is indexed by $x$, but for the
problem at hand reduces to $\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\Theta}$
with $\boldsymbol{X}=\begin{bmatrix}1&-1\\\ 1&1\end{bmatrix}$. The SIP is
solved by
$f^{\text{CoV}}_{\boldsymbol{\Theta}}(\boldsymbol{\theta})=N_{2}\big{(}\boldsymbol{X}^{-1}\boldsymbol{y}_{\text{obs}},\
\boldsymbol{X}^{-1}\boldsymbol{\Sigma}_{y}\boldsymbol{X}^{-\top}\big{)}$ which
simplifies to

$\displaystyle\boldsymbol{\Theta}$ $\displaystyle\sim
N_{2}\left(\begin{bmatrix}0\\\ 1\end{bmatrix},\
\frac{\sigma^{2}}{2}\begin{bmatrix}1&0\\\ 0&1\end{bmatrix}\right)$

When interpolating and extrapolating at some new covariate value $x^{*}$, the
“predictive” distribution is a univariate Gaussian derived from a second CoV
(II):

$[1\ x^{*}]\boldsymbol{\Theta}\sim
N\left(x^{*},\frac{\sigma^{2}}{2}(1+x^{*2})\right)\ .$

The analyst is thus able to use two datapoints to get the distribution on
$\boldsymbol{\Theta}$ and predict (with quantified uncertainty) at any $x^{*}$
without any checks on model assumptions (because there is no falsifiable model
being used!). If more measurements were obtained for new $x$ values (resulting
in $n$ total), the analyst would have to modify their approach to ensure that
the number of parameters was greater. This would be done by 1] reducing the
$n$ measurements to $q=1$ or 2 observables; 2] expanding the
$\boldsymbol{\Theta}$ vector and making the rows of the matrix
$\boldsymbol{X}$ correspond to a higher-order polynomial in $x$, ensuring that
a new forward map will overfit; or 3] augmenting the original forward map to
include $n$ new parameters $\epsilon_{1},\ldots,\epsilon_{n}$, as per Butler
et al. (2018a).

A statistical analyst would model any number of measurements as
$(Y_{i}|\boldsymbol{\theta})=\theta_{1}+\theta_{2}x_{i}+E_{i}$
($i=1,\ldots,n$), _conditional upon_ the covariate values $x_{i}$ and unknown
parameters $\boldsymbol{\theta}$; a marginal, unconditional distribution would
_not_ be specified for $\boldsymbol{Y}$ as in an SIP. The right side of the
model equation above is a random variable (uppercase $Y_{i}$) due to the fact
that the $E_{i}$ term (and only this term) on the left side has a
distribution, say $N(0,\sigma^{2})$. Each given measurement is a _realization_
of a random variable (lowercase $y_{i}$). All these data values will be
present in the likelihood function, but the argument is $\boldsymbol{\theta}$,
not $\boldsymbol{y}$. A Bayesian might specify a flat prior for
$\boldsymbol{\Theta}$ and possibly even take $\sigma^{2}$ as given to complete
the analysis. In this case it is true that if $n=2$, the posterior would be
the same as the CoV solution, but no reasonable Bayesian would be comfortable
with this answer, especially if prediction were required. Additional
measurements beyond $n=2$ pose no philosophical problem for the statistician,
and indeed, the only dimensional consideration on the statistician’s mind is
ensuring $n>>p=2$.

###  8.2 Example: Calibration of Nuclear Reaction Code

Here we give an overview of a real type of analysis that can appear to be
either a stochastic or statistical IP and is thus useful for pedagogical
purposes.

Consider a nuclear reaction model in which a collection of ingoing actinide
isotopes are subjected to a fission and/or fusion environment controlled by
physical fluence parameters and a set of nuclear cross-sections. The
conversion of nuclides within the reactions is governed by a set of
differential equations, and the code output is a vector of the resulting
actinides at the end of the reaction series. The total number of fissions is
also tabulated by the code. In order to compare the outputs to real
measurements, per-fission ratios are formed for each resultant actinide (=
number of atoms / total fissions, so that scaling factors cancel). The goal is
to estimate the physics parameters (and in the case of the true forensics
scenario, ratios of ingoing isotopes) using the observed per-fission
quantities.

Figure 3: Schematic illustrating the forward problem for a simple scenario of
nuclear reactions followed by radioactive decay.

There are multiple actinide per-fission measurements with uncertainties, but
these correspond to radiochemical analysis performed some time after the
nuclear event. To account for this time differential, either the Bateman
equations (describing radioactive decay) must augment the burn code, or the
samples must be decay-corrected back to the end-of-event time ($t=0$) and then
used for the calibration. That is, the forward map for this analysis can
either be thought of as the composition of two maps between physics parameters
and observables
$\boldsymbol{\theta}\stackrel{{\scriptstyle\text{code}}}{{\longrightarrow}}\boldsymbol{Y}_{0}\stackrel{{\scriptstyle\text{decay}}}{{\longrightarrow}}\boldsymbol{Y}_{t}$,
or solely the code. Figure 3 shows this graphically for
$\boldsymbol{Y}^{\top}=({}^{237}$U$,{}^{237}$Np$,{}^{241}$Pu$,{}^{241}$Am).

Inverting the physics of radioactive decay for a measurement vector is
unequivocally a CoV/SIP. It is therefore tempting to consider the entire
calibration process as SIP. It is especially tempting to do so when, as in
many historic radiochemical reports, the original collection of independent
decay-corrected samples $\boldsymbol{y}_{0,1},\ldots,\boldsymbol{y}_{0,n}$ has
been reduced to a single vector with uncertainty.

The most immediate reasons for not treating this nuclear calibration scenario
as an SIP are practical. First, there are often fewer parameters than
observables. This is almost always true when ingoing isotopics are known, and
the goal is to estimate a few fluence parameters. Treating known masses as
unknown to ensure $p>q$ is not appealing. Second, an SIP solution does not
typically exist because the code cannot simultaneously fit all the responses.
After adding discrepancy terms, one is still faced with the possibility of
infinitely many solutions since it is not known where the augmented map is
one-to-one.

The fundamental reason for not treating this IP as an SIP comes from thinking
about what the data represent. The radiochemical measurements are samples from
a population of such quantities, and any reduction to summary values—such as
in the historic reports—does not change this fact. There are unknown
population parameters $\boldsymbol{\theta}$, and any collection of samples
could have been generated by a single $\boldsymbol{\theta}$. Within the SIP
framework, nature’s distribution of $\boldsymbol{\Theta}$ would have to change
between collections of samples in order to produce their respective
variability. Because this is not how the data is generated, an SIP is not
appropriate to infer $\boldsymbol{\theta}$.

##  9 Conclusion

This paper explored so-called “stochastic” inverse problems (SIPs) in great
depth. For the majority of the paper these problems were taken at face value
and various solutions explored. First we provided intuitive solutions derived
from a change-of-variables (CoV) wherein the user explicitly controls $p-q$
degrees of freedom. The two existing types of solutions in the literature (BBE
and BJW) were then shown to be derivable from CoVs. We then showed that any
solution to an SIP must be directly related to a CoV. After not questioning
the SIP framework, we then gave a lengthy discussion wherein we pointed out a
number of fundamental issues inherent to SIPs. This work thus demonstrates
that anyone wanting to treat an estimation or prediction problem as an SIP (as
opposed to one of statistical inference) must answer the following questions:

  * •

Was the data generated in a way consistent with the SIP framework? How will
any possible replicates be treated?

  * •

Does the stipulation $p\geq q$ make sense for this problem?

  * •

Does an SIP solution exist?

  * •

If the answer to the point above is “yes,” then: Given that an infinite number
of SIP solutions are possible, why is one preferable to any other? On the
other hand, given that the SIP solution is almost certainly not the true
distribution, will this be problematic for future tasks, e.g., testing,
filtering, smoothing, and prediction?

##  Appendix A Example: SIP Involving a Two-to-One Map

As a simple demonstration of the method used in the proof of the Proposition
1, consider the function $g(\theta)=\theta^{2}$ and a random variable $Y$ on
$\mathcal{Q}=(0,1)$. If the domain is chosen to be
$\mathcal{P}\stackrel{{\scriptstyle\text{def}}}{{=}}(-1,1)$, then we can take
$\mathcal{P}_{1}\stackrel{{\scriptstyle\text{def}}}{{=}}(-1,0)$ and
$\mathcal{P}_{2}\stackrel{{\scriptstyle\text{def}}}{{=}}(0,1)$. The pullback
of $f_{Y}$ through
$g_{1}^{-1}(y)\stackrel{{\scriptstyle\text{def}}}{{=}}-\sqrt{y}$ and
$g_{2}^{-1}(y)\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{y}$ yields
densities $-2\theta
f_{Y}(\theta^{2})\mathbbm{1}\\{\theta\in\mathcal{P}_{1}\\}$ and $2\theta
f_{Y}(\theta^{2})\mathbbm{1}\\{\theta\in\mathcal{P}_{2}\\}$. Let $0<w<1$ and
consider the continuous mixture of the two pullpack densities:

$\displaystyle f^{\text{CoV}}_{\Theta,w}(\theta)$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}2|\theta|f_{Y}(\theta^{2})\Big{(}w\mathbbm{1}\\{\theta\in\mathcal{P}_{1}\\}+(1-w)\mathbbm{1}\\{\theta\in\mathcal{P}_{2}\\}\Big{)}\
.$

By the CoV Theorem, the distribution of
$\Gamma\stackrel{{\scriptstyle\text{def}}}{{=}}\Theta^{2}$ is

$\displaystyle f_{\Gamma}(\gamma)$
$\displaystyle=f^{\text{CoV}}_{\Theta,w}\big{(}g_{1}^{-1}(\gamma)\big{)}\Big{|}\frac{\partial
g_{1}^{-1}}{\partial\gamma}\Big{|}+f^{\text{CoV}}_{\Theta,w}\big{(}g_{2}^{-1}(\gamma)\big{)}\Big{|}\frac{\partial
g_{2}^{-1}}{\partial\gamma}\Big{|}$
$\displaystyle=\frac{1}{2\sqrt{\gamma}}\left(f^{\text{CoV}}_{\Theta,w}\big{(}-\sqrt{\gamma}\big{)}+f^{\text{CoV}}_{\Theta,w}\big{(}\sqrt{\gamma}\big{)}\right)$
$\displaystyle=f_{Y}(\gamma)\Big{(}w\mathbbm{1}\\{-\sqrt{\gamma}\in\mathcal{P}_{1}\\}+(1-w)\mathbbm{1}\\{-\sqrt{\gamma}\in\mathcal{P}_{2}\\}\Big{)}$
$\displaystyle\ \
+f_{Y}(\gamma)\Big{(}w\mathbbm{1}\\{\sqrt{\gamma}\in\mathcal{P}_{1}\\}+(1-w)\mathbbm{1}\\{\sqrt{\gamma}\in\mathcal{P}_{2}\\}\Big{)}$
$\displaystyle=f_{Y}(\gamma)\Big{(}w\mathbbm{1}\\{-\sqrt{\gamma}\in\mathcal{P}_{1}\\}+(1-w)\mathbbm{1}\\{\sqrt{\gamma}\in\mathcal{P}_{2}\\}\Big{)}\quad=f_{Y}(\gamma)\
.$

To get a glimpse into the more general case, if the domain had been chosen to
be $\mathcal{P}\stackrel{{\scriptstyle\text{def}}}{{=}}(-\epsilon,1)$ for
$0<\epsilon<1$, the domain sets would be
$\mathcal{P}_{1}\stackrel{{\scriptstyle\text{def}}}{{=}}(-\epsilon,0)$ and
$\mathcal{P}_{2}\stackrel{{\scriptstyle\text{def}}}{{=}}(0,\epsilon)$, but now
$\mathcal{P}_{3}$ would be the set $(\epsilon,1)$ where the function is one-
to-one. The infinite set of solutions would take the form

$\displaystyle f^{\text{CoV}}_{\Theta,w}(\theta)$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{\propto}}2|\theta|f_{Y}(\theta^{2})\cdot\left\\{\begin{array}[]{cc}w\mathbbm{1}\\{\theta<0\\}+(1-w)\mathbbm{1}\\{\theta>0\\}&-\epsilon<\theta<\epsilon\\\
1&\phantom{-}\epsilon<\theta<1\end{array}\right.\ .$

Using $Y\sim$ $Unif(0,1)$ with $\epsilon=0.5$, two solutions are shown in 4.

Figure 4:  Two of the infinitely many solutions when $Y\sim$Unif$(0,1)$ and
$g(\theta)=\theta^{2}$ on $(-0.5,1)$ determined by the choice of $w$.

##  Appendix B Additional Derivations

###  B.1 BJW Solution: Linear Transformation of a Multivariate Gaussian Vector

Here we verify the claim in Section 5.1 that the proposed Gaussian
distribution is indeed the BJW solution under the linear map
$\boldsymbol{g}(\boldsymbol{\theta})=\boldsymbol{A}\boldsymbol{\theta}$.

######  Lemma 1.

For $\boldsymbol{\Theta}\sim f^{\text{BJW}}\equiv
N_{p}(\wtilde{\boldsymbol{\mu}},\wtilde{\boldsymbol{\Sigma}})$ with moments
defined by 13, 15, $\boldsymbol{A}\boldsymbol{\Theta}\sim
N_{q}(\boldsymbol{\mu}_{y},\boldsymbol{\Sigma}_{y})$, confirming that the
proposed BJW solution indeed solves the SIP (3).

###### Proof.

We need to show that
$\boldsymbol{A}\wtilde{\boldsymbol{\mu}}=\boldsymbol{\mu}_{y}$ and
$\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{\top}=\boldsymbol{\Sigma}_{y}$.
We will show the second claim about covariances and then use it to show the
first claim about the means.

Starting from 15,

$\wtilde{\boldsymbol{\Sigma}}^{-1}=\boldsymbol{\Sigma}_{\theta}^{-1}+\boldsymbol{A}^{T}\big{\\{}\boldsymbol{\Sigma}_{y}^{-1}-(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T})^{-1}\big{\\}}\boldsymbol{A}\
.$

Let $\boldsymbol{D}=\boldsymbol{\Sigma}_{\theta}^{-1}$,
$\boldsymbol{U}=\boldsymbol{A}^{T}$,
$\boldsymbol{E}=\boldsymbol{\Sigma}_{y}^{-1}-(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T})^{-1}$,
and $\boldsymbol{V}=\boldsymbol{A}$. The Woodbury formula states
$(\boldsymbol{D}+\boldsymbol{U}\boldsymbol{E}\boldsymbol{V})^{-1}=\boldsymbol{D}^{-1}-\boldsymbol{D}^{-1}\boldsymbol{U}(\boldsymbol{E}^{-1}+\boldsymbol{V}\boldsymbol{D}^{-1}\boldsymbol{U})^{-1}\boldsymbol{V}\boldsymbol{D}^{-1}$,
so that

$\displaystyle\wtilde{\boldsymbol{\Sigma}}=\boldsymbol{\Sigma}_{\theta}-\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T}\left\\{(\boldsymbol{\Sigma}_{y}^{-1}-(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T})^{-1})^{-1}+\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T}\right\\}^{-1}\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\
.$

Now consider $\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}$,
and let
$\boldsymbol{B}=\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T}$.
Then,

$\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}=\boldsymbol{B}-\boldsymbol{B}\left\\{(\boldsymbol{\Sigma}_{y}^{-1}-\boldsymbol{B}^{-1})^{-1}+\boldsymbol{B}\right\\}^{-1}\boldsymbol{B}\
.$

Now redefine
$\boldsymbol{D}=(\boldsymbol{\Sigma}_{y}^{-1}-\boldsymbol{B}^{-1})^{-1}$ and
$\boldsymbol{E}=\boldsymbol{B}$ (with
$\boldsymbol{U}=\boldsymbol{V}=\boldsymbol{I}$), and again apply the Woodbury
formula. Then
$(\boldsymbol{D}+\boldsymbol{E})^{-1}=\boldsymbol{B}^{-1}-\boldsymbol{B}^{-1}\boldsymbol{\Sigma}_{y}\boldsymbol{B}^{-1}$.
Thus,

$\displaystyle\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}$
$\displaystyle=\boldsymbol{B}-\boldsymbol{B}(\boldsymbol{B}^{-1}-\boldsymbol{B}^{-1}\boldsymbol{\Sigma}_{y}\boldsymbol{B}^{-1})\boldsymbol{B}$
$\displaystyle=\boldsymbol{\Sigma}_{y}\ .$

Next, to see the first claim about the mean vectors, start from 13 and add and
subtract the term
$\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}\boldsymbol{\mu}_{\theta}$
within the parentheses:

$\displaystyle\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\wtilde{\boldsymbol{\Sigma}}\left(\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{\mu}_{y}-\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}\boldsymbol{\mu}_{\theta}+\big{\\{}\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}\boldsymbol{-}\boldsymbol{A}^{T}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{T})^{-1}\boldsymbol{A}\big{\\}}\boldsymbol{\mu}_{\theta}+\boldsymbol{\Sigma}^{-1}_{\theta}\boldsymbol{\mu}_{\theta}\right)$
$\displaystyle=\wtilde{\boldsymbol{\Sigma}}\left(\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{\mu}_{y}-\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}\boldsymbol{\mu}_{\theta}+\big{\\{}\wtilde{\boldsymbol{\Sigma}}^{-1}-\boldsymbol{\Sigma}^{-1}_{\theta}\big{\\}}\boldsymbol{\mu}_{\theta}+\boldsymbol{\Sigma}^{-1}_{\theta}\boldsymbol{\mu}_{\theta}\right)$
$\displaystyle=\wtilde{\boldsymbol{\Sigma}}\left(\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{\mu}_{y}-\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}\boldsymbol{\mu}_{\theta}\right)+\boldsymbol{\mu}_{\theta}$
$\displaystyle=\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}\boldsymbol{\Sigma}_{y}^{-1}(\boldsymbol{\mu}_{y}-\boldsymbol{A}\boldsymbol{\mu}_{\theta})+\boldsymbol{\mu}_{\theta}\
.$

Now consider $\boldsymbol{A}\wtilde{\boldsymbol{\mu}}$ and use the fact above
that
$\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}=\boldsymbol{\Sigma}_{y}$
so that

$\displaystyle\boldsymbol{A}\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\big{(}\boldsymbol{A}\wtilde{\boldsymbol{\Sigma}}\boldsymbol{A}^{T}\big{)}\boldsymbol{\Sigma}_{y}^{-1}(\boldsymbol{\mu}_{y}-\boldsymbol{A}\boldsymbol{\mu}_{\theta})+\boldsymbol{A}\boldsymbol{\mu}_{\theta}$
$\displaystyle=\boldsymbol{\mu}_{y}\ .$

∎

###  B.2 BJW Solution: Gaussian Mean Estimation as “Stochastic Map” Inversion

The forward map is
$\boldsymbol{g}(\boldsymbol{\theta})\stackrel{{\scriptstyle\text{def}}}{{=}}\boldsymbol{A}\boldsymbol{\theta}$
where
$\boldsymbol{A}\stackrel{{\scriptstyle\text{def}}}{{=}}[\boldsymbol{1}_{n}\vdots\boldsymbol{I}_{n}]$
and
$\boldsymbol{\theta}_{n}\stackrel{{\scriptstyle\text{def}}}{{=}}(\mu,\epsilon_{1},\ldots,\epsilon_{n})^{\top}$.
The observable distribution is
$N_{n}(\boldsymbol{\mu}_{y},\boldsymbol{\Sigma}_{y})$ with
$\boldsymbol{\Sigma}_{y}=\sigma^{2}_{y}\boldsymbol{I}_{n}$. The initial
distribution is $\underaccent{\wtilde}{f}_{\boldsymbol{\Theta}}\equiv
N_{n+1}(\boldsymbol{\mu}_{\theta},\boldsymbol{\Sigma}_{\theta})$ with
$\boldsymbol{\mu}_{\theta}^{\top}=(\mu_{0},\boldsymbol{0}_{n}^{\top})$ and
$\boldsymbol{\Sigma}_{\theta}=\text{diag}(\sigma_{0}^{2},\sigma^{2}_{\epsilon},\ldots,\sigma^{2}_{\epsilon})$.
The BJW solution in this case is
$N_{n+1}(\wtilde{\boldsymbol{\mu}},\wtilde{\boldsymbol{\Sigma}})$ with mean
and covariance given by 13, 15, or equivalently,

$\displaystyle\wtilde{\boldsymbol{\Sigma}}^{-1}\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\boldsymbol{A}^{\top}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{\mu}_{y}-\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})^{-1}\boldsymbol{A}\boldsymbol{\mu}_{\theta}+\boldsymbol{\Sigma}_{\theta}^{-1}\boldsymbol{\mu}_{\theta}$
$\displaystyle\wtilde{\boldsymbol{\Sigma}}^{-1}$
$\displaystyle=\boldsymbol{A}^{\top}\boldsymbol{\Sigma}_{y}^{-1}\boldsymbol{A}-\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})^{-1}\boldsymbol{A}+\boldsymbol{\Sigma}_{\theta}^{-1}\
.$

For any $\sigma^{2}$ term, we will take $\tau=1/\sigma^{2}$ to be the
corresponding precision. Also define
$\boldsymbol{J}_{n}=\boldsymbol{1}_{n}\boldsymbol{1}_{n}^{\top}$. The most
complicated term in the equations above is

$\displaystyle\boldsymbol{A}^{\top}(\boldsymbol{A}\boldsymbol{\Sigma}_{\theta}\boldsymbol{A}^{\top})^{-1}\boldsymbol{A}$
$\displaystyle=\begin{bmatrix}nc_{2}&c_{2}\boldsymbol{1}_{n}^{\top}\\\
c_{2}\boldsymbol{1}_{n}&\tau_{\epsilon}\boldsymbol{I}_{n}-c_{1}\boldsymbol{J}_{n}\end{bmatrix}$
$\displaystyle c_{1}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{\tau^{2}_{\epsilon}}{n\tau_{\epsilon}+\tau_{0}}\quad\quad
c_{2}\stackrel{{\scriptstyle\text{def}}}{{=}}\tau_{\epsilon}-nc_{1}=\frac{\tau_{\epsilon}\tau_{0}}{n\tau_{\epsilon}+\tau_{0}}\
.$

Now the precision matrix is

$\displaystyle\wtilde{\boldsymbol{\Sigma}}^{-1}$
$\displaystyle=\begin{bmatrix}n\tau_{y}&\tau_{y}\boldsymbol{1}_{n}^{\top}\\\
\tau_{y}\boldsymbol{1}_{n}&\tau_{y}\boldsymbol{I}_{n}\end{bmatrix}-\begin{bmatrix}nc_{2}&c_{2}\boldsymbol{1}_{n}^{\top}\\\
c_{2}\boldsymbol{1}_{n}&\tau_{\epsilon}\boldsymbol{I}_{n}-c_{1}\boldsymbol{J}_{n}\end{bmatrix}+\begin{bmatrix}\tau_{0}&\boldsymbol{0}_{n}^{\top}\\\
\boldsymbol{0}_{n}&\tau_{\epsilon}\boldsymbol{I}_{n}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}nc_{3}+\tau_{0}&c_{3}\boldsymbol{1}_{n}^{\top}\\\
c_{3}\boldsymbol{1}_{n}&\tau_{y}\boldsymbol{I}_{n}+c_{1}\boldsymbol{J}_{n}\end{bmatrix}\quad\quad
c_{3}\stackrel{{\scriptstyle\text{def}}}{{=}}\tau_{y}-c_{2}\ .$

The covariance is

$\displaystyle\wtilde{\boldsymbol{\Sigma}}$
$\displaystyle=\begin{bmatrix}a&\boldsymbol{b}^{\top}\\\
\boldsymbol{b}&\boldsymbol{C}\end{bmatrix}^{-1}=\begin{bmatrix}\frac{1}{a}+\frac{1}{a^{2}}\boldsymbol{b}^{\top}(2,2)\boldsymbol{b}&-\frac{1}{a}\boldsymbol{b}^{\top}(2,2)\\\
-\frac{1}{a}(2,2)\boldsymbol{b}&(2,2)\end{bmatrix}$ $\displaystyle(2,2)^{-1}$
$\displaystyle=\boldsymbol{C}-\frac{1}{a}\boldsymbol{b}\boldsymbol{b}^{\top}=\tau_{y}\boldsymbol{I}_{n}+c_{4}\boldsymbol{J}_{n}$
$\displaystyle(2,2)$
$\displaystyle=\sigma^{2}_{y}\boldsymbol{I}_{n}-c_{5}\boldsymbol{J}_{n}$
$\displaystyle c_{4}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{\tau_{\epsilon}^{2}}{n\tau_{\epsilon}+\tau_{0}}-\frac{c_{3}^{2}}{nc_{3}+\tau_{0}}\quad\quad
c_{5}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{\sigma^{4}_{y}c_{4}}{n\sigma^{2}_{y}c_{4}+1}\
.$

and hence

$\displaystyle\wtilde{\boldsymbol{\Sigma}}$
$\displaystyle=\begin{bmatrix}c_{7}&c_{6}\boldsymbol{1}_{n}^{\top}\\\
c_{6}\boldsymbol{1}_{n}&\sigma^{2}_{y}\boldsymbol{I}_{n}-c_{5}\boldsymbol{J}_{n}\end{bmatrix}\quad=\begin{bmatrix}\mathcal{O}(n^{-1})&-\mathcal{O}(n^{-1})\boldsymbol{1}_{n}^{\top}\\\
-\mathcal{O}(n^{-1})\boldsymbol{1}_{n}&\sigma^{2}_{y}\boldsymbol{I}_{n}-\mathcal{O}(n^{-1})\boldsymbol{J}_{n}\end{bmatrix}$
(52) $\displaystyle c_{6}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}-\frac{c_{3}(\sigma^{2}_{y}-nc_{5})}{nc_{3}+\tau_{0}}\quad\quad
c_{7}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{nc_{3}+\tau_{0}}+\frac{nc_{3}^{2}(\sigma^{2}_{y}-nc_{5})}{(nc_{3}+\tau_{0})^{2}}\
.$

Returning to the mean, if
$\overline{\mu}_{y}\stackrel{{\scriptstyle\text{def}}}{{=}}n^{-1}\sum_{i=1}^{n}\mu_{i}$,
then

$\displaystyle\wtilde{\boldsymbol{\Sigma}}^{-1}\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\begin{bmatrix}\tau_{y}\boldsymbol{1}_{n}^{\top}\boldsymbol{\mu}_{y}-nc_{2}\mu_{0}+\tau_{0}\mu_{0}\\\
\tau_{y}\boldsymbol{\mu}_{y}-c_{2}\mu_{0}\boldsymbol{1}_{n}\end{bmatrix}$
$\displaystyle\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\begin{bmatrix}c_{7}&c_{6}\boldsymbol{1}_{n}^{\top}\\\
c_{6}\boldsymbol{1}_{n}&\sigma^{2}_{y}\boldsymbol{I}_{n}-c_{5}\boldsymbol{1}_{n}\boldsymbol{1}_{n}^{\top}\end{bmatrix}\begin{bmatrix}n\tau_{y}\overline{\mu}_{y}+(\tau_{0}-nc_{2})\mu_{0}\\\
\tau_{y}\boldsymbol{\mu}_{y}-c_{2}\mu_{0}\boldsymbol{1}_{n}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}(nc_{6}+nc_{7})\tau_{y}\overline{\mu}_{y}+\left(-nc_{2}c_{6}+c_{7}(\tau_{0}-nc_{2})\right)\mu_{0}\\\
\boldsymbol{\mu}_{y}+(nc_{6}-nc_{5})\tau_{y}\overline{\mu}_{y}\boldsymbol{1}_{n}+\left(-c_{2}\sigma^{2}_{y}+nc_{2}c_{5}+c_{6}(\tau_{0}-nc_{2})\right)\mu_{0}\boldsymbol{1}_{n}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\left(\sigma^{2}_{y}-\mathcal{O}(n^{-1})\right)\tau_{y}\overline{\mu}_{y}+\mathcal{O}(n^{-1})\mu_{0}\\\
\boldsymbol{\mu}_{y}+\left(-\sigma^{2}_{y}+\mathcal{O}(n^{-1})\right)\tau_{y}\overline{\mu}_{y}\boldsymbol{1}_{n}-\mathcal{O}(n^{-1})\mu_{0}\boldsymbol{1}_{n}\end{bmatrix}\
,$

and hence

$\displaystyle\wtilde{\boldsymbol{\mu}}$
$\displaystyle=\begin{bmatrix}\overline{\mu}_{y}+\mathcal{O}(n^{-1})\\\
(\boldsymbol{\mu}_{y}-\overline{\mu}_{y}\boldsymbol{1}_{n})+\mathcal{O}(n^{-1})\boldsymbol{1}_{n}\end{bmatrix}\
.$ (53)

## References

  * Aster et al. (2019) Aster, R. C., Borchers, B., and Thurber, C. H. (2019), Parameter Estimation and Inverse Problems, Amsterdam: Elsevier, 3rd ed.
  * Baptista et al. (2021) Baptista, R., Marzouk, Y., Morrison, R. E., and Zahm, O. (2021), “Learning non-Gaussian graphical models via Hessian scores and triangular transport,” arXiv preprint arXiv:2101.03093.
  * Bernardo and Smith (1994) Bernardo, J. and Smith, A. F. M. (1994), Bayesian Theory, Chichester: John Wiley & Sons.
  * Breidt et al. (2011) Breidt, J., Butler, T., and Estep, D. (2011), “A measure-theoretic computational method for inverse sensitivity problems I: Method and analysis,” SIAM Journal on Numerical Analysis, 49, 1836–1859.
  * Bruder et al. (2020) Bruder, L., Gee, M. W., and Wildey, T. (2020), “Data-consistent solutions to stochastic inverse problems using a probabilistic multi-fidelity method based on conditional densities,” International Journal for Uncertainty Quantification, 10, 399–424.
  * Butler and Estep (2013) Butler, T. and Estep, D. (2013), “A numerical method for solving a stochastic inverse problem for parameters,” Annals of Nuclear Energy, 52, 86–94.
  * Butler et al. (2012) Butler, T., Estep, D., and Sandelin, J. (2012), “A computational measure theoretic approach to inverse sensitivity problems II: A posteriori error analysis,” SIAM Journal on Numerical Analysis, 50, 22–45.
  * Butler et al. (2014) Butler, T., Estep, D., Tavener, S., Dawson, C., and Westerink, J. J. (2014), “A measure-theoretic computational method for inverse sensitivity problems III: Multiple quantities of interest,” SIAM/ASA Journal on Uncertainty Quantification, 2, 174–202.
  * Butler et al. (2015a) Butler, T., Graham, L., Estep, D., Dawson, C., and Westerink, J. J. (2015a), “Definition and solution of a stochastic inverse problem for the Manning’s $n$ parameter field in hydrodynamic models,” Advances in Water Resources, 78, 60–79.
  * Butler et al. (2017) Butler, T., Graham, L., Mattis, S., and Walsh, S. (2017), “A measure-theoretic interpretation of sample based numerical integration with applications to inverse and prediction problems under uncertainty,” SIAM Journal on Scientific Computing, 39, A2072–2098.
  * Butler and Hakula (2020) Butler, T. and Hakula, H. (2020), “What do we hear from a drum? A data-consistent approach to quantifying irreducible uncertainty on model inputs by extracting information from correlated model output data,” Computer Methods in Applied Mechanics and Engineering, 370, 113228.
  * Butler et al. (2015b) Butler, T., Huhtala, A., and Juntunen, M. (2015b), “Quantifying uncertainty in material damage from vibrational data,” Journal of Computational Physics, 283, 414–435.
  * Butler et al. (2018a) Butler, T., Jakeman, J., and Wildey, T. (2018a), “Combining push-forward measures and Bayes’ Rule to construct consistent solutions to stochastic inverse problems,” SIAM Journal on Scientific Computing, 40, A984–A1011.
  * Butler et al. (2018b) — (2018b), “Convergence of probability densities using approximate models for forward and inverse problems in uncertainty quantification,” SIAM Journal on Scientific Computing, 40, A3523–A3548.
  * Butler et al. (2020a) — (2020a), “Optimal experimental design for prediction based on push-forward probability measures,” Journal of Computational Physics, 416, 109518.
  * Butler et al. (2020b) Butler, T., Pilosov, M., and Walsh, S. (2020b), “Simulation-based optimal experimental design: A measure-theoretic perspective,” Tech. rep., University of Colorado.
  * Butler et al. (2020c) Butler, T., Wildey, T., and Yen, T. Y. (2020c), “Data-consistent inversion for stochastic input-to-output maps,” Inverse Problems, 36, 085015.
  * Casella and Berger (2002) Casella, G. and Berger, R. L. (2002), Statistical Inference, Thomas Learning: Duxbury, 2nd ed.
  * Chiachio-Ruano et al. (2022) Chiachio-Ruano, J., Chiachio-Ruano, M., and Sankararaman, S. (eds.) (2022), Bayesian Inverse Problems: Fundamentals and Engineering Applications, Boca Raton: CRC Press.
  * Davison (2003) Davison, A. C. (2003), Statistical Models, Cambridge University Press.
  * Gelman et al. (2014) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2014), Bayesian Data Analysis, Boca Raton: CRC Press, 3rd ed.
  * Graham et al. (2017) Graham, L., Butler, T., Walsh, S., Dawson, C., and Westerink, J. J. (2017), “A measure-theoretic algorithm for estimating bottom friction in a coastal inlet: Case study of Bay St. Louis during Hurricane Gustav (2008),” Monthly Weather Review, 145, 929–954.
  * Groetsch (1993) Groetsch, C. W. (1993), Inverse Problems in the Mathematical Sciences, Springer Fachmedien Wiesbaden.
  * Hannig (2009) Hannig, J. (2009), “On generalized fiducial inference,” Statistica Sinica, 19, 491–544.
  * Hannig et al. (2016) Hannig, J., Iyer, H., Lai, R. C. S., and Lee, T. C. M. (2016), “Generalized fiducial inference: A review and new results,” Journal of the American Statistical Association, 111, 1346–1361.
  * Kaipio and Somersalo (2005) Kaipio, J. P. and Somersalo, E. (2005), Statistical and Computational Inverse Problems, New York: Springer.
  * Kirsch (2011) Kirsch, A. (2011), An Introduction to the Mathematical Theory of Inverse Problems, New York: Springer, 3rd ed.
  * Lee (2013) Lee, J. M. (2013), Introduction to Smooth Manifolds, New York: Springer, 2nd ed.
  * Lesnic (2021) Lesnic, D. (2021), Inverse Problems with Applications in Science and Engineering, Boca Raton: CRC Press.
  * Marzouk et al. (2016) Marzouk, Y., Moselhy, T., Parno, M., and Spantini, A. (2016), “Sampling via measure transport: An introduction,” Handbook of Uncertainty Quantification, 1–41.
  * Mattis and Butler (2019) Mattis, S. A. and Butler, T. (2019), “Enhancing piecewise-defined surrogate response surfaces with adjoints on sets of unstructured samples to solve stochastic inverse problems,” International Journal for Numerical Methods in Engineering, 119, 923–940.
  * Mattis et al. (2015) Mattis, S. A., Butler, T. D., Dawson, C. N., Estep, D., and Vesselinov, V. V. (2015), “Parameter estimation and prediction for groundwater contamination based on measure theory,” Water Resources Research, 51, 7608–7629.
  * Mattis et al. (2022) Mattis, S. A., Steffen, K. R., Butler, T., Dawson, C. N., and Estep, D. (2022), “Learning quantities of interest from dynamical systems for observation-consistent inversion,” Computer Methods in Applied Mechanics and Engineering, 388, 114230.
  * Moore (1977) Moore, R. E. (1977), “A test for existence of solutions to nonlinear systems,” SIAM Journal on Numerical Analysis, 14, 611–615.
  * Pawitan (2001) Pawitan, Y. (2001), In All Likelihood: Statistical Modelling and Inference Using Likelihood, Oxford University Press.
  * Presho et al. (2017) Presho, M., Mattis, S., and Dawson, C. (2017), “Uncertainty quantification of two-phase flow problems via measure theory and the generalized multiscale finite element method,” Computational Geosciences, 21, 187–204.
  * Reid (2010) Reid, N. (2010), “Likelihood inference,” Wiley Interdisciplinary Reviews: Computational Statistics, 2, 517–525.
  * Robert (2007) Robert, C. P. (2007), The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation, New York: Springer, 2nd ed.
  * Shao (2003) Shao, J. (2003), Mathematical Statistics, New York: Springer, 2nd ed.
  * Stuart (2010) Stuart, A. M. (2010), “Inverse problems: A Bayesian perspective,” Acta Numerica, 19, 451–559.
  * Swigon et al. (2019) Swigon, D., Stanhope, S. R., Zenker, S., and Rubin, J. E. (2019), “On the importance of the Jacobian determinant in parameter inference for random parameter and random measurement error models,” SIAM/ASA Journal on Uncertainty Quantification, 7, 975–1006.
  * Tarantola (2005) Tarantola, A. (2005), Inverse Problem Theory and Methods for Model Parameter Estimation, Philadelphia: SIAM.
  * Tikhonov et al. (1995) Tikhonov, A. N., Goncharsky, A. V., Stepanov, V. V., and Yagola, A. G. (1995), Numerical Methods for the Solution of Ill-Posed Problems, Dordrecht: Springer Science+Business Media.
  * Tran and Wildey (2021) Tran, A. and Wildey, T. (2021), “Solving stochastic inverse problems for property-structure linkages using data-consistent inversion and machine learning,” JOM: The Journal of The Minerals, Metals & Materials Society, 73, 72–89.
  * Uy and Grigoriu (2019) Uy, W. I. T. and Grigoriu, M. (2019), “Specification of additional information for solving stochastic inverse problems,” SIAM Journal on Scientific Computing, 41, A2880–A2910.
  * Walsh et al. (2018) Walsh, S. N., Wildey, T. M., and Jakeman, J. D. (2018), “Optimal experimental design using a consistent Bayesian approach,” ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, 4, 1–19.