* [13] Y. Eliashberg and M. Gromov, Convex symplectic manifolds, Proc. of Symp. in Pure Math., 52, Part 2, Several Complex Variables and Complex Geometry, 135-162, Amer. Math. Soc., Providence, RI, 1991.
  * [14] F. Elkhadhra, J-pluripolar subsets and currents on almost complex manifolds, Math. Z., 264 (2010), 399-422.
  * [15] F. Elkhadhra, Poincaré-Lelong formula, J-analytic subsets and Lelong numbers of currents on almost complex manifolds, Bull. Sci. Math. 138 (2014), 393-405.
  * [16] P. Gauduchon, Le théorème de l’excetricit$\acute{e}$ nulle, C. R. Acad. Sci. Paris. Série A, 285 (1977), 387-390.
  * [17] D. Gilbarg and N. Trudinger, Elliptic partial differential equations of second order, 2nd ed., Berlin-Heidelberg-New York, Springer, 1983.
  * [18] P. A. Griffiths and J. Harries, Principles of Algebraic Geometry, New York, Wiley, 1978.
  * [19] M. Gromov, Pseudoholomorphic curves in symplectic manifolds, Invent. Math., 82 (1985), 307-347.
  * [20] R. Harvey and H. B. Lawson, Jr., An intrinsic characterization of Kähler manifolds, Invent. Math., 74 (1983), 169-198.
  * [21] R. Harvey and H. B. Lawson, Jr., Potential theory on almost complex manifolds, Ann. Inst. Fourier, 65 (2015), 171-210.
  * [22] L. Hörmander, $L^{2}$ estimates and existence theorems for the $\bar{\partial}$ operator, Acta Math., 113 (1965), 89-152.
  * [23] L. Hörmander, An introduction to complex analysis in several variables, third edition (revised), D. Van Nostrand Co., Inc., Princeton, N.J.-Toronto, Ont.-London, 1990.
  * [24] A. Lamari, Le cone Kählerien d’une surface, J. Math. Pures Appl., 78 (1999), 249-263.
  * [25] A. Lamari, Courants kaehleriens et surfaces compactes, Ann. Inst. Fourier, 49 (1999), 263-285.
  * [26] M. Lejmi, Strictly nearly Kähler 6-manifolds are not compatible with symplectic forms, C. R. Math. Acad. Sci. Paris, 34 (2006), 759-762.
  * [27] M. Lejmi, Stability under deformations of extremal almost-Kähler metrics in dimension 4, Math. Res. Lett., 17 (2010), 601-612.
  * [28] M. Lejmi, Extremal almost-Kähler metrics, International J. Math., 21 (2010), 1639-1662.
  * [29] T.J. Li and A. K. Liu, Symplectic structures on ruled surfaces and a generalized adjunction formula, Math. Res. Lett., 2 (1995), 453-471.
  * [30] T.-J. Li and W. Zhang, Comparing tamed and compatible symplectic cones and cohomological properties of almost complex manifolds, Comm. Anal. Geom., 17 (2009), 651-684.
  * [31] T.-J. Li and W. Zhang, Almost Kähler forms on rational 4-manifolds, Amer. J. Math., 137 (2015), 1209-1256.
  * [32] T.-J. Li and W. Zhang, J-holomorphic curves in a nef class, IMRN, 2015 (2015), 12070-12104.
  * [33] D. McDuff and D. Salamon, J-holomorphic curves and symplectic topology, Amer. Math. Soc., Providence, 2004.
  * [34] D. McDuff and D. Salamon, Introduction to symplectic topology, Oxford University Press, Third Editon, 2017.
  * [35] J. Morrow and K. Kodaira, complex manifolds, Holt. Rinehart and Winston, New York, 1971.
  * [36] N. Pali, Fonctions plurisousharmoniques et courants positifs de type $(1,1)$ sur une variété complex, Manuscripta Math., 118 (2005), 311-337.
  * [37] J. C. Sikorav, Some properties of holomorphic curves in almost complex manifolds, Holomorphic curves in symplectic geometry, 165-189, Progr. Math., 117, Birkhäuser, Basel, 1994.
  * [38] Y.-T. Siu, Analyticity of sets associated to Lelong numbers and the extension of closed positive currents, Invent. Math., 27 (1974), 53-156.
  * [39] D. Sullivan, Cycles for the dynamical study of foliated manifolds and complex manifolds, Invent. Math., 36 (1976), 225-255.
  * [40] Q. Tan, H. Y. Wang, Y. Zhang and P. Zhu, On cohomology of almost complex 4-manifolds, J. Geom. Anal., 25 (2015), 1431-1443.
  * [41] Q. Tan, H. Y. Wang and J. R. Zhou, Primitive cohomology of almost real degree two on compact symplectic manifolds, Manuscripta Math., 148 (2015), 535-556.
  * [42] Q. Tan, H. Y. Wang and J. R. Zhou, A note on the deformations of almost complex structures on closed four-manifolds, J. Geom. Anal., 27 (2017), 2700-2724.
  * [43] Q. Tan, H.Y. Wang, J.R. Zhou and P. Zhu, On tamed almost complex four manifolds, Peking Mathematical Journal, 5(1) (2022), 37-152.
  * [44] C. H. Taubes, Differential Geommetry: Bundles, Connections, Metrics and Curvatures, Oxford Graduate Texts in Mathematics 23, 2011.
  * [45] C. H. Taubes, Tamed to compatible: symplectic forms via moduli space integration, J. Symplectic Geom., 9(2) (2011), 161-250.
  * [46] K. Yosida, Functional Analysis, Sixth Edition, Springer-Verleg, 1980.
  * [47] W. Zhang, The curve cone of almost complex 4-manifolds, Proc. London Math. Soc., 115(6) (2017), 1227-1275.
  * [48] F. Zheng, Private communications, 2023.

Hongyu Wang,
School of Mathematical Sciences,
Yangzhou University,
Yangzhou, 225002,
Jiangsu, People’s Republic of China
e-mail<EMAIL_ADDRESS>
Ken Wang,
School of Mathematical Sciences,
Fudan University,
Shanghai, 100433,
People’s Republic of China
e-mail<EMAIL_ADDRESS>
Peng Zhu,
School of Mathematics and Physics,
Jiangsu University of Technology,
Changzhou, 213001,
Jiangsu, People’s Republic of China
e-mail<EMAIL_ADDRESS>


# Variational Inference with Coverage Guarantees

Yash Patel
Department of Statistics
University of Michigan
<EMAIL_ADDRESS>
&Declan McNamara
Department of Statistics
University of Michigan
<EMAIL_ADDRESS>
&Jackson Loper
Department of Statistics
University of Michigan
<EMAIL_ADDRESS>
Jeffrey Regier*
Department of Statistics
University of Michigan
<EMAIL_ADDRESS>
&Ambuj Tewari
Department of Statistics
University of Michigan
<EMAIL_ADDRESS>
Senior author order was decided by a coin flip.

###### Abstract

Amortized variational inference produces a posterior approximator that can
compute a posterior approximation given any new observation. Unfortunately,
there are few guarantees about the quality of these approximate posteriors. We
propose Conformalized Amortized Neural Variational Inference (CANVI), a
procedure that is scalable, easily implemented, and provides guaranteed
marginal coverage. Given a collection of candidate amortized posterior
approximators, CANVI constructs conformalized predictors based on each
candidate, compares the predictors using a metric known as predictive
efficiency, and returns the most efficient predictor. CANVI ensures that the
resulting predictor constructs regions that contain the truth with high
probability (exactly how high is prespecified by the user). CANVI is agnostic
to design decisions in formulating the candidate approximators and only
requires access to samples from the forward model, permitting its use in
likelihood-free settings. We prove lower bounds on the predictive efficiency
of the regions produced by CANVI and explore how the quality of a posterior
approximation relates to the predictive efficiency of prediction regions based
on that approximation. Finally, we demonstrate the accurate calibration and
high predictive efficiency of CANVI on a suite of simulation-based inference
benchmark tasks and an important scientific task: analyzing galaxy emission
spectra.

##  1 Introduction

Variational inference (VI) has become a staple in Bayesian inference; however,
it has been repeatedly noted that a major shortcoming of VI is its lack of any
theoretical guarantees and tendency to produce biased posterior estimates [1,
2, 3, 4]. For instance, in the subfield of likelihood-free inference, a meta-
study on widely used algorithms relying on VI revealed that they consistently
produce unfaithful, overconfident posterior approximations [5]. Researchers
have begun to address this shortcoming. For instance, the authors of [6]
proposed a method of correction, but noted that their method “should not be
viewed as a way to obtain conservative posterior estimators with 100%
reliability.” In particular, this method has poor scaling to high-dimensional
problems and only has theoretical guarantees of coverage in expectation, which
are only enforced heuristically.

Variational inference constructs an approximation of the full posterior
distribution; however, in applications where many posteriors need to be
estimated, constructing credible regions with marginally calibrated coverage
can be sufficient for downstream scientific inquiries. For instance, in
astrophysics, there is great interest in constraining the $\Lambda$CDM model,
the current concordance model in cosmology, through large-scale cosmological
observations [7, 8, 9, 10, 11]. A recent work from this community, [12],
leveraged Bayesian inference towards this end, obtaining approximate
posteriors for the parameters of interest on each of 10,000 observations.
Crucially, these posteriors were then solely used to produce credible
intervals with marginally valid frequentist coverage, from which they made
subsequent claims on constraints on the $\Lambda$CDM model.

Figure 1:  CANVI is a wrapper around variational inference requiring minimal
implementation and computational overhead that produces prediction regions
with guaranteed marginal calibration. Among a family of candidate amortized
posterior approximators, CANVI can identify the approximator leading to the
most efficient prediction regions. CANVI can be used in any setting where the
forward model $\mathcal{P}(X\mid\Theta)$ can be sampled.

Our insight is that conformal prediction can be leveraged to provide
variational approximators with marginal coverage guarantees and provide users
new ways to measure the quality of variational approximators. In this
manuscript, we present CANVI (Conformalized Amortized Neural Variational
Inference), a novel, general framework for producing marginally calibrated,
informative prediction regions from a collection of variational approximators.
Such regions can be produced with minimal implementation and computational
overhead, requiring only samples from the prior and
$\mathcal{P}(X\mid\Theta)$, as shown in Figure 1. In Section 3.2, we provide
theoretical analysis of the informativeness of the prediction regions produced
by CANVI using a measure known as “predictive efficiency.” High predictive
efficiency is necessary to draw conclusions in downstream scientific
inquiries. Finally, in Section 4, we show calibration and predictive
efficiency empirically across simulation-based inference benchmark tasks and
an important scientific task: analyzing galaxy emission spectra.

##  2 Background

###  2.1 Variational Inference

Bayesian methods aim to sample the posterior distribution
$\mathcal{P}(\Theta\mid X)$, typically using either MCMC or VI. VI has risen
in popularity recently due to how well it lends itself to amortization. Given
an observation $X$, variational inference transforms the problem of posterior
inference into an optimization problem by seeking a minimizer

$\varphi^{*}(X)=\operatorname*{arg\,min}_{\varphi}D(q_{\varphi}(\Theta)||\mathcal{P}(\Theta\mid
X)),$ (1)

where $D$ is a divergence and $q_{\varphi}$ is a member of a variational
family of distributions $\mathcal{Q}$ indexed by the free parameter $\varphi$.
Normalizing flows have emerged as a particularly apt choice for $\mathcal{Q}$,
as they are highly flexible and perform well empirically [13, 14]. Amortized
variational inference expands on this approach by training a neural network to
approximate $\varphi^{*}(X)$. This leads to a variational posterior
approximator $q(\Theta\mid X)=q_{\varphi^{*}(X)}(\Theta)$ that can be rapidly
computed for any value $X$. The characteristics of $\varphi^{*}$ depend in
part on the variational objective, $D$. For instance, using a reverse-KL
objective, i.e. $D_{KL}(q_{\varphi}(\Theta)||\mathcal{P}(\Theta\mid X))$, is
known to produce mode-seeking posterior approximations, whereas using a
forward-KL objective, i.e. $D_{KL}(\mathcal{P}(\Theta\mid
X)||q_{\varphi}(\Theta))$, encourages mode-covering behavior [15]. Popular
variational objectives include the Forward-Amortized Variational Inference
(FAVI) objective [16] (also known as the “sleep-phase” objective in Reweighted
Wake-Sleep [17]), the Evidence Lower Bound (ELBO), and the Importance Weighted
ELBO (IWBO) [18].

###  2.2 Conformal Prediction

Given a dataset
$\mathcal{D_{C}}=\\{(X_{1},\theta_{1}),\ldots(X_{N_{\mathcal{C}}},\theta_{N_{\mathcal{C}}})\\}$
of i.i.d. observations from a distribution $\mathcal{P}(\Theta,X)$, conformal
prediction [19, 20] produces prediction regions with distribution-free
theoretical guarantees. A prediction region is a mapping from observations of
$X$ to sets of possible values for $\Theta$. A prediction region $\mathcal{C}$
is said to be marginally calibrated at the $1-\alpha$ level if
$\mathcal{P}(\Theta\notin\mathcal{C}(X))\leq\alpha$.

Split conformal is one popular version of conformal prediction. In this
approach, marginally calibrated regions $\mathcal{C}$ are designed using a
“score function” $s(x,\theta)$. Intuitively, the score function should have
the quality that $s(x,\theta)$ is smaller when it is more reasonable to guess
that $\Theta=\theta$ given the observation $X=x$. For example, if one has
access to a function $\hat{f}(x)$ which attempts to predict $\Theta$ from $X$,
one might take $s(x,\theta)=\|\hat{f}(x)-\theta\|$. The score function is
evaluated on each point of the dataset $\mathcal{D_{C}}$, called the
“calibration dataset,” yielding
$\mathcal{S}=\\{s(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}$. Note
that the calibration dataset cannot be used to pick the score function; if
data is used to design the score function, it must independent of
$\mathcal{D_{C}}$. This is how “split conformal” gets its name: in typical
cases, data are split into two parts, one used to design $s$ and the other to
perform calibration. We then define $\widehat{q}(\alpha)$ as the
$\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil/N_{\mathcal{C}}$ quantile of
$\mathcal{S}$. For any future $x$, the set $\mathcal{C}(x)=\\{\theta\mid
s(x,\theta)\leq\widehat{q}(\alpha)\\}$ satisfies
$1-\alpha\leq\mathcal{P}(\Theta\in\mathcal{C}(X))$. This inequality is known
as the coverage guarantee, and it arises from the exchangeability of the score
of a future test point $s(x^{\prime},\theta^{\prime})$ with $\mathcal{S}$.
Those new to conformal inference may be surprised to note that the coverage
guarantee holds regardless of the number of samples $N_{\mathcal{C}}$ used in
calibration; conformal guarantees are not asymptotic results.

As noted in Vovk’s tutorial [20], while the coverage guarantee holds for any
score function, different score functions may lead to more or less informative
prediction regions. For example, the score $s(x,\theta)=1$ leads to the highly
uninformative prediction region of all possible values of $\Theta$. Predictive
efficiency is one way to quantify informativeness [21, 22]. It is defined as
the inverse of the expected Lebesgue measure of the prediction region, i.e.
$\left(\mathbb{E}[|\mathcal{C}(X)|]\right)^{-1}$. Methods employing conformal
prediction often seek to identify prediction regions that are efficient as
well as marginally calibrated.

##  3 Method

We now propose a way to produce efficient prediction regions with coverage
guarantees using a collection of amortized VI approximators
$\\{q^{(1)}(\Theta\mid X),\ldots q^{(T)}(\Theta\mid X)\\}$. CANVI can be
applied whenever it is possible to sample from the joint distribution
$\mathcal{P}(X,\Theta)$. The coverage validity of CANVI is proven in Section
3.1, and analyses of its predictive efficiency follow in Section 3.2.

###  3.1 CANVI

In the simplest case, CANVI takes as input a single amortized posterior
approximator $q(\Theta\mid X)$ and computes an empirical score distribution by
drawing
$\mathcal{D_{C}}=\\{(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}\overset{\mathrm{iid}}{\sim}\mathcal{P}(\Theta)\mathcal{P}(X\mid\Theta)$.
In traditional applications of split conformal, much concern is given to the
loss of accuracy of the predictor $\widehat{f}$ in having to reserve a subset
of the training data for calibration. Here we have no such issues; we can
sample from the joint distribution to produce a calibration dataset that is
arbitrarily large. Given $q(\Theta\mid X)$, we employ the following score, as
used in [19, 23]:

$s(x_{i},\theta_{i})=\frac{1}{q(\theta_{i}\mid x_{i})}.$ (2)

For a desired coverage of $1-\alpha$, we denote the
$\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil/N_{\mathcal{C}}$ quantile of the
empirical score distribution over $\mathcal{D}_{C}$ as
$\widehat{q}(\mathcal{D}_{\mathcal{C}},\alpha)$. The prediction region
$\mathcal{C}(x)=\left\\{\theta:1/q(\theta\mid
x)\leq\widehat{q}(\mathcal{D}_{\mathcal{C}},\alpha)\right\\}$ is then
marginally calibrated. Note that $\mathcal{C}(x)$ may be disjoint if the
posterior is multimodal.

This approach achieves coverage guarantees for any variational approximator
$q(\Theta\mid X)$. However, a poorly chosen approximator may lead to poor
predictive efficiency. To mitigate the risk of producing an inefficient
predictor, it is, therefore, natural to explore multiple posterior
approximations, $\\{q^{(t)}(\Theta\mid X)\\}_{t=1}^{T}$, where these posterior
approximations could, for instance, differ in their choice of training
objective, variational family, or hyperparameters.

CANVI assesses a collection of candidate variational approximators as follows.
For each candidate, CANVI uses a test dataset
$\mathcal{D_{T}}=\\{(x_{i},\theta_{i})\\}_{i=1}^{N_{\mathcal{T}}}\overset{\mathrm{iid}}{\sim}\mathcal{P}(\Theta)\mathcal{P}(X\mid\Theta)$
and calculates the empirical inverse efficiencies for each candidate
$q^{(t)}(\Theta\mid X)$, as follows:

$\widehat{\ell}(q^{(t)},\alpha):=\frac{1}{N_{\mathcal{T}}}\sum_{i=1}^{N_{\mathcal{T}}}\mathcal{L}\left(\left\\{\theta:1/q^{(t)}(\theta\mid
x_{i})\leq\widehat{q}^{(t)}(\mathcal{D}_{\mathcal{C}},\alpha)\right\\}\right),$
(3)

where $\mathcal{L}(\cdot)$ is the Lebesgue measure. CANVI identifies the
$t^{*}$ such that $q^{(t^{*})}$ leads to the lowest empirical inverse
efficiency, which we denote as $q^{(*)}:=q^{(t^{*})}$. For each $x_{i}$,
$\mathcal{L}(\mathcal{C}(x_{i}))$ is empirically estimated for each $q^{(t)}$
using an importance-weighted Monte Carlo estimate over $S$ samples, namely

$\mathbb{E}_{q^{(t)}(\Theta_{1:S}\mid
x_{i})}\left[\frac{1}{S}\sum_{j=1}^{S}\frac{1}{q^{(t)}(\Theta_{j}\mid
x_{i})}\mathbbm{1}\left[1/q^{(t)}(\Theta_{j}\mid
x_{i})\leq\widehat{q}^{(t)}\right]\right].$ (4)

Note that this estimate can alternatively be performed using a grid-
discretization over $\text{Supp}(\Theta\mid x_{i})$. However, doing so is only
feasible in low-dimensional cases or where the support has a known, small
extent. We empirically demonstrate that this estimate captures efficiency
behaviors that parallel those obtained from explicit gridding in Section
4.1.1. After selecting $t^{*}$, an additional recalibration step must then be
performed. CANVI uses an additional calibration dataset $\mathcal{D_{R}}$,
which we take to be the same size as $\mathcal{D_{C}}$, i.e.
$|\mathcal{D_{R}}|=N_{\mathcal{C}}$. CANVI constructs $\mathcal{D_{R}}$ with
i.i.d. draws from $\mathcal{P}(X,\Theta)$ and computes the quantile
$\widehat{q}^{(*)}:=\widehat{q}^{(t^{*})}(\mathcal{D}_{\mathcal{R}},\alpha)$.
Using $(q^{(*)}(\Theta\mid X),\widehat{q}^{(*)})$ to construct prediction
regions produces provably efficient predictions, as we show theoretically in
Section 3.2. The additional recalibration is necessary to retain the
exchangeability between scores of future test points
$s(x^{\prime},\theta^{\prime})$ and $\mathcal{S}$, as exchangeability is lost
in conditioning on $\mathcal{D_{C}}$ for selecting $t^{*}$.

The full CANVI framework is provided in Algorithm 1. In summary, CANVI chooses
a posterior approximator $q^{(*)}(\Theta\mid X)$ from a collection of
posterior approximators $\\{q^{(t)}(\Theta\mid X)\\}_{t=1}^{T}$ and produces a
corresponding quantile $\widehat{q}^{(*)}$ for a target coverage probability
$(1-\alpha)$. CANVI implicitly defines a prediction region map defined by
$\mathcal{C}(x)=\left\\{\theta:1/q^{(*)}(\theta\mid
x)\leq\widehat{q}^{(*)}\right\\}$ that is marginally calibrated at the
$(1-\alpha)$ level. Note that, while CANVI uses Equation 2 as its score, the
use of strictly monotonically increasing transforms of this score produces
identical coverage regions, with a special case being the absolute z-score for
a unimodal Gaussian variational family, proven in Appendix A.

Algorithm 1 CANVI (Conformalized Amortized Neural Variational Inference)

1:procedure CPQuantile

2:Inputs:  Posterior approximation $q(\Theta\mid X)$, Calibration set
$\mathcal{D_{C}}$, Desired coverage $1-\alpha$

3: $\mathcal{S}\leftarrow\\{\frac{1}{q(\theta_{i}\mid
x_{i})}\\}_{i=1}^{N_{\mathcal{C}}}$ $\triangleright$
$N_{\mathcal{C}}=|\mathcal{D_{C}}|$

4: Return
$\frac{\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil}{N_{\mathcal{C}}}\text{
quantile of }\mathcal{S}$

5:end procedure

6:

7:procedure CPRegionSize

8:Inputs:  Posterior approximation $q(\Theta\mid X)$, CP quantile
$\widehat{q}$, Test set $\mathcal{D_{T}}$

9:
$\widehat{\ell}_{i}\leftarrow\left\\{\frac{1}{S}\sum_{j=1}^{S}\frac{1}{q(\theta_{ij}\mid
x_{i}^{\prime})}\mathbbm{1}\left[1/q(\theta_{ij}\mid
x_{i}^{\prime})\leq\widehat{q}\right]\right\\}_{i=1}^{N_{\mathcal{T}}}$
$\triangleright$ $\Theta_{i,j}\sim q_{\varphi}(\Theta\mid x_{i}^{\prime})$

10: Return
$\frac{1}{N_{\mathcal{T}}}\sum_{i=1}^{N_{\mathcal{T}}}\widehat{\ell}_{i}$
$\triangleright$ $N_{\mathcal{T}}=|\mathcal{D_{T}}|$

11:end procedure

12:

13:procedure CANVI

14:Inputs:  Posterior approximations $\\{q^{(t)}(\Theta\mid X)\\}_{t=1}^{T}$,
Prior $\mathcal{P}(\Theta)$, Forward model $\mathcal{P}(X\mid\Theta)$, Desired
coverage $1-\alpha$, Calibration set size $N_{\mathcal{C}}$, Test set size
$N_{\mathcal{T}}$

15:
$\mathcal{D_{C}}\leftarrow\\{\theta_{i}\sim\mathcal{P}(\Theta),x_{i}\sim\mathcal{P}(X\mid\theta_{i})\\}_{i=1}^{N_{\mathcal{C}}}.$

16:
$\mathcal{D_{T}}\leftarrow\\{\theta_{i}\sim\mathcal{P}(\Theta),x_{i}\sim\mathcal{P}(X\mid\theta_{i})\\}_{i=1}^{N_{\mathcal{T}}}.$

17: for $t\in\\{1,\ldots T\\}$ do

18: $\widehat{q}^{(t)}\leftarrow\textsc{CPQuantile}(q^{(t)}(\Theta\mid
X),\mathcal{D_{C}},1-\alpha)$

19: $\widehat{\ell}^{(t)}\leftarrow\textsc{CPRegionSize}(q^{(t)}(\Theta\mid
X),\widehat{q}^{(t)},\mathcal{D_{T}})$.

20: end for

21: $t^{*}\leftarrow\operatorname*{arg\,min}_{t}\widehat{\ell}^{(t)}$

22:
$\mathcal{D_{R}}\leftarrow\\{\theta_{i}\sim\mathcal{P}(\Theta),x_{i}\sim\mathcal{P}(X\mid\theta_{i})\\}_{i=1}^{N_{\mathcal{C}}}.$

23: $\widehat{q}^{(*)}\leftarrow\textsc{CPQuantile}(q^{(t^{*})}(\Theta\mid
X),\mathcal{D_{R}},1-\alpha)$.

24: Return $q^{(*)}(\Theta\mid X)$, $\widehat{q}^{(*)}$ $\triangleright$
$q^{(*)}(\Theta\mid X):=q^{(t^{*})}(\Theta\mid X)$

25:end procedure

######  Lemma 3.1.

Let $q^{(*)}(\Theta\mid
X),\widehat{q}^{(*)}=\mathrm{CANVI}\left(\\{q^{(t)}(\Theta\mid
X)\\}_{t=1}^{T},\mathcal{P}(\Theta,X),1-\alpha,N_{\mathcal{C}},N_{\mathcal{T}}\right).$
Suppose $(x^{\prime},\theta^{\prime})\sim\mathcal{P}(\Theta,X)$ is drawn
independently of
$\mathcal{D}\cup\mathcal{D_{C}}\cup\mathcal{D_{T}}\cup\mathcal{D_{R}}$, where
$\mathcal{D}$ is the aggregate of training data used for
$\\{q^{(t)}(\Theta\mid X)\\}_{t=1}^{T}$. Then
$1-\alpha\leq\mathcal{P}(1/q^{(*)}(\theta^{\prime}\mid
x^{\prime})\leq\widehat{q}^{(*)})$.

###### Proof.

Consider the score function $s^{(*)}(x,\theta):=1/q^{(*)}(\theta\mid x)$.
Observe that $(x^{\prime},\theta^{\prime})\cup\mathcal{D_{R}}$ are jointly
sampled i.i.d. from $\mathcal{P}(\Theta,X)$, independent of the datasets used
to design $s^{(*)}(x,\theta)$, namely
$\mathcal{D}\cup\mathcal{D_{C}}\cup\mathcal{D_{T}}$. Denoting
$\mathcal{S}_{\mathcal{R}}:=\\{s^{(*)}(x_{i},\theta_{i})\\}_{(x_{i},\theta_{i})\in\mathcal{D_{R}}}$,
scores $s^{(*)}(x^{\prime},\theta^{\prime})\cup\mathcal{S}_{\mathcal{R}}$,
thus, too are i.i.d and, hence, exchangeable. The coverage guarantee then
follows from the general theory of conformal prediction, presented in [19]. ∎

###  3.2 Predictive Efficiency of CANVI

We now wish to characterize the efficiency of the prediction regions produced
by CANVI. To do so, we first demonstrate that the choice $q^{(*)}$ is, with
high probability, oracle optimal.

We define the oracle action (for a fixed $\alpha$) as returning the posterior
approximator-quantile pair that minimizes Equation 3 over observed pairs
$\\{(q^{(t)}(\Theta\mid
X),\widehat{q}^{(t)}(\mathcal{D}_{\mathcal{C}},\alpha))\\}_{t=1}^{T}$. Recall
that CANVI performs a recalibration after selecting $t^{*}$, where importantly
$\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{R}},\alpha)$ is defined using
$\mathcal{D}_{\mathcal{R}}$. Unfortunately, $(q^{(*)}(\Theta\mid
X),\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{R}},\alpha))$ may not be the
minimizer of Equation 3 over this extended collection of pairs
$\\{(q^{(t)}(\Theta\mid
X),\widehat{q}^{(t)}(\mathcal{D}_{\mathcal{C}},\alpha))\\}_{t=1}^{T}\cup(q^{(*)}(\Theta\mid
X),\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{R}},\alpha))$. CANVI could return
$(q^{(*)}(\Theta\mid X),\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{C}},\alpha))$
to ensure efficiency, but the recalibration is necessary to ensure coverage
guarantees. This tradeoff between coverage and efficiency was studied in [21],
from which we obtain guarantees on $t^{*}$ remaining oracle optimal with high
probability.

Intuitively, $(q^{(*)}(\Theta\mid
X),\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{R}},\alpha))$ will remain close to
oracle optimal if the score distributions and, hence, quantiles of
$q^{(*)}(\Theta\mid X)$ under $\mathcal{D}_{\mathcal{C}}$ and
$\mathcal{D}_{\mathcal{R}}$ are close, and if, for any deviations in the
quantiles, the predictive efficiency varies smoothly, respectively formalized
from [21] as follows.

######  Assumption 1.

For each $t$, $\exists$ $r^{*},\gamma\in(0,1)$ such that the inverse CDF of
the score under $q^{(t)}$, $F_{t}^{-1}$, is $\gamma$-Hölder continuous on
[$\widehat{q}^{(t)}(\alpha)-r^{*},\widehat{q}^{(t)}(\alpha)+r^{*}$] with
Hölder continuity constant $L_{t}$.

######  Assumption 2.

For each $t$, the map $\alpha\rightarrow\widehat{\ell}(q^{(t)},\alpha)$ is
Lipschitz continuous with constant $L_{W}$.

The proof in [21] is completed by demonstrating that the CDFs of the score
distributions of $q^{(*)}(\Theta\mid X)$ under $\mathcal{D}_{\mathcal{C}}$ and
$\mathcal{D}_{\mathcal{R}}$, and hence
$\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{C}},\alpha)$ and
$\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{R}},\alpha)$ by Assumption 1, are
close and concluding by then bounding resulting deviations in predictive
efficiencies by Assumption 2. To make explicit the distinction in computed
efficiency, we denote the recalibrated efficiency as
$\widehat{\ell}_{\mathcal{R}}(q^{(*)},\alpha)$, where
$\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{R}},\alpha)$ is used in place of the
standard $\widehat{q}^{(*)}(\mathcal{D}_{\mathcal{C}},\alpha)$. The proof of
Theorem 3.2 is explicitly provided in Appendix B.

######  Theorem 3.2.

Let $q^{(*)}(\Theta\mid
X),\widehat{q}^{(*)}=\mathrm{CANVI}\left(\\{q^{(t)}(\Theta\mid
X)\\}_{t=1}^{T},\mathcal{P}(\Theta,X),1-\alpha,N_{\mathcal{C}},N_{\mathcal{T}}\right).$
Let $\mathcal{D},\mathcal{D_{C}},\mathcal{D_{T}}$, and $\mathcal{D_{R}}$ be
drawn i.i.d. from $\mathcal{P}(\Theta,X)$, where $\mathcal{D}$ is the
aggregate of training data used for $\\{q^{(t)}(\Theta\mid X)\\}_{t=1}^{T}$.
Take $\delta\in[0,1]$. If Assumption 2 holds and Assumption 1 holds with

$r^{*}\geq\max\left\\{\sqrt{\frac{\log(4T/\delta)}{2N_{\mathcal{C}}}},\frac{2}{N_{\mathcal{C}}}\right\\},$
(5)

then with probability at least $1-\delta$,

$\widehat{\ell}_{\mathcal{R}}(q^{(*)},\alpha)\leq\min_{1\leq t\leq
T}\widehat{\ell}(q^{(t)},\alpha)+3L_{W}L_{[T]}\left[\left(\frac{\log(4T/\delta)}{N_{\mathcal{C}}}\right)^{\gamma/2}+\left(\frac{2}{N_{\mathcal{C}}}\right)^{\gamma}\right],$
(6)

where $\gamma$, $L_{W}$, and $L_{[T]}=\max_{1\leq t\leq T}L_{t}$ are constants
defined in Assumptions 1 and 2.

Again, in any setting where it is possible to sample from
$\mathcal{P}(X,\Theta)$, $N_{\mathcal{C}}$ can be made arbitrarily large,
achieving a tight bound in Equation 6 as a consequence. The implication of
Theorem 3.2 is that practitioners can focus on obtaining efficient predictors,
as discussed in the following section, knowing that CANVI will make the
optimal selection from this set with high probability.

###  3.3 Optimally Efficient Predictor Characterization

We now wish to study how the qualities of a posterior approximation relate to
the efficiency of its conformalization. Importantly, this characterization is
independent of algorithmic decisions made by CANVI. For this reason, we
introduce a population-level version of the inverse-efficiency, defined as

$\ell(q,\alpha):=\mathbb{E}_{X}\left[\mathcal{L}\left(\left\\{\theta:1/q(\Theta\mid
X)\leq\widehat{q}(\alpha)\right\\}\right)\right],$ (7)

where $\widehat{q}(\alpha)$ is the exact $1-\alpha$ quantile for
$1/q(\Theta|X)$. We begin by observing that $q(\Theta\mid X)$ with incorrect
support can lead to the most inefficient prediction region possible, i.e.,
$\text{Dom}(\Theta)$.

######  Theorem 3.3.

Suppose $\forall x\in\text{Dom}(X)$, there exists a region
$\vartheta\subset\text{Dom}(\Theta)$ s.t. $\mathcal{P}(\vartheta\mid
x)\geq\beta$ but $q(\vartheta\mid x)=0$. Then, for any $\alpha\leq\beta$,
$\ell(q,\alpha)=|\text{Dom}(\Theta)|$.

###### Proof.

We wish to find the distribution of $s(X,\Theta)=1/q(\Theta\mid X)$ jointly
over $X,\Theta$. By assumption, for any draw $\theta\in\vartheta$,
$s(x,\theta)=\infty$. Thus,

$\mathcal{P}(s(X,\Theta)=\infty)=\int\mathcal{P}(s(X,\Theta)=\infty\mid
X=x)\mathcal{P}(X\in dx)\geq\beta\int\mathcal{P}(X\in dx)=\beta.$

Thus, $\forall\alpha\leq\beta$,
$\widehat{q}(\alpha)=\infty\implies\left\\{\theta:1/q(\theta\mid
X)\leq\widehat{q}(\alpha)\right\\}=\text{Dom}(\Theta)$. Since this holds
$\forall x$, it follows that $\ell(q,\alpha)=|\text{Dom}(\Theta)|$. ∎

One implication of Theorem 3.3 is preference should be given to forward-KL
over reverse-KL training objectives to encourage mode coverage, as we
additionally demonstrate empirically in Section 4.2.

We next consider whether the _true_ posterior leads to the greatest
efficiency. Section 4.1.2 empirically suggests that efficiency generally
improves with convergence to the exact posterior. This leads us to conjecture
that the true posterior leads to the greatest average efficiency over all
possible $\alpha$.

######  Conjecture 1.

Let $(X,\Theta)\sim\mathcal{P}(\Theta,X)$, $q^{(*)}(\Theta\mid
X)=\mathcal{P}(\Theta\mid X)$ and $q(\Theta\mid X)$ denote any other posterior
approximator. Then

$\int_{0}^{1}\ell(q^{(*)},\alpha)d\alpha\leq\int_{0}^{1}\ell(q,\alpha)d\alpha.$
(8)

If true, Conjecture 1 suggests practitioners should seek to use more flexible
variational families to aim to recover the exact posterior to maximize
efficiency. However, we note that for any particular $\alpha$, the true
posterior does not always lead to the greatest efficiency, as demonstrated in
Appendix C. We conclude by considering a family of distributions within which
the true posterior leads to the highest efficiency. Future work could be
directed toward extending this analysis to a broader set of families. Details
and experimental verification for this Gaussian case are given in Appendix D.

######  Theorem 3.4.

Let $\Theta$ and $X$ be zero-mean unit-variance Gaussian random variables with
correlation $\rho$. Let $q_{\varphi}(\theta|x)=\mathcal{N}(\theta;\varphi
x,1-\rho^{2})$. Then $q_{\rho}$ gives the true posterior and
$\ell(q_{\rho},\alpha)=\min_{\varphi}\ell(q_{\varphi},\alpha)$.

###### Proof.

We begin by finding the distribution of
$s_{\varphi}(X,\Theta)=1/q_{\varphi}(\Theta\mid X)$. Explicitly computing its
distribution yields a closed form of the threshold as

$\widehat{q}_{\varphi}=\sqrt{2\pi(1-\rho^{2})\exp\left(\frac{\varphi^{2}+1-2\varphi\rho}{1-\rho^{2}}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right)^{2}\right)}$
(9)

Considering an arbitrary fixed $X=x$, prediction intervals are of length

$\ell(q_{\varphi},\alpha)=2\sqrt{\varphi^{2}+1-2\varphi\rho}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right)$
(10)

Observe that $\ell(q_{\varphi},\alpha)$ is an increasing function of
$|\varphi-\rho|$ with a minimum achieved at the exact posterior of
$\varphi=\rho$. Given this holds for any $\alpha$,
$\int_{0}^{1}\ell(q_{\varphi},\alpha)d\alpha$ too is minimized at
$\varphi=\rho$. ∎

##  4 Experiments

We now demonstrate calibration and study the predictive efficiency of CANVI
empirically across several tasks. We investigate predictive efficiency
specifically with respect to convergence of variational posteriors to the
exact posteriors in Section 4.1.2 and training objectives in Section 4.2.

In all experiments, coverage is assessed both by directly using the
variational approximation and its conformalized wrapper. For the former, we
construct the highest density credible region per $x_{i}$, namely by
estimating the $\zeta$ such that $\\{\theta_{j}\mid q(\theta_{j}\mid
x_{i})\geq\zeta\\}$ captures $1-\alpha$ of the probability mass. To estimate
$\widehat{\zeta}$, we draw $\\{\theta_{j}\\}_{j=1}^{N}\sim q(\theta_{j}\mid
x_{i})$ and find the $1-\alpha$ quantile of $\\{q(\theta_{j}\mid
x_{i})\\}_{j=1}^{N}$. Assessing coverage of the true parameter $\theta$ can,
therefore, be done simply by checking if $q(\theta\mid
x_{i})\geq\widehat{\zeta}$. Experimental details are provided in Appendix F,
and code will be made public upon acceptance.

###  4.1 SBI Benchmarks

We follow the precedent set in [5] and evaluate CANVI on the following
simulation-based inference (SBI) benchmark tasks: Two Moons, Lotka-Volterra,
Gaussian Mixture, Susceptible-Infected-Recovered (SIR), Simple Likelihood
Complex Posterior (SLCP) Distractors, and Bernoulli GLM Raw. For full
descriptions of the priors and forward models, refer to Appendix E. Tasks in
this benchmark suite are all likelihood-free, so all posteriors were trained
against $\mathcal{L}_{\text{FAVI}}$ using a Neural Spline Flow variational
family [24].

####  4.1.1 Coverages

Calibration was performed using 10,000 i.i.d. samples and was completed in
under one second for each task. Coverage of $q^{(t)}(\Theta\mid X)$ was
assessed at every 2000 training steps over $\alpha\in[0,1]$ discretized at
steps of .05. Non-conformalized regions ($\widehat{\zeta}$) were estimated
empirically from batches of 100 i.i.d. samples per point. Coverage was
assessed over 10 batches of 10,000 i.i.d. test samples.

Figure 2:  Calibration on the SBI benchmarks with and without
conformalization, respectively the solid and dashed lines. Conformalized lines
are slightly difficult to distinguish, as they all lie along the desired $y=x$
curve. Error bars from coverage assessments across test batches are plotted,
although they are difficult to see due to the low variance between estimates
across batches.

Figure 2 demonstrates the miscalibration of the approximate posteriors and
subsequent conformal correction. As expected, calibration of $q(\Theta\mid X)$
improves with training; however, miscalibration persists even after
convergence. In particular, while $\mathcal{L}_{\text{FAVI}}$-trained
posterior approximations do well on the Two Moons and Bernoulli GLM Raw tasks,
they are noticeably miscalibrated on the Lotka-Volterra, Gaussian Mixture, and
SIR tasks, where conformalizations rectify the situations.
$\mathcal{L}_{\text{FAVI}}$-training completely fails to learn the SLCP
Distractors task, which can nonetheless be calibrated but produces totally
uninformative prediction regions as a result.

####  4.1.2 Predictive Efficiency

We study the predictive efficiencies over the training of $q_{\varphi}$ in
Figure 3, specifically for the marginal posteriors of
$\tilde{\Theta}=(\theta_{1},\theta_{2})$ for the Two Moons, SLCP, Gaussian
Linear Uniform, and Bernoulli GLM Raw tasks, all of which suggest improving
efficiency with convergence to the exact posterior.
$\widehat{\ell}(q_{\varphi},\alpha)$ was estimated for 5 batches of 100 test
points ($|\mathcal{D_{T}}|=100$) at every 500 training steps using 10,000
importance-weighted i.i.d. samples from $q^{(t)}(\Theta\mid X)$ per point. As
previously noted, a grid-discretized estimate is possible only in low-
dimensional cases, which we show has similar trends in Appendix I. Estimation
was performed for a fixed $1-\alpha=0.95$. We additionally visualize the
credible regions over training iterates for multiple tasks in Appendix G.

Figure 3:  Estimation of $\widehat{\ell}(q_{\varphi},\alpha)$ for
$1-\alpha=0.95$, suggesting efficiency improves as $q_{\varphi}(\Theta\mid
X)\rightarrow\mathcal{P}(\Theta\mid X)$. Error bars from estimating
$\widehat{\ell}(q_{\varphi},\alpha)$ across test batches are plotted.

###  4.2 CANVI Over Training Objectives

We consider a lag-one Autoregressive Conditional Heteroskedasticity (ARCH)
model, where the exact likelihood is available [25]. The data are generated as

$y^{(t)}=\theta_{1}y^{(t-1)}+e^{(t)}\text{ and }\quad
e^{(t)}=\xi^{(t)}\sqrt{0.2+\theta_{2}\left(e^{(t-1)}\right)^{2}},$ (11)

where $y^{(0)}=0$, $e^{(0)}=0$, $T=100$, and the $\xi^{(t)}$ are independent
standard normal random variables. For further details, consult Appendix E.9.
Given a set of $n=100$ realizations of $(\theta,y_{1:T})$ from the generative
model, an amortized Neural Spline Flow (NSF) [24] was fit according to the
FAVI, ELBO, and IWBO ($K=10$) objectives. Exact and approximate posteriors are
visualized in Appendix H.7.

We examine how CANVI affects coverage rates both numerically and pictorially
over generated test datasets $\\{(\theta,y_{1:T})\\}$. Table 1 shows that
prior to conformalization, the variational posteriors are generally
miscalibrated: training by the ELBO or IWBO results in significant under-
coverage. As previously mentioned, targeting the ELBO is known to find
solutions that are mode-seeking. While the variational posterior obtained by
FAVI is better calibrated, it still needs correction. After applying CANVI,
Table 1 (right) shows that the conformalized $1-\alpha$ highest density
regions are nearly perfectly calibrated. Importantly, correction by CANVI can
result in either larger or smaller $1-\alpha$ regions, depending on the
direction of miscalibration. In settings where the variational posterior is
overdispersed, applying CANVI results in statistical coverage guarantees and
smaller $1-\alpha$ density regions. We show explicit examples of this in
Appendix H.7. Table 1 (right) additionally demonstrates using the better
calibrated $\mathcal{L}_{\text{FAVI}}$-trained approximation results in higher
predictive efficiency compared to the $\mathcal{L}_{\text{ELBO}}$\- and
$\mathcal{L}_{\text{IWBO}}$-trained counterparts, as anticipated.

Table 1: Coverage rates and standard errors for $\theta$ before (left) and
after (right) conformalization by CANVI, assessed by checking for inclusion of
$\theta$ in the $1-\alpha$ highest density region. Non-conformalized regions
were estimated empirically from batches of 1000 i.i.d. samples per point.

$1-\alpha$ | ELBO | IWBO | FAVI
---|---|---|---
0.50 | 0.1068 (0.0094) | 0.7122 (0.0128) | 0.5134 (0.0157)
0.75 | 0.2444 (0.0125) | 0.8109 (0.0092) | 0.7575 (0.0103)
0.90 | 0.3830 (0.0151) | 0.8668 (0.0053) | 0.9046 (0.0076)
0.95 | 0.4312 (0.0144) | 0.8905 (0.0068) | 0.9541 (0.0031)
$1-\alpha$ | ELBO | IWBO | FAVI
---|---|---|---
0.50 | 0.4930 (0.0229) | 0.4917 (0.0218) | 0.4945 (0.0104)
0.75 | 0.7526 (0.0111) | 0.7445 (0.0154) | 0.7411 (0.0114)
0.90 | 0.9035 (0.0117) | 0.8978 (0.0068) | 0.8995 (0.0081)
0.95 | 0.9509 (0.0084) | 0.9479 (0.0065) | 0.9500 (0.0082)
Eff | 1.9155 | 1.9036 | 0.6749

###  4.3 Galaxy Spectral Energy Distributions

We now present the application of CANVI to an important scientific problem.
The spectrum of an astronomical object is measured via a spectrograph, which
records the flux (light intensity per unit area, time, wavelength) across a
large grid of wavelength values [26, 27]. We construct a model of galaxy
spectra based on the Probabilistic Value-Added Bright Galaxy Survey simulator
(PROVABGS), which maps $\theta\in\mathbb{R}^{11}$ to galaxy spectra.
Additional details of the simulation setup and example draws from the
simulator are available in Appendix J.

A mixture of 20 Gaussian distributions was used as the variational posterior
and trained against the $\mathcal{L}_{\text{FAVI}},\mathcal{L}_{\text{ELBO}},$
and $\mathcal{L}_{\text{IWBO}}$, as in Section 4.2. Table 2 shows that the
ELBO and IWBO tend to be overly concentrated, failing to contain the entire
parameter vector $\theta$ in the $1-\alpha$ highest-density region often.
FAVI, on the other hand, is reasonably well-calibrated. After applying CANVI,
all three methods achieve nearly perfect calibration across a range of desired
confidence levels. Of course, the utility of these corrected regions depends
on the level of information contained in the original model. In the case of
the ELBO or IWBO variational approximations, the corrected regions achieve
statistical validity, but are likely too large to be informative. For FAVI, on
the other hand, application of CANVI results in statistical guarantees on
coverage rates with minimal alterations to the high-density regions.

Table 2: Coverage rates and standard errors for 11-dimensional parameter
$\theta$ before (left) and after (right) conformalization by CANVI, assessed
by inclusion of $\theta$ in the $1-\alpha$ highest density region. Non-
conformalized regions were estimated empirically from batches of 1000 i.i.d.
samples per point.

$1-\alpha$ | ELBO | IWBO | FAVI
---|---|---|---
0.50 | 0.0007 (0.0008) | 0.1031 (0.0110) | 0.5514 (0.0127)
0.75 | 0.0044 (0.0018) | 0.1994 (0.0115) | 0.7534 (0.0127)
0.90 | 0.0195 (0.0044) | 0.3263 (0.0122) | 0.8797 (0.0065)
0.95 | 0.0396 (0.0061) | 0.4074 (0.0186) | 0.9260 (0.0083)
$1-\alpha$ | ELBO | IWBO | FAVI
---|---|---|---
0.50 | 0.4970 (0.0124) | 0.5019 (0.0167) | 0.5086 (0.0144)
0.75 | 0.7488 (0.0139) | 0.7559 (0.0126) | 0.7565 (0.0092)
0.90 | 0.8978 (0.0080) | 0.9036 (0.0111) | 0.9005 (0.0110)
0.95 | 0.9496 (0.0071) | 0.9548 (0.0052) | 0.9487 (0.0081)

##  5 Discussion

We have presented CANVI, a novel framework for producing marginally
calibrated, efficient prediction regions from a collection of variational
approximators with minimal computational and implementation overhead. We view
guaranteeing marginal coverage as an important first step amongst many toward
increasing the utility of such inference algorithms for downstream
applications. This work thus suggests many interesting directions for future
work. Of immediate theoretical interest would be identifying sufficient
conditions for guaranteeing that the exact posterior recovers optimal
efficiency of the CANVI procedure. Additionally, CANVI requires a well-
specified forward model $\mathcal{P}(X\mid\Theta)$ to generate valid
prediction regions; while this is a typical assumption, it generally only
holds approximately in practice. Extensions to the misspecified case could
leverage recent work for using conformal prediction under distribution shift,
such as in [28, 29]. Further, conditional coverage is of interest in certain
scientific applications. While attaining conditional conformal guarantees was
shown to be impossible in the general case in [30], this presumes the lack of
any probabilistic model. We are actively pursuing the extension of CANVI
towards a variant of group-balanced conformal prediction, where the structure
imposed by $q(\Theta\mid X)$ can be exploited to get stronger guarantees.
Finally, leveraging conformal prediction over functional spaces may enable
progress toward achieving guarantees for recovering full posterior
distributions rather than solely for prediction regions.

## References

  * [1] David M Blei, Alp Kucukelbir, and Jon D McAuliffe.  Variational inference: A review for statisticians.  Journal of the American statistical Association, 112(518):859–877, 2017.
  * [2] Kevin P Murphy.  Probabilistic machine learning: an introduction.  MIT press, 2022.
  * [3] Cheng Zhang, Judith Bütepage, Hedvig Kjellström, and Stephan Mandt.  Advances in variational inference.  IEEE transactions on pattern analysis and machine intelligence, 41(8):2008–2026, 2018.
  * [4] Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman.  Yes, but did it work?: Evaluating variational inference.  In International Conference on Machine Learning, pages 5581–5590. PMLR, 2018.
  * [5] Joeri Hermans, Arnaud Delaunoy, François Rozet, Antoine Wehenkel, and Gilles Louppe.  Averting a crisis in simulation-based inference.  arXiv preprint arXiv:2110.06581, 2021.
  * [6] Arnaud Delaunoy, Joeri Hermans, François Rozet, Antoine Wehenkel, and Gilles Louppe.  Towards reliable simulation-based inference with balanced neural ratio estimation.  arXiv preprint arXiv:2208.13624, 2022.
  * [7] Daniel Gilman, Jo Bovy, Tommaso Treu, Anna Nierenberg, Simon Birrer, Andrew Benson, and Omid Sameie.  Strong lensing signatures of self-interacting dark matter in low-mass haloes.  Monthly Notices of the Royal Astronomical Society, 507(2):2432–2447, 2021.
  * [8] Yashar D Hezaveh, Neal Dalal, Daniel P Marrone, Yao-Yuan Mao, Warren Morningstar, Di Wen, Roger D Blandford, John E Carlstrom, Christopher D Fassnacht, Gilbert P Holder, et al.  Detection of lensing substructure using ALMA observations of the dusty galaxy SDP.81.  The Astrophysical Journal, 823(1):37, 2016.
  * [9] S Vegetti, LVE Koopmans, A Bolton, T Treu, and R Gavazzi.  Detection of a dark substructure through gravitational imaging.  Monthly Notices of the Royal Astronomical Society, 408(4):1969–1981, 2010.
  * [10] Simona Vegetti and LVE Koopmans.  Statistics of mass substructure from strong gravitational lensing: quantifying the mass fraction and mass function.  Monthly Notices of the Royal Astronomical Society, 400(3):1583–1592, 2009.
  * [11] David W Hogg and RD Blandford.  The gravitational lens system B1422+231: dark matter, superluminal expansion and the Hubble constant.  Monthly Notices of the Royal Astronomical Society, 268(4):889–893, 1994.
  * [12] Joeri Hermans, Nilanjan Banik, Christoph Weniger, Gianfranco Bertone, and Gilles Louppe.  Towards constraining warm dark matter with stellar streams through neural simulation-based inference.  Monthly Notices of the Royal Astronomical Society, 507(2):1999–2011, 2021.
  * [13] Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  In International conference on machine learning, pages 1530–1538. PMLR, 2015.
  * [14] Abhinav Agrawal, Daniel R Sheldon, and Justin Domke.  Advances in black-box VI: Normalizing flows, importance weighting, and optimization.  Advances in Neural Information Processing Systems, 33:17358–17369, 2020.
  * [15] Kevin P. Murphy.  Probabilistic Machine Learning: Advanced Topics.  MIT Press, 2023.
  * [16] Luca Ambrogioni, Umut Güçlü, Julia Berezutskaya, Eva Borne, Yaǧmur Güçlütürk, Max Hinne, Eric Maris, and Marcel Gerven.  Forward amortized inference for likelihood-free variational marginalization.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 777–786. PMLR, 2019.
  * [17] Jörg Bornschein and Yoshua Bengio.  Reweighted wake-sleep.  arXiv preprint arXiv:1406.2751, 2014.
  * [18] Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov.  Importance weighted autoencoders.  arXiv preprint arXiv:1509.00519, 2015.
  * [19] Anastasios N Angelopoulos and Stephen Bates.  A gentle introduction to conformal prediction and distribution-free uncertainty quantification.  arXiv preprint arXiv:2107.07511, 2021.
  * [20] Glenn Shafer and Vladimir Vovk.  A tutorial on conformal prediction.  Journal of Machine Learning Research, 9(3), 2008.
  * [21] Yachong Yang and Arun Kumar Kuchibhotla.  Finite-sample efficient conformal prediction.  arXiv preprint arXiv:2104.13871, 2021.
  * [22] Matteo Sesia and Emmanuel J Candès.  A comparison of some conformal quantile regression methods.  Stat, 9(1):e261, 2020.
  * [23] Jing Lei, James Robins, and Larry Wasserman.  Distribution-free prediction sets.  Journal of the American Statistical Association, 108(501):278–287, 2013.
  * [24] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows.  Advances in neural information processing systems, 32, 2019.
  * [25] Owen Thomas, Ritabrata Dutta, Jukka Corander, Samuel Kaski, and Michael U. Gutmann.  Likelihood-Free Inference by Ratio Estimation.  Bayesian Analysis, 17(1):1 – 31, 2022.
  * [26] Donald G York, J Adelman, John E Anderson Jr, Scott F Anderson, James Annis, Neta A Bahcall, JA Bakken, Robert Barkhouser, Steven Bastian, Eileen Berman, et al.  The Sloan Digital Sky Survey: Technical summary.  The Astronomical Journal, 120(3):1579, 2000.
  * [27] B Abareshi, J Aguilar, S Ahlen, Shadab Alam, David M Alexander, R Alfarsy, L Allen, C Allende Prieto, O Alves, J Ameel, et al.  Overview of the instrumentation for the Dark Energy Spectroscopic Instrument.  The Astronomical Journal, 164(5):207, 2022.
  * [28] Ryan J Tibshirani, Rina Foygel Barber, Emmanuel Candes, and Aaditya Ramdas.  Conformal prediction under covariate shift.  Advances in neural information processing systems, 32, 2019.
  * [29] Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani.  Conformal prediction beyond exchangeability.  arXiv preprint arXiv:2202.13415, 2022.
  * [30] Vladimir Vovk.  Conditional validity of inductive conformal predictors.  In Asian conference on machine learning, pages 475–490. PMLR, 2012\.
  * [31] Chirag Gupta, Arun K Kuchibhotla, and Aaditya Ramdas.  Nested conformal prediction and quantile out-of-bag ensemble methods.  Pattern Recognition, 127:108496, 2022.
  * [32] Jan-Matthis Lueckmann, Jan Boelts, David Greenberg, Pedro Goncalves, and Jakob Macke.  Benchmarking simulation-based inference.  In International Conference on Artificial Intelligence and Statistics, pages 343–351. PMLR, 2021.
  * [33] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  Advances in neural information processing systems, 32, 2019.
  * [34] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  nflows: normalizing flows in PyTorch, November 2020.
  * [35] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [36] ChangHoon Hahn, KJ Kwon, Rita Tojeiro, Malgorzata Siudek, Rebecca EA Canning, Mar Mezcua, Jeremy L Tinker, David Brooks, Peter Doel, Kevin Fanning, et al.  The DESI PRObabilistic Value-added Bright Galaxy Survey (PROVABGS) mock challenge.  The Astrophysical Journal, 945(1):16, 2023.

##  Appendix A Conformal Score Transforms

######  Lemma A.1.

Let $f:\mathbb{R}\rightarrow\mathbb{R}$ be a strictly monotonically increasing
function,
$\mathcal{D_{C}}=\\{(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}\overset{\mathrm{iid}}{\sim}\mathcal{P}(X,\Theta),\mathcal{S}=\\{s(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}$,
and
$\mathcal{S}^{(f)}=\\{f\left(s(x_{i}^{c},\theta_{i}^{c})\right)\\}_{i=1}^{N_{\mathcal{C}}}$.
For a target coverage $1-\alpha$, denote the
$\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil/N_{\mathcal{C}}$ quantile of
$\mathcal{S}$ as $\widehat{q}(\alpha)$ and that of $\mathcal{S}^{(f)}$ as
$\widehat{q}^{(f)}(\alpha)$. Then, for any $x$,

$\mathcal{C}(x)=\left\\{\theta:s(x,\theta)\leq\widehat{q}(\alpha)\right\\}=\left\\{\theta:f(s(x,\theta))\leq\widehat{q}^{(f)}(\alpha)\right\\}=\mathcal{C}^{(f)}(x).$

###### Proof.

Consider the ordering
$\sigma:\\{1,...,N_{\mathcal{C}}\\}\rightarrow\\{1,...,N_{\mathcal{C}}\\}$ of
$\mathcal{S}$, i.e. $\mathcal{S}_{\sigma(1)}$ is the smallest element of
$\mathcal{S}$ and $\mathcal{S}_{\sigma(N_{\mathcal{C}})}$ the largest. Then,
$\widehat{q}(\alpha)=\mathcal{S}_{\sigma\left(\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil\right)}$.
Notice further that the ordering $\sigma$ is invariant to monotonically
increasing transforms of $s\in\mathcal{S}$. Thus, similarly defining an
ordering $\sigma^{(f)}$ of $\mathcal{S}^{(f)}$, we have $\sigma^{(f)}=\sigma$,
implying

$\widehat{q}^{(f)}(\alpha)=\mathcal{S}^{(f)}_{\sigma^{(f)}\left(\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil\right)}=\mathcal{S}^{(f)}_{\sigma\left(\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil\right)}=f\left(\mathcal{S}_{\sigma\left(\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil\right)}\right)=f(\widehat{q}(\alpha)).$

Further, any strictly monotone $f$ is invertible, from which the equality of
$\mathcal{C}(x)$ and $\mathcal{C}^{(f)}(x)$ follows as

$\displaystyle\mathcal{C}^{(f)}(x)=\left\\{\theta:f(s(x,\theta))\leq\widehat{q}^{(f)}(\alpha)\right\\}=\left\\{\theta:f(s(x,\theta))\leq\widehat{f}(\widehat{q}(\alpha))\right\\}$
$\displaystyle=\left\\{\theta:s(x,\theta)\leq\widehat{q}(\alpha)\right\\}=\mathcal{C}(x),$

completing the proof as desired. ∎

######  Corollary A.2.

Let

$s^{(z)}(x_{i},\theta_{i})=\left|\frac{\mu_{\varphi}(x_{i})-\theta_{i}}{\sqrt{\sigma^{2}_{\varphi}(x_{i})}}\right|\qquad
s^{(p)}(x_{i},\theta_{i})=\frac{1}{\mathcal{N}(\theta_{i};\mu_{\varphi}(x_{i}),\sigma^{2}_{\varphi}(x_{i}))}.$

Let
$\mathcal{D_{C}}=\\{(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}\overset{\mathrm{iid}}{\sim}\mathcal{P}(X,\Theta),\mathcal{S}^{(z)}=\\{s^{(z)}(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}$,
and
$\mathcal{S}^{(p)}=\\{s^{(p)}(x_{i}^{c},\theta_{i}^{c})\\}_{i=1}^{N_{\mathcal{C}}}$.
For a target coverage $1-\alpha$, denote the
$\lceil(N_{\mathcal{C}}+1)(1-\alpha)\rceil/N_{\mathcal{C}}$ quantile of
$\mathcal{S}^{(p)}$ as $\widehat{q}^{(p)}(\alpha)$ and that of
$\mathcal{S}^{(z)}$ as $\widehat{q}^{(z)}(\alpha)$. Then, for any $x$,

$\mathcal{C}^{(p)}(x)=\left\\{\theta:s^{(p)}(x,\theta)\leq\widehat{q}^{(p)}(\alpha)\right\\}=\left\\{\theta:s^{(z)}(x,\theta)\leq\widehat{q}^{(z)}(\alpha)\right\\}=\mathcal{C}^{(z)}(x).$

###### Proof.

Notice there is a strictly monotonically increasing map between
$s^{(z)}(x,\theta)$ and $s^{(p)}(x,\theta)$, given by $f(s):=\phi(s)$, where
$\phi(\cdot)$ denotes the standard normal density. That is:

$s^{(p)}(x,\theta)=\phi(s^{(z)}(x,\theta)).$

The equivalence between $\mathcal{C}^{(p)}(x)$ and $\mathcal{C}^{(z)}(x)$ then
follows from the preceding lemma. ∎

##  Appendix B Predictive Efficiency of CANVI (Detailed)

We now restate the assumptions from the main body, present a proof sketch for
the theorem, and finally present the proof.

######  Assumption 3.

For each $t$, $\exists$ $r^{*},\gamma\in(0,1)$ such that the inverse CDF of
the score under $q^{(t)}$, $F_{t}^{-1}$, is $\gamma$-Hölder continuous on
[$\widehat{q}^{(t)}(\alpha)-r^{*},\widehat{q}^{(t)}(\alpha)+r^{*}$] with
Hölder continuity constant $L_{t}$.

######  Assumption 4.

For each $t$, the map $\alpha\rightarrow\widehat{\ell}(q^{(t)},\alpha)$ is
Lipschitz continuous with constant $L_{W}$.

The proof emerges through a reduction of CANVI to a special case of the
“Validity First Conformal Prediction for the Smallest Prediction Set” (VFCP)
algorithm presented in the [21], which was proven to satisfy Theorem 3.2. The
VFCP algorithm is provided for convenience in Algorithm 2. Note that we made
the appropriate replacements in notations to match those presented in the main
body of our paper to clarify the connection to CANVI.

To further underscore the parallel, we present an immaterially modified
version of VFCP, where it is assumed the data are split prior to the algorithm
execution into disjoint sets $\mathcal{D},\mathcal{D}_{\mathcal{C}}$, and
$\mathcal{D}_{\mathcal{R}}$. We assume input prediction methods
$\mathcal{A}^{(1)},...,\mathcal{A}^{(T)}$ have been trained on $\mathcal{D}$.
These two modifications in the algorithm presentation have no manifest effects
in the proof of [21], meaning all results as presented there apply with the
appropriate choices of inputs in Algorithm 2. We further make explicit that
VFCP is formulated agnostic to the functional choice of inverse efficiency
estimator $\text{Width}(\cdot)$.

We similarly consider an immaterially modified form of CANVI, that is, one in
which $\mathcal{D}_{\mathcal{C}}$ and $\mathcal{D}_{\mathcal{R}}$ are pre-
generated in precisely the fashion outlined Algorithm 1 and provided as input.
Thus, inputs to CANVI parallel those of VFCP, namely $\\{q^{(t)}(\Theta\mid
X)\\}_{t=1}^{T},1-\alpha,D_{\mathcal{C}},$ and $D_{\mathcal{R}}$.

Algorithm 2 Validity First Conformal Prediction (VFCP) [21]

1:Procedure: VFCP

2:Inputs:  Trained methods $\\{\mathcal{A}^{(t)}\\}_{t=1}^{T}$, Desired
coverage $1-\alpha$, Calibration set $D_{\mathcal{C}}$, Recalibration set
$D_{\mathcal{R}}$, Inverse efficiency estimator
$\text{Width}:\mathcal{U}\subset\mathbb{R}^{D}\rightarrow\mathbb{R}$

3:Using $\mathcal{A}^{(t)}$, construct an increasing (nested) sequence of sets
$\\{F^{(t)}_{y}\\}_{y\in\mathcal{Y}}$, where $\mathcal{Y}\subset\mathbb{R}$.

4:Compute conformal prediction set $\mathcal{C}^{(t)}$ based on
$\\{F^{(t)}_{y}\\}_{y\in\mathcal{Y}}$. Specifically, for each
$(x_{i},\theta_{i})\in\mathcal{D}_{\mathcal{C}}$ and $t\in\\{1,2,...,T\\}$,
define the conformal score
$s^{(t)}(x_{i},\theta_{i}):=\inf_{y\in\mathcal{Y}}\\{(x_{i},\theta_{i})\in
F^{(t)}_{y}\\}.$

5:Compute the corresponding conformal prediction set as
$\mathcal{C}^{(t)}:=\\{(x,\theta):s^{(t)}(x,\theta)\leq\widehat{q}^{(t)}(D_{\mathcal{C}},\alpha)\\},$

6:where $\widehat{q}^{(t)}(D_{\mathcal{C}},\alpha)$ is the
$\lceil(|\mathcal{D}_{\mathcal{C}}|+1)(1-\alpha)\rceil$-th largest element of
$\\{s^{(t)}(x_{i},\theta_{i})\\}_{i\in\mathcal{D}_{\mathcal{C}}}$.

7:Set  $t^{(*)}:=\operatorname*{arg\,min}_{1\leq t\leq
T}\text{Width}(\mathcal{C}^{(t)})$

8:For each $(x_{i},\theta_{i})\in\mathcal{D}_{\mathcal{R}}$, define the
conformal score
$s_{\mathcal{R}}^{(*)}(x_{i},\theta_{i}):=\inf_{y\in\mathcal{Y}}\\{(x_{i},\theta_{i})\in
F^{(t^{(*)})}_{y}\\}.$

9:Compute the corresponding conformal prediction set as
$\mathcal{C}^{(*)}:=\\{(x,\theta):s_{\mathcal{R}}^{(*)}(x,\theta)\leq\widehat{q}^{(*)}(D_{\mathcal{R}},\alpha)\\},$

10:where
$\widehat{q}^{(*)}(D_{\mathcal{R}},\alpha):=\lceil(|\mathcal{D}_{\mathcal{R}}|+1)(1-\alpha)\rceil$-th
largest element of
$\\{s_{\mathcal{R}}^{(*)}(x_{i},\theta_{i})\\}_{i\in\mathcal{D}_{\mathcal{R}}}$.

11:Return the prediction set $\mathcal{C}^{(*)}$.

Intuitively, the proof follows from the fact that CANVI is a special case of
Algorithm 2. Algorithm 2 is presented using the so-called “nested set
formulation” of conformal prediction. Briefly, the nested sets formulation
starts by requiring the user to specify a set-valued function $F_{r}$ instead
of a score function. Sets must be nested over increasing $r$. Despite being
seemingly more expressive, the equivalence of the nested set and split score
conformal methods was demonstrated in [31]. In particular, for a chosen score
function $s(x,\theta)$, the corresponding nested set function would be

$F_{\widehat{q}(\alpha)}=\\{\theta\mid s(x,\theta)\leq\widehat{q}(\alpha)\\},$

where $r:=\widehat{q}(\alpha)$. The presentation of Algorithm 2 using the
nested set formulation was, thus, a means to allow for the employment of
relevant proof strategies from classical learning theory. We now proceed
through the formal equivalence of CANVI and a special case of VFCP.

######  Theorem B.1.

Let $q^{(*)}(\Theta\mid
X),\widehat{q}^{(*)}=\mathrm{CANVI}\left(\\{q^{(t)}(\Theta\mid
X)\\}_{t=1}^{T},1-\alpha,D_{\mathcal{C}},D_{\mathcal{R}}\right).$ Let
$\mathcal{D},\mathcal{D_{C}},\mathcal{D_{T}}$, and $\mathcal{D_{R}}$ be drawn
i.i.d. from $\mathcal{P}(\Theta,X)$, where $\mathcal{D}$ is the aggregate of
training data used for $\\{q^{(t)}(\Theta\mid X)\\}_{t=1}^{T}$. Take
$\delta\in[0,1]$. If Assumption 2 holds and Assumption 1 holds with

$r^{*}\geq\max\left\\{\sqrt{\frac{\log(4T/\delta)}{2N_{\mathcal{C}}}},\frac{2}{N_{\mathcal{C}}}\right\\},$
(12)

then with probability at least $1-\delta$,

$\widehat{\ell}_{\mathcal{R}}(q^{(*)},\alpha)\leq\min_{1\leq t\leq
T}\widehat{\ell}(q^{(t)},\alpha)+3L_{W}L_{[T]}\left[\left(\frac{\log(4T/\delta)}{N_{\mathcal{C}}}\right)^{\gamma/2}+\left(\frac{2}{N_{\mathcal{C}}}\right)^{\gamma}\right],$
(13)

where $\gamma$, $L_{W}$, and $L_{[T]}=\max_{1\leq t\leq T}L_{t}$ are constants
defined in Assumptions 1 and 2.

###### Proof.

To demonstrate the reduction, it suffices to illustrate

$\displaystyle\mathcal{C}^{\textsc{CANVI}}:=\left\\{(x,\theta):1/q^{(*)}(\theta\mid
x)\leq\widehat{q}^{(*)}(\alpha)\right\\}$
$\displaystyle=\mathrm{VFCP}\left(\\{q^{(t)}(\Theta\mid
X)\\}_{t=1}^{T},1-\alpha,D_{\mathcal{C}},D_{\mathcal{R}},\widehat{\ell}\right):=\mathcal{C}^{\textsc{VFCP}},\text{
where
}\widehat{\ell}(\mathcal{C}):=\frac{1}{N_{\mathcal{T}}}\sum_{i=1}^{N_{\mathcal{T}}}\mathcal{L}\left(\mathcal{C}(x_{i})\right).$

Importantly, note that $\mathcal{C}^{\textsc{CANVI}}$ is a set over
$\text{Dom}(X)\times\text{Dom}(\Theta)$, unlike
$\mathcal{C}^{\textsc{CANVI}}(x)$ from the main text, which is a set-valued
mapping from $X\rightarrow\vartheta\subset\text{Dom}(\Theta)$. However, since
$\widehat{q}^{(*)}$ is fixed over $x$,

$\mathcal{C}^{\textsc{CANVI}}=\bigcup_{x\in\text{Dom}(X)}\\{x\\}\times\mathcal{C}^{\textsc{CANVI}}(x).$

It is noted in [21] that an appropriate $\text{Width}(\mathcal{C})$ for
prediction tasks with unbounded $\text{Dom}(X)$ is “the average Lebesgue
measure of the x-cross-section of the prediction set,” as is estimated by
$\widehat{\ell}$. We further highlight that the choice of
$\text{Width}(\mathcal{C})=\widehat{\ell}(\mathcal{C})$ is random through the
choice of $\mathcal{D}_{\mathcal{T}}$; however, such randomness is explicitly
noted in [21] as being a supported property of the chosen
$\text{Width}(\mathcal{C})$.

The equivalence of $\mathcal{C}^{\textsc{CANVI}}$ and
$\mathcal{C}^{\textsc{VFCP}}$, therefore, can be shown in demonstrating the
equivalence in corresponding scores, as the resulting empirical score
distributions and hence quantiles over $\mathcal{D}_{\mathcal{C}}$ and
$\mathcal{D}_{\mathcal{R}}$, i.e. $\widehat{q}^{(t)}(\alpha)$ and ultimately
$\widehat{q}^{(*)}(\alpha)$, and future prediction regions follow to then be
equivalent. We employ the notation from Section 3.3, where
$\widehat{q}^{(t)}(\alpha)$ is the exact $1-\alpha$ quantile of
$1/q^{(t)}(\Theta\mid X)$. Then, we can define nested sets as described
previously, namely choosing

$\mathcal{Y}:=\\{\widehat{q}^{(t)}(\alpha)\\}_{\alpha\in[0,1]}\quad\text{ and
}\quad F^{(t)}_{y}:=\left\\{(x,\theta):1/q^{(t)}(\theta\mid x)\leq
y\right\\}.$

Notice that this choice of $\\{F^{(t)}_{y}\\}_{y\in\mathcal{Y}}$ recovers the
original score formulation, as:

$\displaystyle
s^{(t)}(x,\theta)=\inf_{y\in\mathcal{Y}}\left\\{(x,\theta)\in\left\\{(x^{\prime},\theta^{\prime}):1/q^{(t)}(\theta^{\prime}\mid
x^{\prime})\leq y\right\\}\right\\}$
$\displaystyle=\inf_{y\in\mathcal{Y}}\left\\{1/q^{(t)}(\theta\mid x)\leq
y\right\\}=1/q^{(t)}(\theta\mid x).$

The equivalence of $\mathcal{C}^{\textsc{CANVI}}$ and
$\mathcal{C}^{\textsc{VFCP}}$, thus, follows as previously outlined as a
consequence of the equivalence of the corresponding score functions. ∎

##  Appendix C Fixed Coverage Posterior Counterexample

As noted in the main text, we now demonstrate that, for any particular
$\alpha$, the true posterior does not always lead to the greatest efficiency.
Formally, we wish to disprove the following conjecture:

######  Conjecture 2.

Let $(X,\Theta)\sim\mathcal{P}(\Theta,X)$, $q^{(*)}(\Theta\mid
X)=\mathcal{P}(\Theta\mid X)$ and $q(\Theta\mid X)$ denote any other posterior
approximator. Then for any $\alpha\in[0,1]$,

$\ell(q^{(*)},\alpha)\leq\ell(q,\alpha).$ (14)

###### Proof.

Consider the distribution on $X,\Theta$ given by

$\displaystyle\mathcal{P}(X)\sim\text{Bernoulli}(x;1/2)$
$\displaystyle\mathcal{P}(\Theta\mid
x)\sim\begin{cases}\text{Uniform}(\theta;[0,200])&\text{if }x=0\\\
\text{Uniform}(\theta;[200,300])&\text{if }x=1\end{cases}$

We want to make prediction intervals for $\Theta\mid X$. Using the true
posterior distribution, we first wish to find the distribution of
$s(X,\Theta)=1/\mathcal{P}(\Theta\mid X)$. Notice that, for
$X=x,\Theta=\theta$,

$\displaystyle s(x,\theta)=\begin{cases}200&\text{if }x=0\\\ 100&\text{if
}x=1\end{cases}\implies F^{-1}(\alpha)=\begin{cases}100&\text{if
}\alpha<0.5\\\ 200&\text{if }\alpha\geq 0.5\end{cases},$

where $F^{-1}(\alpha)$ is the inverse CDF of $s(X,\Theta)$.
$\widehat{q}(\alpha)$ can now be obtained from $F^{-1}(\alpha)$. Consider now
an $\alpha<0.5$. $\widehat{q}(\alpha)=100$ in this case, meaning prediction
regions will take the form
$\mathcal{C}(x)=\left\\{\theta:1/100\leq\mathcal{P}(\theta\mid x)\right\\}$.
To compute $\ell(\mathcal{P},\alpha)$, we enumerate over the values of $X$,
using that $\forall$ $\theta$, $\mathcal{P}(\theta\mid X=0)=1/200$ and
$\mathcal{P}(\theta\mid X=1)=1/100$:

$\displaystyle\mathcal{C}(x)=\begin{cases}\left\\{\theta:1/100\leq
1/200\right\\}=\emptyset&\text{if }x=0\\\ \left\\{\theta:1/100\leq
1/100\right\\}=[200,300]&\text{if }x=1\end{cases}.$

Computing $\ell(\mathcal{P},\alpha)$ immediately follows as

$\displaystyle\ell(\mathcal{P},\alpha)=\mathbb{E}_{X}\left[\mathcal{L}\left(\mathcal{C}(X)\right)\right]=\mathcal{P}(X=0)\mathcal{L}\left(\mathcal{C}(0)\right)+\mathcal{P}(X=1)\mathcal{L}\left(\mathcal{C}(1)\right)=\frac{1}{2}\left(0+100\right)=50.$

Suppose now we have an incorrect posterior of the form

$\displaystyle q_{b}(\Theta\mid
x)\sim\begin{cases}\text{Uniform}(\theta;[0,b])&\text{if }x=0\\\
\text{Uniform}(\theta;[200,300])&\text{if }x=1\end{cases},$

where $b<100$. We now again find the distribution of the score. Importantly,
remember that the underlying joint distribution of $(X,\Theta)$ has not
changed. Therefore, for $X=x,\Theta=\theta$,

$\displaystyle s(x,\theta)=\begin{cases}b&\text{if }x=0,\theta\leq b\\\
100&\text{if }x=1\\\ \infty&\text{if }x=0,\theta>b\end{cases}\implies
F^{-1}(\alpha)=\begin{cases}b&\text{if }\alpha<b/200\\\ 100&\text{if
}b/200\leq\alpha<b/200+1/2\\\ \infty&\text{if }\alpha\geq
b/200+1/2\end{cases},$

where the difference arises from draws $x=0,\theta>b$ from the true
distribution where $q_{b}(\Theta\mid x)$ ascribes no mass. For the case of
$\alpha<b/200$, we again compute the average length as above:

$\displaystyle\mathcal{C}(x)=\begin{cases}\left\\{\theta:1/b\leq
1/b\right\\}=[0,b]&\text{if }x=0\\\ \left\\{\theta:1/b\leq
1/100\right\\}=\emptyset&\text{if
}x=1\end{cases}\implies\ell(q_{b},\alpha)=\frac{b}{2}<50.$

Thus, $\ell(q_{b},b/200-\epsilon)<\ell(\mathcal{P},b/200-\epsilon)$,
completing the counterexample. ∎

##  Appendix D Gaussian CANVI Efficiency (Detailed)

######  Theorem D.1.

Let $\Theta$ and $X$ be zero-mean unit-variance Gaussian random variables with
correlation $\rho$. Let $q_{\varphi}(\theta|x)=\mathcal{N}(\theta;\varphi
x,1-\rho^{2})$. Then $q_{\rho}$ gives the true posterior and
$\ell(q_{\rho},\alpha)=\min_{\varphi}\ell(q_{\varphi},\alpha)$.

###### Proof.

Notice that in the bivariate Gaussian case, we have closed forms for the
following:

$\displaystyle\Theta\mid x\sim\mathcal{N}(\rho x,1-\rho^{2})\qquad
X\mid\theta\sim\mathcal{N}(\rho\theta,1-\rho^{2})$
$\displaystyle\Theta\sim\mathcal{N}(0,1)\qquad X\sim\mathcal{N}(0,1).$

We wish to find the distribution of $s(X,\Theta)=1/q(\Theta\mid X)$ jointly
over $X,\Theta$ to find $\widehat{q}_{\varphi}$ explicitly. The CDF of this
score can be computed as follows:

$\displaystyle\mathcal{P}\left(1/q_{\varphi}(\Theta\mid X)\leq
q\right)=\mathcal{P}\left(\sqrt{2\pi(1-\rho^{2})}e^{\frac{(\varphi
X-\Theta)^{2}}{2(1-\rho^{2})}}\leq q\right)$
$\displaystyle=\mathcal{P}\left(R^{2}\leq
2(1-\rho^{2})\log\left(\sqrt{\frac{q^{2}}{2\pi(1-\rho^{2})}}\right)\right)=\mathcal{P}\left(R\leq\sqrt{(1-\rho^{2})\log\left(\frac{q^{2}}{2\pi(1-\rho^{2})}\right)}\right),$

where $R:=|\varphi X-\Theta|$. From the above calculation,
$\widehat{q}_{\varphi}$ must satisfy:

$\mathcal{P}\left(R\leq\sqrt{(1-\rho^{2})\log\left(\frac{\widehat{q}_{\varphi}^{2}}{2\pi(1-\rho^{2})}\right)}\right)=1-\alpha.$

Notice now that, since $(X,\Theta)$ are bivariate Gaussian:

$\displaystyle\varphi
X-\Theta\sim\mathcal{N}(0,\varphi^{2}+1-2\varphi\rho)\implies
R\sim\text{HalfNormal}(\varphi^{2}+1-2\varphi\rho).$

Therefore, the $1-\alpha$ quantile of $R$ is
$\sqrt{\varphi^{2}+1-2\varphi\rho}\Phi^{-1}\left(1-\frac{\alpha}{2}\right)$.
Solving for $\widehat{q}_{\varphi}$ in this quantile produces the final
threshold:

$\displaystyle\sqrt{(1-\rho^{2})\log\left(\frac{\widehat{q}_{\varphi}^{2}}{2\pi(1-\rho^{2})}\right)}=\sqrt{\varphi^{2}+1-2\varphi\rho}\Phi^{-1}\left(1-\frac{\alpha}{2}\right)$
$\displaystyle\implies\log\left(\frac{\widehat{q}_{\varphi}^{2}}{2\pi(1-\rho^{2})}\right)=\frac{\varphi^{2}+1-2\varphi\rho}{1-\rho^{2}}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right)^{2}$
$\displaystyle\implies\widehat{q}_{\varphi}=\sqrt{2\pi(1-\rho^{2})\exp\left(\frac{\varphi^{2}+1-2\varphi\rho}{1-\rho^{2}}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right)^{2}\right)}.$

We now wish to see the expected confidence interval size with this score
threshold and demonstrate that this is smaller under closer variational
posterior estimates, i.e. with $\varphi=\rho$. Consider, therefore, an
arbitrary fixed $X=x$. We wish to find explicitly
$\mathcal{L}(\mathcal{C}(x))$. Notice that this immediately follows from the
calculation we performed above, namely:

$\displaystyle\mathcal{C}(x)=\\{\theta\mid 1/q_{\varphi}(\theta\mid
x)\leq\widehat{q}_{\varphi}\\}=\left\\{\theta\mid(\varphi x-\theta)^{2}\leq
2(1-\rho^{2})\log\left(\sqrt{\frac{\widehat{q}_{\varphi}^{2}}{2\pi(1-\rho^{2})}}\right)\right\\}$
$\displaystyle=\left\\{\theta\in\varphi
x\pm\sqrt{(1-\rho^{2})\log\left(\frac{\widehat{q}_{\varphi}^{2}}{2\pi(1-\rho^{2})}\right)}\right\\}.$

Plugging in the explicit form of $\widehat{q}_{\varphi}$, we end up with:

$\displaystyle\mathcal{L}(\mathcal{C}(x))=2\sqrt{(1-\rho^{2})\log\left(\frac{1}{2\pi(1-\rho^{2})}2\pi(1-\rho^{2})\exp\left(\frac{\varphi^{2}+1-2\varphi\rho}{1-\rho^{2}}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right)^{2}\right)\right)}$
$\displaystyle=2\sqrt{\varphi^{2}+1-2\varphi\rho}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right).$

Critically, $\mathcal{L}(\mathcal{C}(x))$ is constant across $x$. Therefore,

$\ell(q_{\varphi},\alpha)=\mathbb{E}_{X}\left[\mathcal{L}\left(\mathcal{C}(X)\right)\right]=2\sqrt{\varphi^{2}+1-2\varphi\rho}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right).$
(15)

Now observe that length is minimized with

$\displaystyle\frac{d\ell}{d\varphi}=0\implies\frac{d}{d\varphi}\left(2\sqrt{\varphi^{2}+1-2\varphi\rho}\left(\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\right)\right)$
$\displaystyle=C\left(\frac{\varphi-\rho}{\sqrt{\varphi^{2}+1-2\varphi\rho}}\right)=0\implies\varphi=\rho$

Notice that $\frac{d\ell}{d\varphi}<0$ for $\varphi<\rho$ and
$\frac{d\ell}{d\varphi}>0$ for $\varphi>\rho$, which demonstrates set size is
an increasing function of $|\varphi-\rho|$ with a minimum achieved at the true
variational posterior of $\varphi=\rho$. ∎

###  D.1 Empirical Validation

We validate the result obtained in Equation 15. To do so, we assume a setup
of:

$X,\Theta\sim\mathcal{N}\left(\begin{pmatrix}10\\\
-3\end{pmatrix},\begin{pmatrix}1&0.3\\\ 0.3&1\end{pmatrix}\right)$

A fixed target coverage of $1-\alpha=0.95$ and, thus,
$\Phi^{-1}(1-\alpha/2)=1.96$ was used. Empirical lengths
$\widehat{\ell}(q_{\phi},\alpha)$ were computed at each $\phi\in[-1,1]$ at
discretizations of 0.01.

$\widehat{\ell}(q_{\phi},\alpha)$ estimated for a fixed test batch of 1000
points ($|\mathcal{D_{T}}|=1000$). Due to the light-tailed nature of the
marginal $\Theta$ distribution, we employ an explicit grid-discretization to
compute $\mathcal{L}(\mathcal{C}(x^{\prime}))$ for each
$x^{\prime}\in\mathcal{D_{T}}$, unlike the importance sampled Monte Carlo
estimate employed in the main text. Grid discretization was taken for
$\theta\in[-10,10]$ in steps of $0.1$. This empirical coverage is plotted
against the theoretically prediction of Equation 15 in Figure 4, validating
the above finding.

Figure 4:  The set size curve $\widehat{\ell}(q_{\varphi},\alpha)$ with
$\rho=0.3$ as a function of the variational approximation $\varphi$. The
theoretical prediction from Equation 15 is drawn in orange and corresponding
empirical confidence interval lengths in blue.

##  Appendix E Simulation-Based Inference Benchmarks

The benchmark tasks are a subset of those provided by [32]. For convenience,
we provide brief descriptions of the tasks curated by this library; however, a
more comprehensive description of these tasks can be found in their
manuscript.

###  E.1 Gaussian Linear

10-dimensional Gaussian model with a Gaussian prior:

$\displaystyle\text{{Prior}: }\mathcal{N}(0,0.1\odot I)$
$\displaystyle\text{{Simulator}:
}x\mid\theta\sim\mathcal{N}(x\mid\theta,0.1\odot I)$

###  E.2 Gaussian Linear Uniform

10-dimensional Gaussian model with a uniform prior:

$\displaystyle\text{{Prior}: }\mathcal{U}(-1,1)$
$\displaystyle\text{{Simulator}:
}x\mid\theta\sim\mathcal{N}(x\mid\theta,0.1\odot I)$

###  E.3 SLCP with Distractors

Simple Likelihood Complex Posterior (SLCP) with Distractors has uninformative
dimensions in the observation over the standard SLCP task:

$\displaystyle\text{{Prior}: }\mathcal{U}(-3,3)$
$\displaystyle\text{{Simulator}: }x\mid\theta=p(y)\text{ where }p\text{
reorders }y\text{ with a fixed random order }$ $\displaystyle
y_{[1:8]}\sim\mathcal{N}\left(\begin{bmatrix}\theta_{1}\\\
\theta_{2}\end{bmatrix},\begin{bmatrix}\theta_{3}^{4}&\theta_{3}^{2}\theta_{4}^{2}\tanh(\theta_{5})\\\
\theta_{3}^{2}\theta_{4}^{2}\tanh(\theta_{5})&\theta_{4}^{4}\end{bmatrix}\right),\qquad
y_{9:100}\sim\frac{1}{20}\sum_{i=1}^{20}t_{2}(\mu^{i},\Sigma^{i}),$
$\displaystyle\mu^{i}\sim\mathcal{N}(0,15^{2}I),\Sigma^{i}_{j,k}\sim\mathcal{N}(0,9),\Sigma^{i}_{j,j}=3e^{a},a\sim\mathcal{N}(0,1),$

###  E.4 Bernoulli GLM Raw

10-parameter GLM with Bernoulli observations and Gaussian prior. Observations
are not sufficient statistics, unlike the standard “Bernoulli GLM” task:

$\displaystyle\text{{Prior}:
}\beta\sim\mathcal{N}(0,2),f\sim\mathcal{N}(0,(F^{T}F)^{-1})$
$\displaystyle\qquad
F_{i,i-2}=1,F_{i,i-1}=-2,F_{i,i}=1+\sqrt{\frac{i-1}{9}},F_{i,j}=0;i\leq j$
$\displaystyle\text{{Simulator}:
}x_{i}\mid\theta\sim\text{Bern}(\eta(v_{i}^{T}f+\beta)),\eta(\odot)=\exp(\odot)/(1+\exp(\odot))$

###  E.5 Gaussian Mixture

A mixture of two Gaussians, with one having a much broader covariance
structure:

$\displaystyle\text{{Prior}: }\beta\sim\mathcal{U}(-10,10)$
$\displaystyle\text{{Simulator}: }x\mid\theta\sim
0.5\mathcal{N}(x\mid\theta,I)+0.5\mathcal{N}(x\mid\theta,.01I)$

###  E.6 Two Moons

Task with a posterior that has both global (bimodal) and local (crescent-
shaped) structure:

$\displaystyle\text{{Prior}: }\beta\sim\mathcal{U}(-1,1)$
$\displaystyle\text{{Simulator}:
}x\mid\theta=\begin{bmatrix}r\cos(\alpha)+0.25\\\
r\sin(\alpha)\end{bmatrix}+\begin{bmatrix}-|\theta_{1}+\theta_{2}|/\sqrt{2}\\\
(-\theta_{1}+\theta_{2})/\sqrt{2}\end{bmatrix}$
$\displaystyle\qquad\alpha\sim\mathcal{U}(-\pi/2,\pi/2),r\sim\mathcal{N}(0.1,0.01^{2})$

###  E.7 SIR

Epidemiology model with $S$ (susceptible), $I$ (infected), and $R$
(recovered). A contact rate $\beta$ and mean recovery rate of $\gamma$ are
used as follows:

$\displaystyle\text{{Prior}:
}\beta\sim\text{LogNormal}(\log(0.4),0.5),\gamma\sim\text{LogNormal}(\log(1/8),0.2)$
$\displaystyle\text{{Simulator}:
}x=(x_{i})_{i=1}^{10};x_{i}\mid\theta\sim\text{Bin}(1000,\frac{I}{N}),\text{
where }I\text{ is simulated from: }$
$\displaystyle\frac{dS}{dt}=-\beta\frac{SI}{N},\qquad\frac{dI}{dt}=\beta\frac{SI}{N}-\gamma
I,\qquad\frac{dR}{dt}=\gamma I$

###  E.8 Lotka-Volterra

An ecological model commonly used in describing dynamics of competing species.
$\theta$ parameterizes this interaction as
$\theta=(\alpha,\beta,\gamma,\delta)$:

$\displaystyle\text{{Prior}:
}\alpha\sim\text{LogNormal}(-.125,0.5),\qquad\beta\sim\text{LogNormal}(-3,0.5)$
$\displaystyle\gamma\sim\text{LogNormal}(-.125,0.5),\qquad\delta\sim\text{LogNormal}(-3,0.5)$
$\displaystyle\text{{Simulator}:
}x=(x_{i})_{i=1}^{10};x_{1,i}\mid\theta\sim\text{LogNormal}(\log(X),0.1),x_{2,i}\mid\theta\sim\text{LogNormal}(\log(Y),0.1)$
$\displaystyle\text{ where }X,Y\text{ is simulated from: }\frac{dX}{dt}=\alpha
X-\beta XY,\qquad\frac{dY}{dt}=-\gamma Y+\delta XY$

###  E.9 ARCH

The two-dimensional parameter $\theta=(\theta_{1},\theta_{2})$ includes both
an autoregressive component ($\theta_{1}$) and a component controlling the
level of conditional noise ($\theta_{2}$). Given a full realization of the
time series $y_{1:T}$, we aim to amortize inference over $\theta$. Priors are
taken to be $\theta_{1}\sim\textrm{Unif}(-1,1)$ and
$\theta_{2}\sim\textrm{Unif}(0,1)$. One important change from the model of
[25] is that we fix $e^{(0)}=0$ rather than drawing this quantity from a
standard Gaussian.

The Adam optimizer was used with a learning rate of 0.0001 for 25,000 training
steps for each of the three methods: IWBO ($K=10$), ELBO, FAVI. Due to
constraints arising from the uniform priors on the parameters, the IWBO and
ELBO implementations rely on logit transformations of the latent random
variables $\theta$ to avoid zero-density regions that result in undefined
gradients. The encoder network is trained to learn distributions on the
unconstrained space of the transformed random variable $\theta^{\prime}$, and
visualizations are produced by performing the inverse transformation. While
FAVI avoids these issues because it is likelihood-free, for an apples-to-
apples comparison we also implement FAVI on the unconstrained latent space as
well.

##  Appendix F Training Details

All encoders were implemented in PyTorch [33] with a Neural Spline Flow
architecture. The NSF was built using code from [34]. Specific architecture
hyperparameter choices were taken to be the defaults from [34] and are
available in the code. Optimization was done using Adam [35] with a learning
rate of $10^{-3}$ over 5,000 training steps. Minibatches were drawn from the
corresponding prior $\mathcal{P}(\Theta)$ and simulator
$\mathcal{P}(X\mid\Theta)$ as specified per task in the preceding section.
Training these models required between 10 minutes and two hours using an
Nvidia RTX 2080 Ti GPUs for each of the SBI tasks.

##  Appendix G Prediction Regions

Credible regions for $q_{\varphi}(\Theta\mid x)$ (for a single
$x\sim\mathcal{P}(X)$) on a subset of the SBI benchmark tasks are plotted
below over varying degrees of training, as indicated in the figures. As
expected, the efficiency of the prediction regions improves over training
across all tasks, resulting in smaller regions for each target coverage
$1-\alpha$.

###  G.1 Gaussian Linear

###  G.2 Gaussian Linear Uniform

###  G.3 SLCP

###  G.4 Bernoulli GLM Raw

###  G.5 Gaussian Mixture

###  G.6 Two Moons

##  Appendix H Posteriors

We provide visualizations of approximate and reference posteriors (produced
with MCMC from [32]) to justify the overdispersion claims made on the
variational approximation procedure.

###  H.1 Gaussian Linear

###  H.2 Gaussian Mixture

###  H.3 Gaussian Linear Uniform

###  H.4 Two Moons

###  H.5 SLCP

###  H.6 Bernoulli GLM Raw

###  H.7 ARCH

For selected points from the training sets, we show the exact and approximate
posteriors obtained by training according to the ELBO, IWBO, and FAVI
objectives. The ground-truth parameter value is noted in red.

Depending on the miscalibration of the variational posterior on either per-
objective or per-point basis, the conformalized $1-\alpha$ high-density region
(HDR) either shrinks or expands, and can be visualized in two dimensions.
Figure 5, Figure 6 show examples where the 50% prediction regions obtained
from the FAVI-learned variational posterior shrinks and grows after applying
CANVI.

Figure 5: 50% prediction region, observation 7, before and after applying
CANVI. Figure 6: 50% prediction region, observation 8, before and after
applying CANVI.

We visualize the efficiency of the iterates $q_{\varphi}$ across training
iterations in Figure 7. Recalling the form of Equation 3, it is unsurprising
that FAVI tends to be more efficient, as it trains using simulated data (over
which Equation 3 is computed). To estimate Equation 3, at every 500 training
steps, we simulate 20 $\theta,x$ pairs from the forward model. A larger number
of Monte Carlo samples can be used but at an increased cost. As the resulting
estimates are noisy, we smooth the resulting series with a Savitzky-Golay
filter with a window length of 10 and third-degree polynomial order for better
visualization.

Figure 7: Efficiency estimates for variational posteriors trained by ELBO,
IWBO, and FAVI across training iterations.

##  Appendix I Predictive Efficiency Over Training Objectives

As an idealized “ground truth” computation, we additionally provide the
predictive efficiency as computed on a $100\times 100$ discretized grid of the
support of the prior, $[\theta_{\min},\theta_{\max}]$. Notice that graphs are
understandably less noisy than their counterparts as estimated using
importance-weighted Monte Carlo sampling. The drawback, however, is both the
computational intensity of such sampling and the restriction to low-
dimensional cases.

Figure 8:  Estimation was performed for a fixed $1-\alpha=0.95$, suggesting
efficiency improves as $q^{(t)}(\Theta\mid X)\rightarrow\mathcal{P}(\Theta\mid
X)$. Error bars from estimating $\widehat{\ell}(q_{\varphi},\alpha)$ across
test batches are plotted.

##  Appendix J SED Experimental Details

The PROVABGS emulator (Section 4.3) was trained to minimize the MSE using
normalized simulated PROVABGS outputs with fixed log stellar mass parameter
[36]. Training data were generated from PROVABGS with a fixed magnitude
parameter ($\theta_{0}=10.5$), resampled onto a 5 Angstrom grid, and
normalized to integrate to one. After training, forward passes through the
emulator are significantly faster than the base simulator. We provide two
simulated draws from our emulator in Figure 9. We use the recommended priors
from [36] on the remaining eleven parameters. As these are highly constrained
(uniform priors, vastly different scales, and a 4-dimensional vector on the
simplex), we similarly operate on an unconstrained, 10-dimensional space by
invertible transformations. All three methods were trained for 10,000 steps
using the Adam optimizer with learning rate 0.0001. $K=1000$ was used for the
IWBO.

Figure 9: Two example draws from our neural network emulator of PROVABGS.

We let the output of our simulator be a mean parameter
$\mu\in\mathbb{R}^{1400}$, and generate observed data independently binwise
via $x_{i}\sim\mathcal{N}(\mu_{i},|\mu_{i}|^{2}\sigma^{2})$ for a fixed
hyperparameter $\sigma=.1$ and all $i=1,\dots,1400$. We adopt this noise model
for simplicity, as it imposes a fixed signal-to-noise ratio (SNR) of
$1/\sigma=10$ across all spectra and wavelength bins.