# Generative AI for Bayesian Computation

Nicholas G. Polson111Nick Polson is a Professor in Econometrics and Statistics
at the U Chicago Booth, email<EMAIL_ADDRESS>We thank the participants
of the GDRR 13 meeting in Madrid, May 24-26
Booth School of Business
University of Chicago
Vadim Sokolov222Vadim Sokolov is an Associate Professor in Operations Research
at George Mason University, email<EMAIL_ADDRESS>
Department of Systems Engineering and Operations Research
George Mason University

(First Draft: Dec 10, 2022
This Draft: February 12, 2024)

###### Abstract

Bayesian Generative AI (BayesGen-AI) methods are developed and applied to
Bayesian computation. BayesGen-AI reconstructs the posterior distribution by
directly modeling the parameter of interest as a mapping (a.k.a. deep learner)
from a large simulated dataset. This provides a generator that we can evaluate
at the observed data and provide draws from the posterior distribution. This
method applies to all forms of Bayesian inference including parametric models,
likelihood-free models, prediction and maximum expected utility problems.
Bayesian computation is then equivalent to high dimensional non-parametric
regression. Bayes Gen-AI main advantage is that it is density-free and
therefore provides an alternative to Markov Chain Monte Carlo. It has a number
of advantages over vanilla generative adversarial networks (GAN) and
approximate Bayesian computation (ABC) methods due to the fact that the
generator is simpler to learn than a GAN architecture and is more flexible
than kernel smoothing implicit in ABC methods. Design of the Network
Architecture requires careful selection of features (a.k.a. dimension
reduction) and nonlinear architecture for inference. As a generic
architecture, we propose a deep quantile neural network and a uniform base
distribution at which to evaluate the generator. To illustrate our
methodology, we provide two real data examples, the first in traffic flow
prediction and the second in building a surrogate for satellite drag data-set.
Finally, we conclude with directions for future research.

##  1 Introduction

Our goal is to develop Bayesian Generative AI (BayesGen-AI) algorithms. Bayes
Gen-AI calculates the posterior distributions and functionals of interest that
arise in machine learning task and statistical inference by directly writing
the parameter of interest as a stochastic mapping (a.k.a deep learner) and
thus avoids the use of densities and Markov Chain Monte Carlo (MCMC)
simulation, by directly learning the inverse posterior mapping from parameters
$\theta$ to data $y$, our method doesn’t rely on the use of iterative
simulation methods or objective functions that involve density calculations.
BayesGen-AI applies to all forms of Bayesian inference including parametric
models, likelihood-free models, prediction and maximum expected utility
problems. To illustrate our methodology, we use two real world data examples:
first, a nonlinear traffic flow prediction problem based on Polson and
Sokolov, (2017, 2015) and, second, a satellite drag dataset.

BayesGen-AI starts with a large sample from a joint distribution of
observables, parameters and an independent base distribution, denoted by
$(y^{(i)},\theta^{(i)},\tau^{(i)})_{i=1}^{N}$ where $N$ is large. For the data
generating process we allow for both parameter likelihood inference
$y^{(i)}|\theta^{(i)}\sim p(y|\theta^{(i)})$ and likelihood-free models where
$y^{(i)}=f(\theta^{(i)})$ is a given forward map. Samples are generated from
the prior $\theta^{(i)}\sim\pi(\theta)$, and a base distribution
$\tau^{(i)}\sim p(\tau)$, typically a vector of uniforms. We also allow for
the addition of a latent variable which is stochastic $z^{(i)}\sim p(z)$ or
deterministic regressor variables where $z^{(i)}=x$ for some given $x$. Hence,
our methodology also applies to conditional density estimation.

When the base distribution is uniform, the inverse Bayes map is simply given
by the multivariate inverse posterior cdf where we write
$\theta=F^{-1}_{\theta\mid y}(\tau)$, where $\tau\sim p(\tau)$. Hence, given a
training data set, we can directly find the mapping
$\theta^{(i)}=H\left(y^{(i)},\tau^{(i)}\right)$ by learning $H$ as a deep
neural network from a large simulated training dataset. The deep learner
requires the specification of a sequence of layers including feature
extraction, dimension reduction and nonlinearity to pattern match the training
data. Designing the NN architecture for $H$ is the main challenge in Bayes
Gen-AI. For a generic choice, we propose the use of quantile neural networks
(QNNs). Our work builds on Jiang et al., (2017) who were the first to propose
deep learners for dimension reduction methods and to provide asymptotic
theoretical results. We also build on the insight by Dabney et al., 2018a ;
Ostrovski et al., (2018) that implicit quantile neural networks can be used to
approximated posterior distributions that arise in decision theory. Dabney et
al., (2017) also show the connection between the widely used $1$-Wasserstein
distance and quantile regression.

Given an observable data point $y$ and a new base draw $\tau$, the posterior
$p(\theta|y)$ is simply determined from the map $\theta=H(y,\tau)$ where
$\tau\sim U(0,1)$. From this we can calculate any quantiles of densities that
we need. From this we also show how to calculate any posterior functionals
$E(f(\theta)|y)$ by integration the appropriate quantile function using the
trapezoidal rule. We provide an $O(N^{-4})$ approximation bound showing the
advantage of quantile methods in estimating functionals. If the base measure
is chosen to be uniform, deep Quantile NNs provide a general framework for
training $H$ and for Bayesian decision-making. Our method also provides an
alternative to invertible NN approaches such as normalizing flows, to Gaussian
Process surrogates, and we show how ABC methods can be viewed within our
framework.

Architecture design to learn an inverse CDF (quantile function), namely
$F^{-1}(\tau,y)=f_{\theta}(\tau,y)$, will use a kernel embedding by augmenting
the input space to the network. The quantile function is represented as a
superposition of two other functions
$F^{-1}(\tau,y)=f_{\theta}(\tau,y)=g(\psi(y)\circ\phi(\tau))$ where $\circ$ is
the element-wise multiplication operator. Both functions $g$ and $\psi$ are
feed-forward neural networks and $\phi$ is a cosine embedding transform. To
avoid over-fitting, we use a sufficiently large training dataset, see Dabney
et al., 2018a  in a reinforcement learning context. Dimension reduction
(a.k.a. feature extraction) will draw on methods used in ABC, see Jiang et
al., (2017) for an approach using deep learning.

Generative AI also allow the researcher to learn a dimensionality-reduced
summary (a.k.a. sufficient) statistics along with a non-linear map (Jiang et
al.,, 2017; Albert et al.,, 2022). A useful interpretation of the sufficient
statistic as a posterior mean, which also allows us to view posterior inputs
as one of the inputs to the posterior mean. One can also view a NN is an
approximate nonlinear Bayes filter to perform such tasks (Müller et al.,,
1997). Our framework provides a natural link for black box methods and
stochastic methods, as commonly known in the machine learning literature
(Bhadra et al.,, 2021; Breiman,, 2001). Quantile deep neural networks and
their ReLu/$\tanh$ counterparts provide a natural architecture. Approximation
properties of those networks are discussed in Polson and Ročková, (2018).
Dimensionality reduction can be performed using auto-encoders and partial
least-squares (Polson et al.,, 2021) due to the result by Brillinger, (2012);
Bhadra et al., (2021), see survey by Blum et al., (2013) and kernel embeddings
approach discussed by Park et al., (2016). Generative AI circumvents the need
for methods like MCMC that require the density evaluations.

ABC methods can also be viewed as a generative approach. Here we show that $H$
can be viewed as a nearest neighbor network. ABC relies on comparing summary
statistic $S$ calculated from the observed data to the summary statistics
calculated from the simulated data. Jiang et al., (2017) show that a natural
choice of $S$ is via the posterior mean. Papamakarios and Murray, (2016) shows
how to use mixture density networks (Bishop,, 1994) to approximate the
posterior for ABC calculations. In an ABC framework with a parametric
exponential family, see Beaumont et al., (2002) and Nunes and Balding, (2010)
for the optimal choice of summary statistics. A local smoothing version of ABC
is given in Jiang et al., (2018); Bernton et al., (2019); Fearnhead and
Prangle, (2012), Longstaff and Schwartz, (2001) take a basis function
approach. Pastorello et al., (2003) provide an estimation procedure when
latent variables are present.

In the statistical and engineering literature, generative models often arise
in the context of inverse problems Baker et al., (2022) and decision-making
Dabney et al., (2017). In the context of inverse problems, the prediction of a
mean is rarely is an option, since average of several correct values is not
necessarily a correct value and might not even be feasible from the physical
point of view. The two main approaches are surrogate-based modeling and
approximate Bayes computations (ABC) Park et al., (2016); Blum et al., (2013);
Beaumont et al., (2002). Surrogate-based modeling Gramacy, (2020) is a general
approach to solve inverse problems, which relies on the availability of a
forward model $y=f(\theta)$, which is a deterministic function of parameters
$\theta$. The forward model is used to generate a large sample of pairs
$(y,\theta)$, which is then used to train a surrogate, typically a Gaussian
Process, which can be used to calculate the inverse map. For a recent review
of the surrogate-based approach see Baker et al., (2022). There are multiple
papers that address different aspects of surrogate-based modeling.

Bayes Gen-AI provides an alternative to non-parametric Gaussian process-based
surrogates which heavily rely on the informational contribution of each sample
point and quickly becomes ineffective when faced with significant increases in
dimensionality (Shan and Wang,, 2010; Donoho,, 2000). Further, the homogeneous
GPs models predict poorly (Binois et al.,, 2018). Unfortunately, the
consideration of each input location to handle these heteroskedastic cases
result in analytically intractable predictive density and marginal likelihoods
(Lázaro-Gredilla et al.,, 2010). Furthermore, the smoothness assumption made
by GP models hinders capturing rapid changes and discontinuities in the input-
output relations. Popular attempts to overcome these issues include relying on
the selection of kernel functions using prior knowledge about the target
process (Cortes et al.,, 2004); splitting the input space into subregions so
that inside each of those smaller subregions the target function is smooth
enough and can be approximated with a GP model (Gramacy and Lee,, 2008;
Gramacy and Apley,, 2015; Chang et al.,, 2014); and learning spatial basis
functions (Bayarri et al.,, 2007; Wilson et al.,, 2014; Higdon,, 2002).

For low-dimensional $\theta$, the simplest approach is to discretize the
parameter space and the data space and use a lookup table to approximate
$p(\theta\mid y)$. This is the approach taken by Jiang et al., (2017).
However, this approach is not scalable to high-dimensional $\theta$. For
practical cases, when the dimension of $\theta$ is high, we can use
conditional independence structure present in the data to decompose the joint
distribution into a product of lower-dimensional functions (Papamakarios et
al.,, 2017). In machine learning literature a few approaches were proposed
that rely on such a decomposition (van den Oord and Kalchbrenner,, 2016;
Germain et al.,, 2015; Papamakarios et al.,, 2017). Most of those approaches
use KL divergence as a metric of closeness between the target distribution and
the learned distribution.

A natural approach to model the posterior distribution using a neural network
is to assume that parameters of a neural network are random variables to use
either approximation or MCMC techniques to model the posterior distribution
over the parameters (Izmailov et al.,, 2021). A slightly different approach is
to assume that only weights of the last output layer of a neural network are
stochastic (Wang et al.,, 2022; Schultz et al.,, 2022). BayesGen-AI provides a
natural alternative to these methods.

The paper is outlined as follows. Section 1 provides a review of the existing
literature. Sections 2 describes our GenAI-Bayes Bayesian algorithm. Section 3
describes quantile neural network architectures. Sections 4 provides two real
world applications from traffic flow prediction problem and a satellite drag
data-set. Finally, Section 5 concludes with directions for future research.

##  2 Generative AI for Bayesian Computation

In this section we provide a general description of the Bayes Gen-AI
algorithm. Then we describe architecture design and quantile neural networks.
We begin by establishing the following notations that will be used.

$\displaystyle y=$ $\displaystyle\;\;{\rm outcome\;of\;interest}$
$\displaystyle\theta=$ $\displaystyle\;\;{\rm parameters}$ $\displaystyle z=$
$\displaystyle\;\;{\rm latent\;variables}$ $\displaystyle\tau=$
$\displaystyle\;\;{\rm base\;distribution}$

Let $(y,\theta)\in\Re^{n+k}$ be observable data and parameters. Write
$\theta=(\theta_{1},\ldots,\theta_{k})\in\Re^{k}$. Let $\tau\in\Re^{K}$ be a
base measure and $z\in\Re^{D}$ be a vector of latent variables. Suppose
$\theta\sim\pi(\theta)$, $z\sim p(z)$ and $\tau\sim p(\tau)$. A natural choice
will be $\tau$ a vector of uniforms and $K=k$. In the presence of latent
variables, the likelihood is given by $p(y|\theta)=\int p(y,z|\theta)dz$.

We allow for the special case where the latent variables are observed
regressor variable and $p(z)=\delta_{\\{z=x\\}}$ where $\delta$ denotes a
Krönecker delta function. In this case, we can find the posterior
$p(\theta|x,y)$. Moreover, we allow for likelihood-free inference where
$y=f(\theta)$ for some deterministic map. The key property is that we will
never have to know the densities to learn the inverse Bayes map.

The goal is to compute the posterior distribution $p(\theta\mid y)$. The
underlying assumptions are that $\theta\sim\pi(\theta)$ a prior distribution.
Our framework allows for many forms of stochastic data generating processes.
The dynamics of the data generating process are such that it is
straightforward to simulate from a so-called forward model or traditional
stochastic model, namely

$y=f(\theta)\;\;{\rm or}\;\;y|\theta\sim p(y|\theta).$ (1)

The idea is quite straightforward, if we could perform high dimensional non-
parametric regression, we could simulate a large training dataset of
observable parameter, data and base pairs, denoted by
$(y^{(i)},\theta^{(i)},\tau^{(i)})_{i=1}^{N}$. If latent variables are part of
the model specification we simply attach $z^{(i)}$ to the triple. The inverse
Bayes map is then given by

$\theta\stackrel{{\scriptstyle D}}{{=}}H(y,\tau),$ (2)

where $\tau$ is the vector with elements from the baseline distribution, such
as Gaussian, are simulated training data. Here $\tau$ is a vector of standard
uniforms. The function $H:\Re^{k}\times\Re^{D}\rightarrow\Re^{k}$ is a deep
neural network. The function $H$ is again trained using the simulated data
$(y^{(i)},\theta^{(i)})_{i=1}^{N}$, via quantile regression using the
simulated data

$\theta^{(i)}=H(y^{(i)},\tau^{(i)}),~{}i=1,\ldots,N.$

The main construction is the specification of the network structure implicit
in $H$. We will provide guidelines in the following sections.

In many practical applications $y$ is high-dimensional, and we can improve the
performance of the deep neural network by using a summary statistic $S(y)$.
Here $S:\Re^{N}\rightarrow\Re^{k}$ is a $k$-dimensional sufficient statistic,
which is a lower-dimensional representation of $y$. The summary statistic is a
sufficient statistic in the Bayes sense (Kolmogorov,, 1942). The main idea is
to use a deep neural network to approximate the posterior mean
$E_{\pi}(\theta|y)$ as the optimal estimate of $S(y)$.

Assuming that we have fitted the deep neural network, $H$, to the training
data, we can use the estimated inverse map to evaluate at new $y$ and $\tau$
to obtain a set of posterior samples for any new $y$ by interpolation and
evaluating the map

$\theta\stackrel{{\scriptstyle D}}{{=}}H_{N}(S(y),\tau),$ (3)

where $H_{N}$ denotes the estimated map.

The caveat being is to how to choose $H$ and how well the deep neural network
interpolates for the new inputs. We also have flexibility in choosing the
distribution of $\tau$, for example, we can also for $\tau$ to be a high-
dimensional vector of Gaussians, and essentially provide a mixture-Gaussian
approximation for the set of posterior. Gen-AI in a simple way is using
pattern matching to provide a look-up table for the map from $y$ to $\theta$.
Bayesian computation has then being replaced by the optimisation performed by
Stochastic Gradient Descent (SGD). Thus, we can sequentially update $H$ using
SGD step as new data arrives. MCMC, in comparison, is computationally
expensive and needs to be re-run for any new data point. In our examples, we
discuss choices of architectures for $H$ and $S$. Specifically, we propose
cosine-embedding for transforming $\tau$.

#### Bayes Gen-AI Algorithm:

A necessary condition is the ability to simulate from the parameters, latent
variables, and data process. This generates a (potentially large) triple

$\left\\{y^{(i)},\theta^{(i)},\tau^{(i)}\right\\}_{i=1}^{N},$

By construction, the posterior distribution can be characterized by the von
Neumann inverse CDF map. For $\theta=(\theta_{1},\ldots,\theta_{k})$ we have

$\theta\stackrel{{\scriptstyle D}}{{=}}F_{\theta|y}^{-1}(\tau),\;\;{\rm
where}\;\;\tau\sim U(0,1)$

is a vector of standard uniforms.

Given the observed data $y$ and new base draw $\tau$, we then simply evaluate
the following posterior map

$\theta\stackrel{{\scriptstyle D}}{{=}}H(S(y),\tau)$

When the base draw $\tau$ is uniform of same dimension as the parameter
vector, the map $H$ is precisely the inverse posterior cdf. This characterises
the posterior distribution $p(\theta|y)$ of interest. As simulation and
fitting deep neural networks is computationally cheap, we typically choose $N$
to be of order $10^{6}$ or more.

In many cases, we can estimate a summary statistic, $S$ also via deep
learning, and we write the inverse cdf as a composite map, and we have to fit

$\theta^{(i)}=H(S(y^{(i)}),\tau^{(i)})\;\;{\rm
where}\;\;F_{\theta|y}^{-1}=H\circ S$

There are many variations in the architecture design and construction. It is
useful to interpret $S$ as a summary statistic providing dimension reduction
in the $y$ space.

There is flexibility in the choice of base distribution. For example, we can
replace the $\tau$ with a different distribution that we can easily sample
from. One example is a multivariate Normal, proposed for diffusion models
Sohl-Dickstein et al., (2015). The dimensionality of the normal can be large.
The main insight is that you can solve a high-dimensional least squares
problem with nonlinearity using stochastic gradient descent. Deep quantile NNs
provide a natural candidate of deep learners. Dabney et al., 2018b  show that
implicit quantile networks can be used to for distributional reinforcement
learning and has proposed the use of quantile regression as a method for
minimizing the 1-Wasserstein in the univariate case when approximating using a
mixture of Dirac functions. Figure 1 illustrates the different in using a
quantile objective versus pure density simulation based on latent variable
approach.

The 1-Wasserstein distance is the $\ell_{1}$ metric on the inverse
distribution function. It is also known as earth mover’s distance and can be
calculated using order statistics (Levina and Bickel,, 2001). For quantile
functions $F^{-1}_{U}$ and $F^{-1}_{V}$ the 1-Wasserstein distance is given by

$W_{1}(F^{-1}_{U},F^{-1}_{V})=\int_{0}^{1}|F^{-1}_{U}(\tau)-F^{-1}_{V}(\tau)|d\tau.$

Weng, (2019) shows Wasserstein GAN networks out-perform vanilla GAN due to the
improved quantile metric $q=F^{-1}_{U}(\tau)$ minimize the expected quantile
loss

$E_{U}[\rho_{\tau}(u-q)]$

Quantile regression can be shown to minimize the 1-Wasserstein metric. A
related loss is the quantile divergence,

$q(U,V)=\int_{0}^{1}\int_{F^{-1}_{U}(q)}^{F^{-1}_{V}(q)}(F_{U}(\tau)-q)dqd\tau.$

Figure 1: Deep Quantile Neural Network. Source: Ostrovski et al., (2018)

Hence, even if new $y$ is not in the simulated training dataset we can still
learn the posterior map of interest. The Kolmogorov-Arnold theorem says any
multivariate function can be expressed this way. So in principle if sample
size is large enough we can learn the manifold structure in the parameters for
any arbitrary nonlinearity. As the dimension of the data $y$ is large, in
practice, this requires providing an efficient architecture.

The following provides a summary of Bayes Gen-AI algorithms:

Algorithm 1 Gen-AI for Bayesian Computation (GenAI-Bayes)

Simulate $\theta^{(i)}\sim\pi(\theta)$. Simulate
$(y^{(i)}\mid\theta^{(i)})_{1\leq i\leq N}\sim p(y\mid\theta)$ or
$y^{(i)}=f(\theta^{(i)})$ and $\tau^{(i)}\sim p(\tau)$

Train $H$ using the simulated dataset for $i=1,\ldots N$, via
$\theta^{(i)}=H(y^{(i)},\tau^{(i)})$

For the observed $y$, calculate a sample from $p(\theta\mid y)$ with a new
base draw $\tau\sim p(\tau)$ by  $\theta\stackrel{{\scriptstyle
D}}{{=}}H(y,\tau)$

In the case when $y$ is high-dimensional, and we can learn a summary statistic
$S(y)$, we can use the following algorithm:

Algorithm 2 GenAI-Bayes with Dimension Reduction

Simulate $\theta^{(i)}\sim\pi(\theta)$. Simulate
$(y^{(i)}\mid\theta^{(i)})_{1\leq i\leq N}\sim p(y\mid\theta)$ or
$y^{(i)}=f(\theta^{(i)})$ and $\tau^{(i)}\sim p(\tau)$

Learn a summary statistic $S(y)$ using a deep neural network

Train $H$ using the simulated dataset for $i=1,\ldots N$, via
$\theta^{(i)}=H(S(y^{(i)}),\tau^{(i)})$

For the observed $y$, calculate a sample from $p(\theta\mid y)$ with a new
base draw $\tau\sim p(\tau)$ by  $\theta\stackrel{{\scriptstyle
D}}{{=}}H(S(y),\tau)$

Consider now the case where there exist latent variables, $z$, related to the
data density via

$p(y|\theta)=\int p(y,z|\theta)dz$

Again assume that and when $z|\theta$ and $y|z,\theta$ are straightforward to
simulate either from a density of forward model. Hence, we can simulate
$y|\theta$ by composition.

Algorithm 3 Gen-AI with Latent Variables (GenAI-Bayes-Latent)

Simulate $\theta^{(i)}\sim\pi(\theta)$. Simulate
$y^{(i)},z^{(i)}\mid\theta^{(i)}\sim p(y,z\mid\theta)$, $i=1,\ldots,N$ or
$y^{(i)}=f(\theta^{(i)})$.

Let $S$ denotes feature variables: Train $H$ and $S$, using
$\theta^{(i)}=H(S(z^{(i)},y^{(i)}),\tau^{(i)})$, where $\tau^{(i)}\sim
p(\tau)$

For a given $y$, calculate a sample from $p(\theta\mid y)$ with $\tau$ a new
base draw  $\theta\stackrel{{\scriptstyle D}}{{=}}H\left(S((z,y),\tau\right)$

To learn the feature selection variables $S(y)$, Jiang et al., (2017) propose
to use a deep learner to approximate the posterior mean $E_{\pi}(\theta|y)$ as
the optimal estimate of $S(y)$. Specifically, they construct a minimum squared
error estimator $\hat{\theta}(y)$ from a large simulated dataset and calculate

$\min_{\psi}\frac{1}{N}\sum_{i=1}^{N}\|f_{\psi}(y^{(i)})-\theta^{(i)}\|^{2}_{2}$

where $f_{\psi}$ denotes a DNN with parameter $\psi$. The resulting estimator
$\hat{\theta}(y)=f_{\hat{\psi}}(y^{(i)})$ approximates the feature vector
$S(y)$.

Together with the following architecture for the summary statistic neural
network

$\displaystyle H^{(1)}=$
$\displaystyle\tanh\left(W^{(0)}H^{(0)}+b^{(0)}\right)$ $\displaystyle
H^{(2)}=$ $\displaystyle\tanh\left(W^{(1)}H^{(1)}+b^{(1)}\right)$
$\displaystyle\vdots$ $\displaystyle H^{(L)}=$
$\displaystyle\tanh\left(W^{(L-1)}H^{(L-1)}+b^{(L-1)}\right)$
$\displaystyle\hat{\theta}=$ $\displaystyle W^{(L)}H^{(L)}+b^{(L)},$

where $H^{(0)}=\theta$ is the input, and $\hat{\theta}$ is the summary
statistic output. ReLU activation function can be used instead of $\tanh$. In
our applications we use ReLU.

There is still the question of approximation and the interpolation properties
of a DNN.

#### Folklore of Deep Learning:

Shallow Deep Learners provide good representations of multivariate functions
and are good interpolators.

A number of authors have provided theoretical approximation results that apply
to BayesGen-AI. For example, interpolation properties of quantile neural
networks were recently studied by Padilla et al., (2022) and Shen et al.,
(2021), Schmidt-Hieber, (2020)

###  2.1 Dimension Reduction

Given $y$, the posterior density is denoted by $p(\theta\mid y)$. Here
$y=(y_{1},\ldots,y_{n})$ is high dimensional. Moreover, we need the set of
posterior probabilities $\pi_{\theta\mid y}(\theta\in B\mid y)$ for all Borel
sets $B$. Hence, we need two things, dimension reduction for $y$. The whole
idea is to estimate ”maps” (a.k.a. transformations/feature extraction) of the
output data $y$, thus, it is reduced to uniformity.

There is a nice connection between the posterior mean and the sufficient
statistics, especially minimal sufficient statistics in the exponential
family. If there exists a sufficient statistic $S^{*}$ for $\theta$, then
Kolmogorov, (1942) shows that for almost every $y$, $\pi(\theta\mid
y)=\pi(\theta\mid S^{*}(y))$, and further $S(y)=E_{\pi}(\theta\mid
y)=E_{\pi}(\theta\mid S^{*}(y))$ is a function of $S^{*}(y)$. In the special
case of an exponential family with minimal sufficient statistic $S^{*}$ and
parameter $\theta$, the posterior mean $S(y)=E_{\pi}(\theta\mid y)$ is a one-
to-one function of $S^{*}(y)$, and thus is a minimal sufficient statistic.

Deep learners are good interpolators (one of the folklore theorems of machine
learning). Hence, the set of posteriors $p(\theta\mid y)$ is characterized by
the distributional identity

$\theta\stackrel{{\scriptstyle D}}{{=}}H(S(y),\tau_{k}),\;\;{\rm
where}\;\;\tau_{k}\sim U(0,1),$

where $y=(y_{1},\ldots,y_{n})$ and $\tau_{k}$ is a vector of standard
uniforms.

#### Summary Statistic:

Let $S(y)$ is sufficient summary statistic in the Bayes sense (Kolmogorov,,
1942), if for every prior $\pi$

$f_{B}(y):=\pi_{\theta\mid y}(\theta\in B\mid y)=\pi_{\theta\mid
s(y)}(\theta\in B\mid s(y)).$

Then we need to use our pattern matching dataset $(y^{(i)},\theta^{(i)})$
which is simulated from the prior and forward model to ”train” the set of
functions $f_{B}(y)$, where we pick the sets $B=(-\infty,q]$ for a quantile
$q$. Hence, we can then interpolate the mapping to the observed data.

The notion of a summary statistic is prevalent in the ABC literature and is
tightly related to the notion of a Bayesian sufficient statistic $S^{*}$ for
$\theta$, then (Kolmogorov 1942), for almost every $y$,

$\pi(\theta\mid Y=y)=\pi(\theta\mid S^{*}(Y)=S^{*}(y))$

Furthermore, $S(y)=\mathrm{E}\left(\theta\mid
Y=y\right)=\mathrm{E}_{\pi}\left(\theta\mid S^{*}(Y)=S^{*}(y)\right)$ is a
function of $S^{*}(y)$. In the case of exponential family, we have
$S(Y)=\mathrm{E}_{\pi}\left(\theta|Y\right)$ is a one-to-one function of
$S^{*}(Y)$, and thus is a minimal sufficient statistic.

Sufficient statistics are generally kept for parametric exponential families,
where $S(\cdot)$ is given by the specification of the probabilistic model.
However, many forward models have an implicit likelihood and no such
structures. The generalization of sufficiency is a summary statistics (a.k.a.
feature extraction/selection in a neural network). Hence, we make the
assumption that there exists a set of features such that the dimensionality of
the problem is reduced

Learning $S$ can be achieved in a number of ways. First, $S$ is of fixed
dimension $S\in\Re^{s}$ even though $y=(y_{1},\ldots y_{n})$. Typical
architectures include Auto-encoders and traditional dimension reduction
methods. Polson et al., (2021) propose to use a theoretical result of
Brillinger methods to perform a linear mapping $S(y)=Wy$ and learn $W$ using
PLS. Nareklishvili et al., (2022) extend this to IV regression and casual
inference problems.

#### PLS

Another architecture for finding summary statistics is to use PLS. Given the
parameters and data, the map is

$\theta^{(i)}=H\left(By^{(i)},\epsilon\right)$

We can find a set of linear maps $S(y)=By$. This rule also provide dimension
reduction. Moreover, due to orthogonality of $y,\epsilon$, we can simply
consistently estimate $\hat{S}(y)=\hat{B}y$ via $\theta^{(i)}=H(\hat{B}y)$. A
key result of Brillinger, (2012) shows that we can use linear SGD methods and
partial least squares to find $\hat{B}$.

#### ABC

Approximate Bayesian Computations (ABC) is a common approach in cases when
likelihood is not available, but samples can be generated from some model,
e.g. epidemiological simulator. The ABC rely on comparing summary statistic
calculated from data $s_{y}$ and of from observed output $s(\theta)$ and
approximates

$p(\theta,s(\theta)\mid
s_{y},\epsilon)\propto\pi(\theta)f(s(\theta)\mid\theta)K_{\epsilon}\left(\|s(\theta)-s_{y}\|\right).$

Then the approximation to the posterior is simply $p(\theta\mid s_{y})=\int
p_{\epsilon}(\theta,s(\theta)\mid s_{y})ds$. Then the ABC algorithm simply
samples from $p(\theta)$, then generates summary statistic $s(\theta)$ and
rejects the sample with probability proportional to
$K_{\epsilon}(\|s(\theta)-s_{y}\|)$. The Kernel function $K$ can be a simple
indicator function

$K_{\epsilon}(\|s(\theta)-s_{y}\|)=\mathds{1}(\|s(\theta)-s_{y}\|<\epsilon)$

The use of deep neural networks to select $s(\theta)$ has been proposed by
Jiang et al., (2017).

Conventional ABC methods suffers from the main drawback that the samples do
not come from the true posterior, but an approximate one, based on the
$\epsilon$-ball approximation of the likelihood, which is a non-parametric
local smoother. Theoretically, as $\epsilon$ goes to zero, you can guarantee
samples from the true posterior. However, the number of sample required is
prohibitive. Our method circumvents this by replacing the ball with a deep
learning generator and directly model the relationship between the posterior
and a baseline uniform Gaussian distribution. Our method is also not a density
approximation, as many authors have proposed. Rather, we directly use $L_{2}$
methods and Stochastic Gradient Descent to find transport map from $\theta$ to
a uniform or a Gaussian random variable. The equivalent to the mixture of
Gaussian approximation is to assume that our baseline distribution is high-
dimensional Gaussian. Such models are called the diffusion models in
literature. Full Bayesian computations can then be reduced to high-dimensional
$L_{2}$ optimization problems with a carefully chosen neural network.

#### Invertible NNs

In some cases, we can model $H:\Re^{n}\rightarrow\Re^{n}$ as an invertible map
where $x=g(z)$ is generated from a base density $p(z)$. When $g$ is chosen to
be an invertible neural network with a structured diagonal Jacobian that is
easy to compute we can calculate densities. This makes the mapping $g$ far
easier to learn. Our objective function can then be standard log-likelihood
based on

$p(x)=p(z)|{\rm det}J_{z}|^{-1},\;\;{\rm where}\;\;J_{z}=\frac{\partial
g(z)}{\partial z}$

Bayes Flow provides an example of this type of architecture.

The second class of methods proposed on machine leaning literature involves
using deep learners to approximate an inverse CDF function or a more general
approach that represents the target distribution over $\theta$ as a
marginalization over a nuance random variable $z$ (Kingma and Welling,, 2022).
In the case of inverse CDF, the latent variable $z$ is simply uniform on
$(0,1)$ (Bond-Taylor et al.,, 2022). One of the approaches of this type is
called Normalizing flows. Normalizing flows provide an alternative approach of
defining a deterministic map $\theta\mid x,\theta_{g}=G(z,x,\theta_{g})$ that
transforms a univariate random variable $z\sim p(z)$ to a sample from the
target distribution $G(z,x,\theta_{g})=x\sim F(\theta)$. If transformation $G$
is invertible ($G^{-1}$ exists) and differentiable, then the relation between
target density $F$ and the latent density $p(z)$ is given by Rezende and
Mohamed, (2015):

$F(\theta)=p(z)\left|\det\frac{\partial G^{-1}}{\partial
z}\right|=p(z)\left|\det\frac{\partial G}{\partial z}\right|^{-1},$ (4)

where $z=G^{-1}(y)$. A typical procedure for estimating the parameters of the
map $G$ relies on maximizing the log-likelihood

$\log p(z)+\log\left|\det\frac{\partial G^{-1}}{\partial z}\right|$ (5)

The normalizing flow model requires constructing map $G$ that have tractable
inverse and Jacobian determinant. It is achieved by representing $G$ as a
composite map

$G=T_{k}\circ\cdots\circ T_{1},$ (6)

and to use simple building block transformations $T_{i}$ that have tractable
inverse and Jacobian determinant.

We now describe our Quantile Neural Network approach.

##  3 Quantile Neural Networks

###  3.1 Learning Quantiles

The quantile regression likelihood function as an asymmetric function that
penalizes overestimation errors with weight $\tau$ and underestimation errors
with weight $1-\tau$. For a given input-output pair $(x,y)$, and the quantile
function $f(x,\theta)$, parametrized by $\theta$, the quantile loss is
$\rho_{\tau}(u)=u(\tau-I(u<0))$, where $u=y-f(x)$. From the implementation
point of view, a more convenient form of this function is

$\rho_{\tau}(u)=\max(u\tau,u(\tau-1)).$

Given a training data $\\{x_{i},y_{i}\\}_{i=1}^{N}$, and given quantile
$\tau$, the loss is

$L_{\tau}(\theta)=\sum_{i=1}^{N}\rho_{\tau}(y_{i}-f(\tau,x_{i},\theta)).$

Further, we empirically found that adding a means-squared loss to this
objective function, improves the predictive power of the model, thus the loss
function, we use is

$\alpha
L_{\tau}(\theta)+\frac{1}{N}\sum_{i=1}^{N}(y_{i}-f(x_{i},\theta))^{2}.$

One approach to learn the quantile function is to use a set of quantiles
$0<\tau_{1}<\tau_{2},\ldots,\tau_{K}<1$ and then learn $K$ quantile functions
simultaneously by minimizing

$L(\theta,\tau_{1},\ldots,\tau_{K})=\frac{1}{NK}\sum_{i=1}^{N}\sum_{k=1}^{K}\rho_{\tau_{k}}(y_{i}-f_{\tau_{k}}(x_{i},\theta_{k})).$

The corresponding optimization problem of minimizing $L(\theta)$ can be
augmented by adding a non-crossing constraint

$f_{\tau_{i}}(x,\theta_{i})<f_{\tau_{j}}(x,\theta_{j}),~{}\forall X,~{}i<j.$

The non-crossing constraint has been considered by several authors, including
Chernozhukov et al., (2010); Cannon, (2018).

#### Cosine Embedding for $\tau$

To learn an inverse CDF (quantile function)
$F^{-1}(\tau,y)=f_{\theta}(\tau,y)$ we will use a kernel embedding trick and
augment the predictor space. We then represent the quantile function is a
function of superposition for two other functions

$F^{-1}(\tau,y)=f_{\theta}(\tau,y)=g(\psi(y)\circ\phi(\tau))$

where $\circ$ is the element-wise multiplication operator. Both functions $g$
and $\psi$ are feed-forward neural networks. $\phi$ is a cosine embedding. To
avoid over-fitting, we use a sufficiently large training dataset, see Dabney
et al., 2018a  in a reinforcement learning context.

Let $g$ and $\psi$ be feed-forward neural networks and $\phi$ a cosine
embedding given by

$\phi_{j}(\tau)=\mathrm{ReLU}\left(\sum_{i=0}^{n-1}\cos(\pi
i\tau)w_{ij}+b_{j}\right).$

We now illustrate our approach with a simple synthetic dataset.

###  3.2 Synthetic Data

Consider a synthetic data generated from the model

$\displaystyle x$ $\displaystyle\sim U(-1,1)$ $\displaystyle y$
$\displaystyle\sim N(\sin(\pi x)/(\pi x),\exp(1-x)/10).$

The true quantile function is given by

$f_{\tau}(x)=\sin(\pi x)/(\pi x)+\Phi^{-1}(\tau)\sqrt{\exp(1-x)/10},$

|
---|---
(a) Implicit Quantile Network | (b) Explicit Quantile Network
Figure 2: We trained both implicit and explicit networks on the synthetic data
set. The explicit network was trained for three fixed quantiles
(0.05,0.5,0.95). We see no empirical difference between the two

We train two quantile networks, one implicit and one explicit. The explicit
network is trained for three fixed quantiles (0.05,0.5,0.95). Figure 2 shows
fits by both of the networks, we see no empirical difference between the two.

###  3.3 Quantiles as Deep Learners

One can show that quantile models are direct alternatives to other Bayes
computations. Specifically, given $F(y)$, a non-decreasing and continuous from
right function. We define

$Q_{\theta|y}(u):=F^{-1}_{\theta|y}(u)=\inf\left(y:F_{\theta|y}(y)\geq
u\right)$

which is non-decreasing, continuous from left. Now let $g(y)$ to be a non-
decreasing and continuous from left with

$g^{-1}(z)=\sup\left(y:g(y)\leq z\right)$

Then, the transformed quantile has a compositional nature, namely

$Q_{g(Y)}(u)=g(Q(u))$

Hence, quantiles act as superposition (a.k.a. deep Learner).

This is best illustrated in the Bayes learning model. We have the following
result updating prior to posterior quantiles known as the conditional quantile
representation

$Q_{\theta|Y=y}(u)=Q_{Y}(s)\;\;{\rm where}\;\;s=Q_{F(\theta)|Y=y}(u)$

To compute $s$ we use

$u=F_{F(\theta)|Y=y}(s)=P(F(\theta)\leq s|Y=y)=P(\theta\leq
Q_{\theta}(s)|Y=y)=F_{\theta|Y=y}(Q_{\theta}(s))$

Parzen, (2004) also shows the following probabilistic property of quantiles

$\theta=Q_{\theta}(F_{\theta}(\theta))$

Hence, we can increase the efficiency by ordering the samples of $\theta$ and
the baseline distribution as the mapping being the inverse CDF is monotonic.

###  3.4 Normal-Normal Learning

Consider the normal learning model $y|\theta\sim N(\theta,\sigma^{2})\;\;{\rm
and}\;\;\theta\sim N(\mu,\tau^{2})$. Given observations
$(y_{1},\ldots,y_{n})$, we can reduce the likelihood using sufficiency to give

$\bar{y}|\theta,\sigma^{2}\sim
N\left(\theta,\frac{\sigma^{2}}{n}\right)\;\;{\rm and}\;\;\theta\sim
N(\mu,\tau^{2})$

The posterior mean is

$\hat{\theta}(y)=E(\theta|y)=\frac{\sigma^{2}/n}{\sigma^{2}/n+\tau^{2}}\mu+\frac{\tau^{2}}{\sigma^{2}/n+\tau^{2}}\bar{y}$

The prior to posterior quantile map is given by

$1-F_{\theta|y}(\theta)=g(1-F_{\theta}(\theta))$

where $g:[0,1]\rightarrow[0,1]$ is a concentration (aka distortion) function
given by

$g(x)=\Phi\left(w\Phi^{-1}(x)+b\right)\;\;,\;w=\tau\sqrt{\frac{1}{\sigma^{2}/n+\tau^{2}}},b=\frac{1}{\sigma/\sqrt{n}}\sqrt{\frac{\tau^{2}}{\sigma^{2}/n+\tau^{2}}}$

where $\Phi$ is standard normal CDF function.

###  3.5 Calculating Functionals of Interest

Another advantage of quantile methods is that we can find tight error bounds
of $O(N^{-4})$ for estimating functionals of interest. Another important
property of quantiles is that they can be used to calculate functionals of
interest

$E_{\theta|y}(\theta)=\int_{-\infty}^{\infty}\theta\pi(\theta|y)d\theta=\int_{0}^{1}Q_{\theta|y}(\tau)d\tau$

We can then use the trapezoidal rule and obtain an efficient approximation due
lemma of Yakowitz et al., (1978).

######  Proposition 1.

Assume $Q(x)$ is a function with a continuous second derivative on $[0,1]$ and

$\theta\equiv\int_{0}^{1}Q(x)dx.$

If $Y_{0}\equiv 0,Y_{n+1}\equiv 1$, and $\\{Y_{i}\\}_{i=1}^{n}$ is the ordered
sample associated with n independent uniform observations on [0, 1], then for
the estimator

$\theta_{n}=\frac{1}{2}\left[\sum_{i=0}^{n}(Q(Y_{i})+Q(Y_{i+1}))(Y_{i+1}-Y_{i})\right]$

We have that for some constant M,

$E[(\theta_{n}-\theta)^{2}]\leq M/n^{4},\qquad\text{for all }n\geq 1.$

###### Proof.

For any particular sample $Y_{1},Y_{2},...,Y_{n}$, we have

$\theta-\theta_{n}=\sum_{j=1}^{n}\left(\int_{Y_{j}}^{Y_{j+1}}Q(t)dt-\frac{1}{2}[(Q(Y_{j})+Q(Y_{j+1}))(Y_{j+1}-Y_{j})]\right).$

A well-know error bound for the trapezoidal rule is that for any numbers a, b
such that $0\leq a<b\leq 1$,

$\int_{a}^{b}Q(t)dt-\frac{b-1}{2}[Q(b)+Q(a)]=-\frac{(b-a)^{3}}{12}Q^{\prime\prime}(\xi),$

where $a\leq\xi\leq b$. Using this result we have that if
$C\geq|Q^{\prime\prime}(x)|,$ $0\leq x\leq 1$, then

$\displaystyle|\theta-\theta_{n}|$ $\displaystyle\leq$
$\displaystyle\left|\int_{0}^{Y_{1}}Q(t)dt-\frac{1}{2}Y_{1}(Q(0)+Q(1))\right|$
$\displaystyle+\left|\left[\sum_{i=1}^{n-1}\int_{Y_{i}}^{Y_{i+1}}Q(t)dt-\frac{1}{2}[(Q(Y_{i})+Q(Y_{i+1}))(Y_{i+1}-Y_{i})\right]\right|$
$\displaystyle+\left|\int_{Y_{n}}^{1}Q(t)dt-\frac{1}{2}[(Q(Y_{n})+Q(Y_{1}))(1-Y_{n})\right|$
$\displaystyle\leq$ $\displaystyle\frac{C}{12}\sum_{i=1}^{n+1}Z_{i}^{3},$

where the $Z_{i}$’s are as defined as follows.

Let $\\{Y_{i}\\}_{i=1}^{n}$ be the ordered sample constructed from $n$
independent, uniform observations on the unit interval. Define

$Z_{1}=Y_{1},\qquad Z_{j}=Y_{j}-Y_{j-1},\qquad 2\leq j\leq n,\text{ and
}Z_{n+1}=1-Y_{n}$

Yakowitz et al., (1978) provides the following result. Let
$\\{Y_{i}\\}_{i=1}^{n}$ be the ordered sample constructed from $n$
independent, uniform observations on the unit interval, then

$E[Z_{i}^{6}]=6!n!/(n+6)!,\qquad i=1,2,...,n+1,$
$E[Z_{i}^{3}Z_{j}^{3}]=(3!)^{2}n!/(n+1)!,\qquad i,j=1,2,...,n+1,\qquad i\neq
j.$

Consequently, letting $i$ and $j$ range from $1$ to $n+1$, we have

$E[(\theta_{n}-\theta)^{2}]\leq
E\left[\left(\frac{C}{12}\sum_{i}Z_{i}^{3}\right)^{2}\right]=\frac{C^{2}}{144}\left[\sum_{i}E[Z_{i}^{6}]+\sum_{i\neq
j}E[Z_{i}^{3}Z_{j}^{3}]\right]=O\left(n^{-4}\right).$

∎

##  4 Applications

###  4.1 Traffic Data

We further illustrate our methodology, using data from a sensor on interstate
highway I-55. The sensor is located eight miles from the Chicago downtown on
I-55 north bound (near Cicero Ave), which is part of a route used by many
morning commuters to travel from southwest suburbs to the city. As shown on
Figure 3, the sensor is located 840 meters downstream of an off-ramp and 970
meters upstream from an on-ramp.

Figure 3: Implicit neural network for traffic speed observed on I-55 north-
bound towards Chicago

In a typical day traffic flow pattern on Chicago’s I-55 highway, where sudden
breakdowns are followed by a recovery to free flow regime. We can see a
breakdown in traffic flow speed during the morning peak period followed by
speed recovery. The free flow regimes are usually of little interest to
traffic managers. We also, see that variance is low during the free flow
regime and high during the breakdown and recovery regimes.

We use the following architecture to model the implicit quantile function

$\displaystyle\tau_{1}=$
$\displaystyle\mathrm{ReLU}\left(w^{(1)}_{0}+\sum_{i=1}^{64}w^{(1)}_{i}\cos(i\pi\tau)\right)$
$\displaystyle x_{1}=$
$\displaystyle\mathrm{ReLU}\left(w^{(2)}_{0}+\sum_{i=1}^{64}w^{(2)}_{i}x_{i}\right)$
$\displaystyle z=$
$\displaystyle\mathrm{ReLU}\left(w^{(3)}_{0}+\sum_{i=1}^{64}w^{(3)}_{i}\tau_{i}x_{i}\right)$
$\displaystyle y_{1}=$
$\displaystyle\mathrm{ReLU}\left(w^{(4)}_{0}+\sum_{i=1}^{32}w^{(4)}_{i}x_{i}\right)$
$\displaystyle\hat{y}=$ $\displaystyle w^{(5)}_{0}+w^{(5)}_{i}x_{i}.$

We now consider a satellite drag example.

###  4.2 Satellite Drag

Accurate estimation of satellite drag coefficients in low Earth orbit (LEO) is
vital for various purposes such as precise positioning (e.g., to plan
maneuvering and determine visibility) and collision avoidance.

Recently, 38 out of 49 Starlink satellites launched by SpaceX on Feb 3, 2022,
experienced an early atmospheric re-entry caused by unexpectedly elevated
atmospheric drag, an estimated $100MM loss in assets. The launch of the SpaceX
Starlink satellites coincided with a geomagnetic storm, which heightened the
density of Earth’s ionosphere. This, in turn, led to an elevated drag
coefficient for the satellites, ultimately causing the majority of the cluster
to re-enter the atmosphere and burn up. This recent accident shows the
importance of accurate estimation of drag coefficients in commercial and
scientific applications (Berger et al.,, 2023).

Accurate determination of drag coefficients is crucial for assessing and
maintaining orbital dynamics by accounting for the drag force. Atmospheric
drag is the primary source of uncertainty for objects in LEO. This uncertainty
arises partially due to inadequate modeling of the interaction between the
satellite and the atmosphere. Drag is influenced by various factors, including
geometry, orientation, ambient and surface temperatures, and atmospheric
chemical composition, all of which are dependent on the satellite’s position
(latitude, longitude, and altitude).

Los Alamos National Laboratory developed the Test Particle Monte Carlo
simulator to predict the movement of satellites in low earth orbit (Mehta et
al.,, 2014). The simulator takes two inputs, the geometry of the satellite,
given by the mesh approximation and seven parameters, which we list in Table 1
below. The simulator takes about a minute to run one scenario, and we use a
dataset of one million scenarios for the Hubble Space Telescope (Gramacy,,
2020). The simulator outputs estimates of the drag coefficient based on these
inputs, while considering uncertainties associated with atmospheric and gas-
surface interaction models (GSI).

Parameter | Range
---|---
velocity [m/s] | [5500, 9500]
surface temperature [K] | [100, 500]
atmospheric temperature [K] | [200, 2000]
yaw [radians] | $[-\pi,\pi]$
pitch [radians] | $[-\pi/2,\pi/2]$
normal energy AC [unitless] | [0,1]
tangential momentum AC [unitless] | [0,1]
Table 1: Input parameters for the satellite drag simulator.

We use the data set of 1 million simulation runs provided by Sun et al.,
(2019). The data set has 1 million observations, and we use 20% for training
and 80% for testing out-of-sample performance. The model architecture is given
below. We use the Adam optimizer and a batch size of 2048, and train the model
for 200 epochs.

Sauer et al., (2023) provides a survey of modern Gaussian Process based models
for prediction and uncertainty quantification tasks. They compare five
different models, and apply them to the same Hubble data set we use in this
section. We use two metrics to assess the quality of the model, namely RMSE,
which captures predictive accuracy, and continuous rank probability score
(CRPS; Gneiting and Raftery, (2007); Zamo and Naveau, (2018)). Essentially
CRPS is the absolute difference between the predicted and observed cumulative
distribution function (CDF). We use the degenerative distribution with the
entire mass on the observed value (Dirac delta) as to get the observed CDF.
The lower CRPSis, the better.

Their best performing model is treed-GP has the RMSE of 0.08 and CRPS of 0.04,
the worst performing model is the deep GP with approximate “doubly stochastic”
variational inference has RMSE of 0.23 and CRPS of 0.16. The best performing
model in our experiments is the quantile neural network with RMSE of 0.098 and
CRPS of 0.05, which is comparable to the top performer from the survey.

Figure 4 plots of the out-of-sample predictions for forty randomly selected
responses (green crosses) and compares those to 50th quantile predictions
(orange line) and 95% credible prediction intervals (grey region).

Figure 4: Randomly selected 40 out-of-sample observations from the satellite
drag dataset. Green crosses are observed values, orange line is predicted 50%
quantile and the grey region is the 95% credible prediction interval.

Figure 5(a) compares the out-of-sample predictions (50% quantiles) $\hat{y}$
and observed drag coefficients $y$. We can see that histogram resembles a
normal distribution centered at zero, with some “heaviness” on the left tail,
meaning that for some observations, our model under-estimates. The scatterplot
in Figure 5(b) shows that the model is more accurate for smaller values of $y$
and less accurate for larger values of $y$ and values of $y$ at around three.

|
---|---
(a) Histogram of errors $\hat{y}-y$ | (b) $y$ vs $\hat{y}$
Figure 5: Comparison of out-of-sample predictions (50% quantiles) $\hat{y}$
and observed drag coefficients $y$.

Finally, we show histograms of the posterior predictive distribution for four
randomly chosen out-of-sample response values in Figure 6. We can see that the
model concentrates the distribution of around the true values of the response.

|
---|---
(a) Observation 948127 | (b) Observation 722309
|
(a) Observation 608936 | (b) Observation 988391
Figure 6: Posterior histograms for four randomly chosen out-of-sample response
values. The vertical red line is the observed value of the response.

Overall, our model provides a competitive performance to the state-of-the art
techniques used for predicting and UQ analysis of complex models, such as
satellite drag. The model is able to capture the distribution of the response
and provide accurate predictions. The model is also able to provide
uncertainty quantification in the form of credible prediction intervals.

##  5 Discussion

Generative AI is a simulation-based methodology, that takes joint samples of
observables and parameters as an input and then applies nonparametric
regression in a form of deep neural network by regressing $\theta$ on a non-
linear function $h$ which is a function of dimensionality-reduced sufficient
statistics of $\theta$ and a randomly generated stochastically uniform error.
In its simplest form, $h$, can be identified with its inverse CDF.

One solution to the multi-variate case is to use auto-regressive quantile
neural networks. There are also many alternatives to the architecture design
that we propose here. For example, autoencoders Albert et al., (2022); Akesson
et al., (2021) or implicit models, see Diggle and Gratton, (1984); Baker et
al., (2022); Schultz et al., (2022) There is also a link with indirect
inference methods developed in Pastorello et al., (2003); Stroud et al.,
(2003); Drovandi et al., (2011, 2015)

There are many challenging future problems. The method can easily handle high-
dimensional latent variables. But designing the architecture for fixed high-
dimensional parameters can be challenging. We leave it for the future
research. Having learned the nonlinear map, when faced with the observed data
$y_{\mathrm{obs}}$, one simply evaluates the nonlinear map at newly generated
uniform random values. Generative AI circumvents the need for methods like
MCMC that require the density evaluations.

We also think that over-parameterisation of the problem might lead to a useful
model, although it might also lead to non-identifiabaility of the weights in
the regression. Two interesting approaches include, the case when $K>k$ and
mixture models with mixtures of Gaussian for $\tau$.

## References

  * Akesson et al., (2021) Akesson, M., Singh, P., Wrede, F., and Hellander, A. (2021).  Convolutional neural networks as summary statistics for approximate bayesian computation.  IEEE/ACM Transactions on Computational Biology and Bioinformatics.
  * Albert et al., (2022) Albert, C., Ulzega, S., Ozdemir, F., Perez-Cruz, F., and Mira, A. (2022).  Learning Summary Statistics for Bayesian Inference with Autoencoders.
  * Baker et al., (2022) Baker, E., Barbillon, P., Fadikar, A., Gramacy, R. B., Herbei, R., Higdon, D., Huang, J., Johnson, L. R., Ma, P., Mondal, A., et al. (2022).  Analyzing stochastic computer models: A review with opportunities.  Statistical Science, 37(1):64–89.
  * Bayarri et al., (2007) Bayarri, M. J., Berger, J. O., Paulo, R., Sacks, J., Cafeo, J. A., Cavendish, J., Lin, C.-H., and Tu, J. (2007).  A Framework for Validation of Computer Models.  Technometrics, 49(2):138–154.
  * Beaumont et al., (2002) Beaumont, M. A., Zhang, W., and Balding, D. J. (2002).  Approximate Bayesian computation in population genetics.  Genetics, 162(4):2025–2035.
  * Berger et al., (2023) Berger, T. E., Dominique, M., Lucas, G., Pilinski, M., Ray, V., Sewell, R., Sutton, E. K., Thayer, J. P., and Thiemann, E. (2023).  The thermosphere is a drag: The 2022 starlink incident and the threat of geomagnetic storms to low earth orbit space operations.  Space Weather, 21(3).
  * Bernton et al., (2019) Bernton, E., Jacob, P. E., Gerber, M., and Robert, C. P. (2019).  Approximate Bayesian computation with the Wasserstein distance.  Journal of the Royal Statistical Society: Series B, 81(2):235–269.
  * Bhadra et al., (2021) Bhadra, A., Datta, J., Polson, N., Sokolov, V., and Xu, J. (2021).  Merging two cultures: Deep and statistical learning.  arXiv preprint arXiv:2110.11561.
  * Binois et al., (2018) Binois, M., Gramacy, R. B., and Ludkovski, M. (2018).  Practical heteroskedastic Gaussian process modeling for large simulation experiments.  Journal of Computational and Graphical Statistics, 27(4):808–821.
  * Bishop, (1994) Bishop, CM. (1994).  Mixture density networks.  Technical Report NCRG/94/004, Aston University.
  * Blum et al., (2013) Blum, M. G. B., Nunes, M. A., Prangle, D., and Sisson, S. A. (2013).  A Comparative Review of Dimension Reduction Methods in Approximate Bayesian Computation.  Statistical Science, 28(2):189–208.
  * Bond-Taylor et al., (2022) Bond-Taylor, S., Leach, A., Long, Y., and Willcocks, C. G. (2022).  Deep Generative Modelling: A Comparative Review of VAEs, GANs, Normalizing Flows, Energy-Based and Autoregressive Models.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(11):7327–7347.
  * Breiman, (2001) Breiman, L. (2001).  Statistical modeling: The two cultures (with comments and a rejoinder by the author).  Statistical science, 16(3):199–231.
  * Brillinger, (2012) Brillinger, D. R. (2012).  A Generalized Linear Model With “Gaussian” Regressor Variables.  In Guttorp, P. and Brillinger, D., editors, Selected Works of David Brillinger, Selected Works in Probability and Statistics, pages 589–606. Springer, New York, NY.
  * Cannon, (2018) Cannon, A. J. (2018).  Non-crossing nonlinear regression quantiles by monotone composite quantile regression neural network, with application to rainfall extremes.  Stochastic Environmental Research and Risk Assessment, 32(11):3207–3225.
  * Chang et al., (2014) Chang, W., Haran, M., Olson, R., and Keller, K. (2014).  Fast dimension-reduced climate model calibration and the effect of data aggregation.  The Annals of Applied Statistics, 8(2):649–673.
  * Chernozhukov et al., (2010) Chernozhukov, V., Fernández-Val, I., and Galichon, A. (2010).  Quantile and Probability Curves Without Crossing.  Econometrica, 78(3):1093–1125.
  * Cortes et al., (2004) Cortes, C., Haffner, P., and Mohri, M. (2004).  Rational kernels: Theory and algorithms.  Journal of Machine Learning Research, 5(Aug):1035–1062.
  * (19) Dabney, W., Ostrovski, G., Silver, D., and Munos, R. (2018a).  Implicit Quantile Networks for Distributional Reinforcement Learning.
  * (20) Dabney, W., Ostrovski, G., Silver, D., and Munos, R. (2018b).  Implicit quantile networks for distributional reinforcement learning.  In International conference on machine learning, pages 1096–1105. PMLR.
  * Dabney et al., (2017) Dabney, W., Rowland, M., Bellemare, M. G., and Munos, R. (2017).  Distributional Reinforcement Learning with Quantile Regression.
  * Diggle and Gratton, (1984) Diggle, P. J. and Gratton, R. J. (1984).  Monte Carlo Methods of Inference for Implicit Statistical Models.  Journal of the Royal Statistical Society. Series B (Methodological), 46(2):193–227.
  * Donoho, (2000) Donoho, D. L. (2000).  High-dimensional data analysis: The curses and blessings of dimensionality.  In Ams Conference on Math Challenges of the 21st Century.
  * Drovandi et al., (2011) Drovandi, C. C., Pettitt, A. N., and Faddy, M. J. (2011).  Approximate Bayesian computation using indirect inference.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 60(3):317–337.
  * Drovandi et al., (2015) Drovandi, C. C., Pettitt, A. N., and Lee, A. (2015).  Bayesian Indirect Inference Using a Parametric Auxiliary Model.  Statistical Science, 30(1):72–95.
  * Fearnhead and Prangle, (2012) Fearnhead, P. and Prangle, D. (2012).  Constructing summary statistics for approximate Bayesian computation: Semi-automatic approximate Bayesian computation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):419–474.
  * Germain et al., (2015) Germain, M., Gregor, K., Murray, I., and Larochelle, H. (2015).  MADE: Masked Autoencoder for Distribution Estimation.  In Proceedings of the 32nd International Conference on Machine Learning, pages 881–889. PMLR.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction, and estimation.  Journal of the American statistical Association, 102(477):359–378.
  * Gramacy, (2020) Gramacy, R. B. (2020).  Surrogates: Gaussian process modeling, design, and optimization for the applied sciences.  CRC press.
  * Gramacy and Apley, (2015) Gramacy, R. B. and Apley, D. W. (2015).  Local Gaussian Process Approximation for Large Computer Experiments.  Journal of Computational and Graphical Statistics, 24(2):561–578.
  * Gramacy and Lee, (2008) Gramacy, R. B. and Lee, H. K. H. (2008).  Bayesian Treed Gaussian Process Models With an Application to Computer Modeling.  Journal of the American Statistical Association, 103(483):1119–1130.
  * Higdon, (2002) Higdon, D. (2002).  Space and space-time modeling using process convolutions.  In Quantitative methods for current environmental issues, pages 37–56. Springer.
  * Izmailov et al., (2021) Izmailov, P., Vikram, S., Hoffman, M. D., and Wilson, A. G. G. (2021).  What are Bayesian neural network posteriors really like?  In International conference on machine learning, pages 4629–4640. PMLR.
  * Jiang et al., (2017) Jiang, B., Wu, T.-Y., Zheng, C., and Wong, W. H. (2017).  Learning Summary Statistic For Approximate Bayesian Computation Via Deep Neural Network.  Statistica Sinica, 27(4):1595–1618.
  * Jiang et al., (2018) Jiang, B., Wu, Tung-Yu, and Wing Hung Wong (2018).  Approximate Bayesian Computation with Kullback-Leibler Divergence as Data Discrepancy.  In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, pages 1711–1721. PMLR.
  * Kingma and Welling, (2022) Kingma, D. P. and Welling, M. (2022).  Auto-Encoding Variational Bayes.
  * Kolmogorov, (1942) Kolmogorov, AN. (1942).  Definition of center of dispersion and measure of accuracy from a finite number of observations (in Russian).  Izv. Akad. Nauk SSSR Ser. Mat., 6:3–32.
  * Lázaro-Gredilla et al., (2010) Lázaro-Gredilla, M., Quinonero-Candela, J., Rasmussen, C. E., and Figueiras-Vidal, A. R. (2010).  Sparse spectrum gaussian process regression.  The Journal of Machine Learning Research, 11:1865–1881.
  * Levina and Bickel, (2001) Levina, E. and Bickel, P. (2001).  The earth mover’s distance is the mallows distance: Some insights from statistics.  In Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001, volume 2, pages 251–256. IEEE.
  * Longstaff and Schwartz, (2001) Longstaff, F. A. and Schwartz, E. S. (2001).  Valuing American options by simulation: A simple least-squares approach.  The review of financial studies, 14(1):113–147.
  * Mehta et al., (2014) Mehta, P. M., Walker, A., Lawrence, E., Linares, R., Higdon, D., and Koller, J. (2014).  Modeling satellite drag coefficients with response surfaces.  Advances in Space Research, 54(8):1590–1607.
  * Müller et al., (1997) Müller, P., West, M., and MacEachern, S. N. (1997).  Bayesian models for non-linear auto-regressions.  Journal of Time Series Analysis, 18:593–614.
  * Nareklishvili et al., (2022) Nareklishvili, M., Polson, N., and Sokolov, V. (2022).  Deep partial least squares for IV regression.  arXiv preprint arXiv:2207.02612.
  * Nunes and Balding, (2010) Nunes, M. A. and Balding, D. J. (2010).  On Optimal Selection of Summary Statistics for Approximate Bayesian Computation.  Statistical Applications in Genetics and Molecular Biology, 9(1).
  * Ostrovski et al., (2018) Ostrovski, G., Dabney, W., and Munos, R. (2018).  Autoregressive Quantile Networks for Generative Modeling.
  * Padilla et al., (2022) Padilla, O. H. M., Tansey, W., and Chen, Y. (2022).  Quantile regression with relu networks: Estimators and minimax rates.  The Journal of Machine Learning Research, 23(1):11251–11292.
  * Papamakarios and Murray, (2016) Papamakarios, G. and Murray, I. (2016).  Fast \epsilon -free Inference of Simulation Models with Bayesian Conditional Density Estimation.  In Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc.
  * Papamakarios et al., (2017) Papamakarios, G., Pavlakou, T., and Murray, I. (2017).  Masked Autoregressive Flow for Density Estimation.  In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Park et al., (2016) Park, M., Jitkrittum, W., and Sejdinovic, D. (2016).  K2-ABC: Approximate Bayesian Computation with Kernel Embeddings.  In Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, pages 398–407. PMLR.
  * Parzen, (2004) Parzen, E. (2004).  Quantile Probability and Statistical Data Modeling.  Statistical Science, 19(4):652–662.
  * Pastorello et al., (2003) Pastorello, S., Patilea, V., and Renault, E. (2003).  Iterative and recursive estimation in structural nonadaptive models.  Journal of Business & Economic Statistics, 21(4):449–509.
  * Polson and Sokolov, (2015) Polson, N. and Sokolov, V. (2015).  Bayesian analysis of traffic flow on interstate I-55: The LWR model.  The Annals of Applied Statistics, 9(4):1864 – 1888.
  * Polson et al., (2021) Polson, N., Sokolov, V., and Xu, J. (2021).  Deep learning partial least squares.  arXiv preprint arXiv:2106.14085.
  * Polson and Ročková, (2018) Polson, N. G. and Ročková, V. (2018).  Posterior concentration for sparse deep learning.  Advances in Neural Information Processing Systems, 31.
  * Polson and Sokolov, (2017) Polson, N. G. and Sokolov, V. O. (2017).  Deep learning for short-term traffic flow prediction.  Transportation Research Part C: Emerging Technologies, 79:1–17.
  * Rezende and Mohamed, (2015) Rezende, D. J. and Mohamed, S. (2015).  Variational inference with normalizing flows.  arXiv preprint arXiv:1505.05770.
  * Sauer et al., (2023) Sauer, A., Cooper, A., and Gramacy, R. B. (2023).  Non-stationary gaussian process surrogates.
  * Schmidt-Hieber, (2020) Schmidt-Hieber, J. (2020).  Nonparametric regression using deep neural networks with relu activation function.  The Annals of Statistics, 48(4):1875–1897.
  * Schultz et al., (2022) Schultz, L., Auld, J., and Sokolov, V. (2022).  Bayesian Calibration for Activity Based Models.  arXiv preprint arXiv:2203.04414.
  * Shan and Wang, (2010) Shan, S. and Wang, G. G. (2010).  Survey of modeling and optimization strategies to solve high-dimensional design problems with computationally-expensive black-box functions.  Structural and Multidisciplinary Optimization, 41(2):219–241.
  * Shen et al., (2021) Shen, G., Jiao, Y., Lin, Y., Horowitz, J. L., and Huang, J. (2021).  Deep quantile regression: Mitigating the curse of dimensionality through composition.  arXiv preprint arXiv:2107.04907.
  * Sohl-Dickstein et al., (2015) Sohl-Dickstein, J., Weiss, E. A., Maheswaranathan, N., and Ganguli, S. (2015).  Deep Unsupervised Learning using Nonequilibrium Thermodynamics.
  * Stroud et al., (2003) Stroud, J. R., Müller, P., and Polson, N. G. (2003).  Nonlinear state-space models with state-dependent variances.  Journal of the American Statistical Association, 98(462):377–386.
  * Sun et al., (2019) Sun, F., Gramacy, R. B., Haaland, B., Lawrence, E., and Walker, A. (2019).  Emulating satellite drag from large simulation experiments.  SIAM/ASA Journal on Uncertainty Quantification, 7(2):720–759.
  * van den Oord and Kalchbrenner, (2016) van den Oord, A. and Kalchbrenner, N. (2016).  Pixel RNN.  In ICML.
  * Wang et al., (2022) Wang, Y., Polson, N., and Sokolov, V. O. (2022).  Data augmentation for bayesian deep learning.  Bayesian Analysis, 1(1):1–29.
  * Weng, (2019) Weng, L. (2019).  From gan to wgan.
  * Wilson et al., (2014) Wilson, A. G., Gilboa, E., Nehorai, A., and Cunningham, J. P. (2014).  Fast kernel learning for multidimensional pattern extrapolation.  Advances in neural information processing systems, 27.
  * Yakowitz et al., (1978) Yakowitz, S., Krimmel, J., and Szidarovszky, F. (1978).  Weighted monte carlo integration.  SIAM Journal on Numerical Analysis, 15(6):1289–1300.
  * Zamo and Naveau, (2018) Zamo, M. and Naveau, P. (2018).  Estimation of the continuous ranked probability score with limited information and applications to ensemble weather forecasts.  Mathematical Geosciences, 50(2):209–234.