Understanding the difficulties of posterior predictive estimation
Authors: Abhinav Agrawal
Justin Domke
Venue: ICML 2025 poster (2025)
Keywords: Posterior predictive distribution
Variational Inference
Bayesian prediction
Primary area: probabilistic_methods->variational_inference

Abstract: Predictive posterior densities (PPDs) are essential in approximate inference for quantifying predictive uncertainty and comparing inference methods. Typically, PPDs are estimated by simple Monte Carlo (MC) averages. In this paper, we expose a critical under-recognized issue: the signal-to-noise ratio (SNR) of the simple MC estimator can sometimes be extremely low, leading to unreliable estimates. Our main contribution is a theoretical analysis demonstrating that even with exact inference, SNR can decay rapidly with an increase in (a) the mismatch between training and test data, (b) the dimensionality of the latent space, or (c) the size of test data relative to training data. Through several examples, we empirically verify these claims and show that these factors indeed lead to poor SNR and unreliable PPD estimates (sometimes, estimates are off by hundreds of nats even with a million samples). While not the primary focus, we also explore an adaptive importance sampling approach as an illustrative way to mitigate the problem, where we learn the proposal distribution by maximizing a variational proxy to the SNR. Taken together, our findings highlight an important challenge and provide essential insights for reliable estimation.

Understanding the Difficulties of Posterior Predictive Estimation
Abhinav Agrawal
Justin Domke'
Abstract
PPD, is ing predictions and forecasts (Filos et al., 2019; Kompa al., 2021; Immer et al., 2021; Petropoulos et al., 2022; Tyralis & Papacharalampous, 2022; Martin et al., 2024) and evalu- ating inference methods where higher PPD), values indicate better results (Wu et al., 2017: Wenzel et al., 2020; Izmailov et al., 2021b; Kim et al., 2022; Reichelt al., 2022b; Yu & Zhang, 2023; Zimmermann al.. 2023; Sendera et al., 2024; Cheng et al., 2024). (Note, PPD, in eq. (1) is different from the PPD under the true posterior, see section A for discussion.) The integral in eq. (1) is typically estimated via the simple Monte Carlo (MC) estimator
used extensively in approximate inference for mak-
Predictive posterior densities (PPDs) are essential in approximate inference for quantifying predic- tive uncertainty and comparing inference methods. Typically, PPDs are estimated by simple Monte Carlo (MC) averages. In this paper, we expose critical under-recognized issue: the signal-to- noise ratio (SNR) of the simple MC estimator can sometimes be extremely low, leading un- reliable estimates. Our main contribution is theoretical analysis demonstrating that even with exact inference, SNR can decay rapidly with an increase in (a) the mismatch between training and test data, (b) the dimensionality of the latent space, or (c) the size of test data relative training data. Through several examples, we empirically verify these claims and show that these factors indeed lead to poor SNR and unreliable PPD estimates (sometimes, estimates are off by hundreds of nats even with a million samples). While not the pri- mary focus, we also explore an adaptive impor- tance sampling approach as an illustrative way mitigate the problem, where we learn the proposal distribution by maximizing a variational proxy the SNR. Taken together, our findings highlight an important challenge and provide essential insights for reliable estimation.
K Σ k=1
R
p (D
where 2x ~ 9D-
K
It is typical to work in log-space and use log Rx to estimate log PPDg. From Jensen's 's incquality, we know log Rx is a biased estimator of log PPD, (with a bias related to Rx 's variance, see section A). This paper is motivated by the observation that this bias can be extremely large, leading to unreliable PPD, estimates.
Estimated PPD (nats)
820
-840
-660
Gaussion
Flow
1. Introduction
Given model with prior p(z) and likelihood p(D[=s), train- ing data D. and an approximate posterior (D(z), the predic- tive posterior density (PPD) of some test data set D* under qp is defined as
10P
10
10#
# of samples (K)
Figure 1: Unreliable PPD, estimates. log PPD, estimates for user-preference model on the MovieLens-25M dataset (Harper & Konstan, 2015). with approximate posterior qo produced from variational inference (VI) with either a Gaussian or flow-basee bas-daeed family. Lines show the mean log Rx and shaded regions show 95% intervals (uses 1000 repetitions.) Even million samples (K 10"), we do not appear to have accurate estimates.
PPD, == f ((D))2)))()))d.
(1)
'Manning College of Information and Computer Sciences, Uni- versity of Massachusetts, Amherst, MA, USA. Correspondence to: Abhinav Agrawal <aagrawal@coo @cs uda.s.
Figure shows log PPD,, estimates for user-preference model on the MovieLens-255 25M dataset (Harper & Konstan, 2015), with approximate posterior qp produced from varia- tional inference (VT) with either Gaussian or flow-based
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Understanding the difficulties of posterior predictive estimation
family (see section 5.4 for setup). Even with a million sam- ples, the curves are still increasing, indicating bias remains and comparisons between the two methods are unreliable. Such failures raise the question: Why and when is simple MC unreliable for pPD, estimation?
2. Analysis With Exact Inference
This section focuses on understanding the SNR of the naive estimator from eq. (2) under exact inference; that is, 9D(2) = p(z[D). Working with exact inference removes the confounding effects of the posterior approximations. We start by deriving SNR expressions in any model and then analyze these expressions to identify the influential factors.
The main contribution of this paper is to develop a better understanding of the PPD, estimation problem. To do this, we focus our analysis on the signal-to-noiss ratio (SNR) of Rx defined as SNR (Rx) E // / [/[Rx] This definition captures how high or low is the variance relative to the target value, and becomes crucial when the target value is small. SNR has been used in the past to study gradient estimators when gradient values can be extremely small (Roberts & Tedrake, 2008; Rainforth et al., 2018; Tucker et al., 2019; Finke & Thiery, 2019; Liévin et al., 2020; Geffner & Domke, 2021; Rudner et al., 2021). Since PPD, values are (generally) numerically small, we focus on studying SNR (see section A for more discussion),
Let D = {yı, V2..... D/D be the training data, D = [vi 3/2-. -. } be a test dataset, where y can be dis- crete, continuous, or mixed. Then, the following result gives two equivalent expressions for SNR (Rx). (Note: We use multiset notation, so D + DD is the multiset addition with 2D = 'D + DD (Costa, 2021)).
Theorem 2.1. Let Rx be as eq. (2) with exact inference:
that is, (2(2) = p(=1D). Let p(z, D) p(
Then, SNR (Rx) ) = v K /KKK/Vexl/e exp(6)² for
In section 2, we first provide a general expression for the SNR in any model assuming exact inference. Our subse- quent analyses (eq. (5), theorem 2.2, and corollary 2.3) finds that SNR decays quickly with increase in (1) the degree of "mismatch" between the test and training data, (2) the dimensionality of the latent variable space, and (3) the size of the test data relative to training data. Several illustrative examples (figs. 4 to 6) empirically validate these findings (see fig. 2 for demonstration on a linear regression model). In section 3, we also extend our analysis to approximate inference where qD # p(z[D).
8 = KL (L(p22/D (D(2z|D+ + D*) p(2
2
1
(p(2) ( +
||
2D"))
(3)
2
V(D) - V(D+ 2D") V(D + +D)
(4)
2
where V is the log- nog-orrmaiization function V(D) log ʃ pD))))p(zddz.
=
Proof sketch. A simple calculation gives SNR(R1) = vK/ / /eeplo)=- where 8 = 1 loog [[]] - log TD R1 (lemma C.1). Using p(D]2) = II p(ul=) and sim- ple algebra, gives E. R exp V(D + cD*) exp V(D) for all non-negative integers and V(D) log p(D)=)p(z)d (lemma C.3). Plugging E[R]]] with c = 1 and = 2 in 8, and simplifying gives eq. (4). Then, some simple observations give an expression for KI-divergence between two poste- riors in terms of likelihood ratios and log-normalization constants (see lemma C.4). Applying this to each KL diver- gence in eq. (3) and averaging gives eq. (4). (For formal proof, see section C.)
10a.
Exal Error (nans)
109_
Baseline
Mismatch
10 2.
Dim of
Size of D
100
102
10
106
# of samples (K)
Figure 2: Error worsens with the three factors. PPD, estima- tion error for a linear regression model. For baseline. none of the three factors influencing the SNR are high and the estimation error low. Error increases dramatically when either of three factors is high (see section 5.2 for setup).
To understand this result, note that if 8 is reasonably large, then SNR (Rx) R K cxp(-5) (see fig. 3).
SNR(R:) y = exp( -5)
SN
So, we can focus on d to un- derstand SNR. The rest of
We also explore the use of adaptive importance sampling (IS) an illustrative way to resolve poor estimation, where we learn parameterized proposal by optimizing varia- tional proxy to the SNR (section 4). This strategy, referred to as leamned IS, vastly improves estimates, including a five- fold improvement in the estimated difference between the VI methods from fig. 1 (see table 4).
this section analyzes ō as in eq. (3) to determine the in- fluential factors.
2
Figure 3: SNR rapidly de- cays with 8.
KL-divergences in eq. (3) show that SNR is determined by how different the pos- terior p(2|D + D²) is from the posteriors p(=[D) and
2

---

Understanding the difficulties of posterior predictive estimation
p(=|D 2D*). This gives us the first factor: mismatch be- tween training and test data. Intuitively, if datasets mismatch such that adding or subtracting D significantly changes the posterior p(z|D + D"), then the SNR will be small.
BI: The test features XD- consist of m copies of X and
responses yD- consist of m copies of yD + A where A is the mismatch vector of size (D).
B2: For sequence of increasingly large training datasets, posterior is dominated by data. Mathematically, as
It is important to also analyze the case when datasets match; that is, when adding or subtracting D* would not have a significant effect on p(z D D. Intuitively, if the training and test data are similar and large such that the posteriors involving D are concentrated similar to p(z[D), then the divergence in eq. (3) will be small. To understand what re- maining factors influence SNR in such conditions, we make some simplifying assumptions and provide an informal re- sult. Assume the following hold.
[D] 00, (xxXp) 5 0.
Assume B1 and B2. Then, lim 18 .timl 0 and |P|
d
1+ m
ôlim
log
2
v1 + 2m
m
+ 2
A A (XXXp) XTA. (6)
2m² + 3m +
Al: The datasets D' and D are large enough that the KL-divergence between the posteriors eq. (4) is well-approximated by the KL-divergence between the corresponding Bayesian CLT approximations.
We discuss the proof for this result in section B1 nec- essary to simplify the expressions and B2 essentially means the posteriors in eq. (3) are dominated by data. This is analogous to (but weaker than) AI. In fact, if there is no dataset mismatch (A = 0), then the right-hand side of eq. (6) reduces exactly to that of eq. (5).
A2: The datasets D. D + D", and D + 2D* are similar enough that MLE and Hessian of the average log- likelihood are the same for all three.
The expressions in eq. (6) explicitly combine the three fac- tors: mismatch (A), dimensionality (d), and relative size of D' (m). With some manipulation, the following bounds hold that bring out these relationships even more clearly (see section E. for derivation of the bounds).
Then, if d is the number of dimensions of 2, we have
d log 1 1+ 1D111//1 [D]
ô ~ 2
(5)
V1+ 2 D-1/
In section D. we provide detailed version of this result (alongside the discussion for why we need A1 than simply applying Bayesian CLT). Here, we focus on the implica- tion of eq. (5). The -att-haan--sde of the eq. (5) is well approximated by d log (/D))/1D) when |D* is large. Therefore, when datasets match, the SNR still depends on di- mensionality of the latent space and the relative sizes of the datasets: increases linearly in the number of dimensions and logarithmically in |D*|/ /D]. This completes identifi- cation of the three factors that influence the SNR: dataset mismatch, dimensionality, and relative size of D*
d m log
2 S 6pm S
$
d
log
m
+
A (7) + 40
1
2
Overall, eqs. (6) and (7) capture the strength of the three factors: 8 is affected quadratically by mismatch, linearly by the dimensionality, and logarithmically by the relative size.
(1)
10
Data mismatch (A)
100
10
2.1. Example: Gaussian Linear Regression
While theorem 2.1 provides a general expression for SNR, it required some intuition to see the impact of the three factors. With the aim to make these relationships more explicit, this section considers a linear regression model with analytically tractable Gaussian posteriors. The closed-form 6 expression for this model also remain terse and we make some simplifying assumptions (like, D* contains copies of training data) to gather more intuition. The following result provides a single expression capturing the the impact of the three factors on SNR.
150
-10
50
100
Dimensionality of (d)
Figure 4: SNR isocontours. SNR for the Gaussian linear regres- sion model where the mismatch A and the dimensionality of are varied, keeping relative size of D" fixed (m 1).
Theorem 2.2. Let p(z) = -((lpo, Ea) and p(yD|2) = (VDD X22,o21), where yD E R!D are the responses, XD E RD are the features, 2 E Rd are the weights, and 2 the known variance. Let 8 be as in eq. (3).
We consider an instance of the linear regression model from theorem 2.2 in figs. 2 and 4. Figure shows SNR decays quickly with an increase in mismatch and dimensionality. Figure 2 shows the average difference between log PPD,
3

---

Understanding the difficulties of posterior predictive estimation
and log Rg when the three factors influencing SNR are in- creased independently. As any factor increases, estimation error remains high and even using many a million sam- ples does not help much (see section 5.2 for setup details). Moreover, the steepness of the slopes in fig. correlates inversely with the relative impact of the three factors. For instance, ò scales quadratically with the mismatch, and con- sequently, the corresponding error slope is the least steep, indicating this error is harder to reduce with more samples.
that the second form for ò in terms of log partition functions (eq. (12)) allows additional insight. Figure 5 plots the SNR for three conjugate models using eq. (12).
To understand eq. (12). note that ED+D- 2ED + fD+2D-). Since is convex, 8 is the looseness in Jensen's 's inequality: the mean of ((ED) and B(ED42D-) versus B applied the mean of ED and {D+2D-
Jensen's inequality is tight when the function is nearly linear in the range evaluated. Imagine evaluating B ((aE) for > O, i.e., along ray emanating from the origin. B has "log- sum-exp" form, so as a becomes large, (a becomes nearly linear along this ray (Boyd & Vandenberghe, 2004). So, when eD and {D+2D= are large and lie along the ray, 8 is small. In summary, 5 is small (and the SNR large) when:
2.2. Analysis With Exact Inference And Conjugacy
This section considers theorem 2.1 for conjugate models. These models enjoy nice analytical properties, allowing further insights into the low SNR problem. First, we discuss some conjugate model notation and then present the result.
1. ED is large (so that B is locally "flat" near {p).
Consider an exponential family
2. ED and ED+2D* lie close to ray emanating from the origin. This happens when the statistics T(D) and T(D*) are similar, and the prior parameters Eo are either small or nearly proportional to ED.
ply z) hly) exp(T(y)) 6(2) A(2)),
A(z) log / h(v) exp((())) (Tp())](](())) d()))dy,
(8)
where h(y) is the base measure, Tyy) is the sufficient statis- tic, Φ is one-to-one parameter map, and A is the log- partition function ensuring normalization. The correspond- ing conjugate family is
Overall, corollary 2.3 corroborates the results from theo- rem 2.1 and eq. (5). The log- peg-partiiion function B behaves like a "soft-max" function or, very informally, a "rounded cone." We know it is "rounded" near the origin, but if you follow any ray away from the origin, it becomes "flat." " In these flat regions, Jensen's 's inequality is tight and SNR is high. Therefore, the impact of relative dataset size and mis- match of data statistics on SNR can be understood in terms of how they position points along the surface of the "cone" small datasets position you near the origin where the log partition function is rounded, while different moments mean the points do not lie on a ray pointing near the origin. Either of these cases lead to a poor SNR. Below, we provide an example to visually demonstrate this phenomenon.
6(2) 1 B(E)) A((=)]
x(([g) exp
T
((2) AA())
B{E) log
exp
ET
dz,
where & is the conjugate family parameter. s is "conjugate" because if the prior is p(z) = >=)(Sn) and the likelihood is p(D)=) = IL p(y|2). then posterior is within the family and given by p(z (D= = s(2(ED) where
ED = Eo + EVED T(y)
(10)
Example: A Gaussian Conjugate Model
D
Take the conjugate normal model with p(z) (2 1) and p(yl=) (v 1). Let T(D) denote the mean sufficient statistics of y E D. Take training dataset ID with 7(D) 10, a "matching" test dataset with T(D;) = 10 and a "mismatched" test dataset with T(D;) = 5. All datasets have 100 examples
Corollary 2.3. Take a model with likelihood (D)2z) 2) in an exponential family (eq. (8)) and prior p(z) = s(2(5)) in the conjugate family (eq. (9)). Let Rk be as in eq. (2). Then,
SNR (Rx) =VK/ exp(6)2 ffor
(ID) = [D]| = DD = 100).
KL (s(eEE++)) s(((Ex))
B(
B(E) along Ep and Ep
1 +2 KL (1
5(2(E
(11)
15000
0.15
2000
12200
200.6
B (ED) + B(ED+2D=) (ED+D*),
Socoud dim
(12)
1000
54
12000
2
where &D is as in eq. (10) and B is as in eq. (9).
1100
011E2)
BiCo og)
-200
This result is similar to theorem 2.1 (see section F for a
First dim of
standalone proof). The main advantage of this new result is

---

Understanding the difficulties of posterior predictive estimation
SNR(R)) for Binonial
SNR(R)) for Exp
SNR(R)) for Normal
20
0-
15
1.0
--0 (1)-20....
10
10
10-10 to
30
10
10
10
00.1
-10"
10.0
-10"
110 10
.10-
-10
50
100
30 D
100 0 50
100
Figure 5: SNR isocontours. See section for model details. For each model, we derive the log- normalizationn constant and use expressions from eq. (12) to calculate the SNR in closed form. We set the mean sufficient statistics T(D) to 10 (denoted in red dotted line) and the number of data points [D] to 1000. Setting these values theoretically is sufficient calculate the SNR since eq. (12) requires the natural parameter& & under different datasets, which can be evaluated using eq. (10). Dataset mismatch increases as we move away from the red dotted line, and the relative size of "D" increases as we move along the horizontal axis. Either way, SNR decreases exponentially.
when using (z(2) as a prior and conditioning on D When inference is exact, (D(2|D*) p(z + DD) and eq. (13) reduces to eq. (3).
Figure 6: Left: Isocontours of the log partition function B(E) (eq. (9). Right. The values of B(E) along the lines joining ED to SD42; and ED+
Similarly, eq. (14) is generalization of eq. (4). To see this, write 6 = 1/2 (Z) + Zp(2D") - (p(D") where ( log (D())d2 = O. When inference is exact, sim- ple manipulations make the two expressions equal.
Figure 6 shows B(E) along with the values of natural pa- rameter & for each dataset. Notice how {D. SD+D; and ED+20; are equidistant on "ray" pointing to near the origin (left panel), meaning Jensen's 's inequality is nearly tight (right panel), but the line joining &D and CD+D; does not point towards the origin, meaning Jensen's 's in- equality is not tight, resulting in Rs 200.56 and an astronomically small SNR of SNR (R1) Re 7.9 x 10 88,
Next corollary specializes theorem 3.1 to the case of conju- gate models. For simplicity, we assume that approximate distribution lies in the conjugate family.
Corollary 3.2. Let p(D|2) and p(2) be as in corollary 2.3. Let (p(z) = s(2(n) be the conjugate family (eq. (9)) with parameters n, and let Rg be as in eq. (2). Then, (vR(Rx) = vK/ /exp(8)2 for
3. Analysis With Approximate Inference
In the previous section, we assumed that inference was exact, i.e., (2) P(=[D). This allowed us to uncover the factors that affect the SNR. However, in practice, methods like VI or Laplace's s approximation are not exact. This section generalizes our analysis to cases where 4p may not be the same as p)z We start by generalizing theorem 2.1 and then specialize to conjugate models.
8 1 (L( + U(')) s
+ (L( + | 0(11++U000 U((2")
(15)
(n) + B(n + U(2D") B (1 UD')) (16)
2
where B is as in eq. (9) and UD) [T(D),
Theorem 3.1. Let Rx be in eq. (2). Then, SNR (Rx)= vK / exp forr
See section H for a proof. This result has the same func- tional form as corollary 2.3 and differs only in the canonical parameters. Now, 17 are the parameters of qp and n + (D*) are the parameters of the posterior obtained by condition- ing on D* with qp as prior. When the inference is exact, 1) = ED. the above expressions reduce to corollary 2.3.
1
8- KL (gD(>|D") l1 [D(2))
2
+ KL (aD( (z ( qD(2|2DD))
(13)
2
2 (2D") - (D"),,
(14)
Note ö as in eq. (16) is again the looseness in Jensen's inequality: mean of B(n) and B(n + U(DD") versus B applied to mean of and + U(DD)). So, when test data statistics UD*) are such that 1) and 17 + U(2D") lie along a ray originating at the origin, the Jensen's 's inequality in eq. (16) is tight and 8 is small. Otherwise, we expect 8 to be large and the SNR to be low.
where (D*)) = log f | (2)g0 and (=(2 (D*) x (()2))))
For a standalone proof, check section G. As before, eq. (13) determines 8 in terms of divergences, but these distributions are unusual. One may think of (p<z|D*) as the posterior
5

---

Understanding the difficulties of posterior predictive estimation
Overall, theorem 3.1 and corollary 3.2 tell us that when approximate posterior closely resembles the true posterior, the relationships from section 2 hold as is. However, for arbitrary approximations, we can only reason about the SNR in terms of how much the "posteriors qD (z(2") and QD (=(222D") differ from the "prior (z(2)" as making more precise statements in terms of the datasets, like we did in section will require specific assumptions on qD.
or equivalently, maximizing SNR (R!) (Maddison et al., 2017; Dieng et al., 2017; Rainforth et al., 2018; Domke & Sheldon, 2018). More formally,
lim M (log PPD, - IW-ELBO,r) M-
-V[P9)/(PPD))) (19)
So, optimizing the IW-ELBO in eq. (18) can be thought of as surrogate for optimizing the SNR of the IS estimator. The naive gradient estimator of IW-ELBO also has poor SNR (Rainforth et al., 2018; 8: Finke & Thiery, 2019). Fortunately, recently proposed re-parameterized gradient estimator circumvents this issue (Tucker et al., 2019; Finke & Thiery, 2019) and offers practical option (Agrawal et al., 2020).
4. Learned Importance Sampling
The previous sections analyzed the SNR expressions for the simple MC estimator and found it decays quickly with increase in mismatch, dimensionality, and relative size of D' When the SNR is poor, even with a large number of samples the estimation error remains high (fig. 2). This section considers an instance of the adaptive importance sampling idea (Owen, 2013, Chapter 10) as an illustrative way to mitigate the poor estimation error.
We refer to our specific instance of adaptive IS as the Learned IS (IIS).. Figure provides the pseudocode.
TearnedTs(D", K)
In general, when an MC estimator has high variance, a standard solution is to replace it with importance sampling (IS). For a valid proposal r, the IS estimator for PPD, can be written
w <- Op-ipimmmii ze (IW-ELBO)
26 Tw vk E {1, (1..... K}
R K
Σ p(D [x (i)
PDD ((D*)=x))D(zz) r(zx)
Figure 7: Evaluating PPD, with Learned IS.
R
where 2 ~r r. (17)
1=
The choice of the proposal is delicate. Setting r(2) = (2) reduces eq. (17) to the simple MC estimator, while using Ops x p (2) the optimal IS -esulls in infinite SNR (Rainforth et al., 2020). The idea of adaptive IS is to somehow learn a good proposal (since Opt is often intractable). Converting this idea into an algorithm is non- trivial as it requires making several decisions, like choosing the family of proposal distributions, an update scheme, and a stopping criterion (Owen, 2013, Algorithm 10.1).
5. Examples
We consider four models: exponential family, linear re- gression, logistic regression, and a hierarchical model (see section B for more details of the setups). For each model, we use the simple MC estimator and demonstrate whenever either of the three factors is high, SNR is poor and PPD, estimates are biased even after using large number of sam- ples. We then use the learned IS estimator from fig. 7 and show that it enjoys sigmificantly higher SNR and provides accurate estimates with fewer samples.
One straightforward way to learn r, based on our analysis from sections 2 and 3, is to learn a parameterized proposal Tw. with parameters w, by optimizing the SNR of the es- timator in eq. (17). However, maximizing SNR (R() is equivalent to minimizing the variance of RIS. which in turn is equivalent to minimizing the x² divveregee between Opt and rw (Dieng et al., 2017), and recent research suggests gradients for x2 themselves suffer from poor SNR, making this direct optimization difficult (Geffner & Domke, 2021).
5.1. Exponential Family Models
Tables 1 and 2 show the results of PPD, estimation under exact inference and approximate inference, respectively. In both the cases, simple MC estimator suffers from low SNR and high estimation error even after using a million samples. The empirical SNR of R is much higher than Rx. Under exact inference, R is deterministically equal to PPD,, and under approximate inference, the learned IS estimates are hundreds of nats higher than the simple MC estimates.
We take an alternative approach and learn Fur by optimizing the importance weighted evidence lower-bound ([W-ELBO)
(Burda et al., 2016). Let 2m Tw. Then
For the proposal and the variational families, we use full- rank Gaussians. For learned IS, we set M = 16 for the IW-ELBO and optimize for 1000 iterations with Adam (Kingma & Ba, 2015) and a learning rate of 0.001. (See section for details, and table 6 for computation of log R and SNR (R) in tables and 2).
IW-ELBOM m (rw p(D" ((DD())) == log
1
Σ A
Em=1
p(D" 2l (1D))mm)9 [9p(mm)
(18)
M
Puo (3mm)
is known that maximizing IW-ELBO in eq. (18) is asymp- totically equivalent to minimizing the variance of RiS.
6

---

Understanding the difficulties of posterior predictive estimation
Table I: Results for log , estimation under exact inference (see table 5 for dataset details). We use K = 10 for simple
MC and K = 10 for learned IS. Mean and standard deviation reported over ten runs. (See section B for model details.)
Model
Ground truth log PPD,
Lower bound on log PPD, (1)
Signal-to-noise ratio (1)
Simple MC.
Learndd IS
Simple MC
Learned IS (SNR(R)))
(Ellog Rx))
(Ellog R
(SNR(Rx))
Normal
-774.64
-1183.98 0.34
-774.64 + 0.00
0.35 + 0.02
76.5 + 23.43
Exp
-527.44
-559.98 0.16
-527.44 + 0.00
0.34 + 0.01
222.76 + 147.59
Binomial
-327.42
-487.13 1.29
-327.41 + 0.00
0.03 + 0.00
173.06 + 104.2
Table 2: Results for log pPD, estimation under approximate inference (see table for dataset details). We use K 10 for simple MC and K = 10 for learned IS. Mean and standard deviation reported over ten runs.
Model
Ground fruth log PPD,
Lower bound on log PPD, (t)
Signal-to-noise ratio (1)
Simple MC
Leamed IS
Simple MC (SNR(Rx))
Learned IS (SNR(R)))
(Ellog Rx))
(Ellog R
Normal Exp Binomial
1194.3 0.41
-775.23 + 0.00
0.35 + 0.02
238.79 172.46
-57..27 + 0.14
-542.34 + 0.00
0.00
0.36 + 0.01
215.09 140.52
332.46 + 0.74
-322.66 0.00
0.34 + 0.02
70.13 + 35.29
Table 3: Results for log PPD, estimation for logistic regression (section 5.3). We use K 10 for simple MC and K 103 for IS estimators. Mean and standard deviation reported over ten runs.
Model
Ground truth log PPD
Lower bound on log PPD, (t)
Signal-to-noise ratio (t)
Simple MC
Learned
Simple MC
Learned IS
(Bllog Rx))
(Ellog R
(SNR(Rx)
(SNR (R
Baseline
-525.25 + 0.01
-525.12 +: 0.00
1.35 + 0.20
645.67 + 14.99
More dimension
-702.07 = 0.28
-543.00 + 0.00
0.04 0.01
57.78 1.37
More mismatch
-1687.98 = 0.96
-734.32 + 0.00
0.03 0.00
728.63 + 16.01
More test data
-5143.69 + 0.34
-5097.60 + 0.00
0.04 + 0.01
802.34 + 18.37
Table Results for log PPD, estimation for Moviel 25M dataset (section 5.4). Mean and standard deviation reported
over ten runs.
Method
No. of samples (K)
Lower bound on log PPD, (t)
Signal-to-noise ratio (t)
Simple MC (E[log Rx])
Learned IS
Simple MC (SNR(Rx))
Learned IS
(E[log R())
(SNR(RS)
K = 103
-796.24 + 0.13
-779.39
0.02
0.05 + 0.02
0.11 0.04
Flow VI
K = 106
-787.27 + 0.08
-777.73 = 0.01
0.04 + 0.01
0.48 0.29
K 103
-828.22 + 0.17
-783.89 + 0.03
0.04 + 0.00
0.12 + 0.01
Gaussian VI
K = 106
-811.61 + 0.13
-781.88 = 0.02
0.04
0.32 + 0.13

---

Understanding the difficulties of posterior predictive estimation
Baseline
More Mismatch
More test points
High dimensions
Eval Error [nats)
109
10
Naive MC LIS
10-3.
10-6
10
10a
10
10
10
10a
10a
10
10P
10a
10
10
10
10a
10
10
K
K
K
K
Figure 8: Estimation Error. We plot the estimation error, log PPD, - log Rx, against the number of samples K for the linear regression model (section 5.2). The fifth and the ninety-fifth percentiles are represented as the filled regions. Across the scenarios, the estimation error worsens as either of the three factors increase. LIS significantly reduces the error. See section K for more details on scenarios.
5.2. Linear Regression
SNR of these estimates is extremely low. With learned IS, flow VI is only 4 nats higher than the Gaussian VI, with much higher SNR (fourth column). So, while flow VI may be better, the difference is not as large as it seems initially, (See section M for setup details.)
We take the linear regression model with exact inference as described in theorem 2.2 and construct four scenarios. We start with a baseline where none of the three factors influ- encing SNR are too high. We then independently increase mismatch, dimensionality, and relative size (see section K).
6. Discussion
In fig. 8, plot the error in evaluating log PPD, under exact inference using simple MC and learned IS. For the baseline (first panel), the simple MC has high enough SNR, and the estimation is accurate for K- = 10 This error in- creases dramatically when either of the factors increase (see red curves in last three panels). The learned IS consistently evaluates accurately across all scenarios (see blue curves). (See section K for details of learned IS settings).
Related Work. Wu et al. (2017) explored the use of An- nealed Importance Sampling (AIS) (Neal, 2001) for estimat- ing the posterior predictive density in decoder based models. In particular, they used AIS for estimating the normalization constant of the unnormalized density p(y, |=))g(z)) (D) for each data point yı in the test data set D Different from them, we focus on black-box treatment of probabilistic mod- els (Ranganath et al., 2014; Kucukelbir et al., 2017) and exploit BBVI schemes (Kucukelbir et al., 2017; Agrawal et al., 2020) for estimating the posterior predictive densi- ties over D* Recent theoretical advances (Domke, 2019; 2020; Domke et al., 2024; Kim et al., 2023; 2024b:a) make BBVT general purpose inference method that is reliably applicable to a wide range of problems (Carpenter et al., 2017).
5.3. Logistic Regression
We construct four scenarios similar to section 5.2 -th base- line where none of the three factors are too high and addi- tional scenarios where each factor is increased individually.
Table 3 reports the results from estimating PPD, using the simple MC estimator and the learned IS estimator. For baseline, the simple MC estimator has reasonably high SNR but for other scenarios it suffers from low SNR and poor estimation error. The learned IS consistently provides higher SNRs and better estimates. See section L. for setup details.
Other research has explored learning approximate posterior distributions ID to calibrate for test-time utilities (Stoy- anov et al., 2011; Lopez et al., 2020; Lacoste-Julien et al., 2011: Morais & Pillow, 2022; Kušmierezyk et al., 2019; Knoblauch et al., 2022; Kušmierczyk et al., 2020; Vadera et al., 2021). Such methods aim to learn distribution q' that is different from ypD and optimizes the expectation of some utility function under gD at test-time. We focus on the problem of estimating the posterior predictive density for a given go at test-time, and do not change the given posterior; we simply focus on accurate estimation.
5.4. Hierarchical Model
We consider hierarchical user- uee--reffeeenc model for the MovieLens 25M dataset. Here, we have dataset shift and the relative size is small ((D*)/[D) 0.1). However, the di- mensionality is moderately high (d 1065) and simple MC estimator can still suffer from low SNR. We consider two posterior variational approximations-fuus full-rank Gaussians and normalizing flows-for learning (p.
Ruiz et al. (2016) explored learning an importance sampling estimator for estimating the gradients for BBVI (Ranganath et al., 2014). They learn a proposal distribution r while learning the parameters of the variational distribution qp. and rely on exponential families for closed-form updates. We do not focus on learning the variational distribution JD. and use BBVI methods for learning the proposal that
Table reports log PPD, estimates using simple MC and learned IS for different values of K. When using the simple MC with K- 10%, flow VI reports PPD, more than 20 nats higher than Gaussian VI (second column). However, the
8

---

Understanding the difficulties of posterior predictive estimation
can be in any suitable family of distributions (Rezende & Mohamed, 2015; Papamakarios et al., 2021; Webb al., 2019; Agrawal et al., 2020).
Kristiadi et al. (2022) hypothesize that the inaccuracy of the posterior predictive estimation is related to the inaceuracy of the approximate posterior. The authors explore this empiri- cally and consider the use of normalizing flows to get better posterior approximations, We find this work orthogonal to ours. We lay out the conditions that lead to poor posterior predictive estimations (including when inference is exact) and provide a post hoc method for improving the estimation without changing the approximate posterior.
Vehtari et al. (2016) evaluate predictive accuracy using met- rics that involve leave-one-out "point wise" predictive den- sity of the type pDu))D_)) over the training data D. To estimate p (D-)) = J' ((Nl=))((>
the authors consider using the full posterior distribution p(z)D) as the proposal. However, p(=(D) can have thinner tails than p(=|D=,) leading to large importance weights. To remedy this, the authors fit Pareto distribution to the importance weights and then use statistics from the fitted distribution for final estimation. While the PSIS-LOO setting differs from our focus, one can use PSIS ideas to improve LIS estimates if ( is suspected of thin tails.
Conclusions. This paper develops theoretical understand- ing of why the SNR of the naive PPD estimator can be extremely poor. We demonstrate that even if the inference is exact, SNR suffers when either of the following quanti- ties increases: mismatch between the training and test data, dimensionality of the latent space, or the size of the test data relative training data. We hope work serves as strong cautionary note and more practitioners monitor (and report) the SNR value of their PPD estimates. Moreover, whenever the SNR values of the naive estimator is low, our work shows that this does not necessarily reflect on the accu- racy of the approximate posterior but simply on our ability to accurately estimate PPD. In cases of low SNR, we sug- gest practitioners use approaches like LIS improve the reliability of the estimates.
Rainforth al. (2020) propose framework for target-aware Bayesian inference (TABI) in which they decompose the posterior expectations into three components. Each of the three components is then computed as an instance of impor- tance sampling using the Annealed Importance Sampling (AIS) or Nested Importance Sampling (NIS). One can ap- ply the TABI framework for PPD, estimation; however, after some simple observations, this reduces estimating I' [(D ())))=)))))))))))))))))))>)))))) with AIS or NIS (and is same as the ap- proach from Wu et al. (2017)). In recent work, Llorente et al. (2023) extend the TABI framework by employing the generalized thermodynamic integration scheme (GIS) for solving the posterior expectations. When placing these TABI approaches in context, it is crucial to note that we focus approximate inference problems. Running MCMC procedures like AIS or thermodynamin integrations proce- dures like GIS is often infeasible or extremely slow on such problems (due to a large number of data points or dimen- sions), and view the MCMC procedures as orthogonal to our variational approach.
Limitations. The primary contribution of this work is the- oretical. Our explorative strategy performs favorably and can be developed further in the future work. The proposal learning can be computationally expensive when the test data large. Further, learned IS requires evaluating the approximate distribution and does not immediately extend to Markov chain Monte Carlo methods.
A Aeknowledgement
The authors would like to thank Edmond Cunnigham, Mohit Yadav, and Miguel Fuentes for useful discussions and com- ments, and Michelle Hedlund and Ankit Patel for helping improve the clarity of the manuscript.
Reichelt al. (2022a) propose the concept of expectation programming, where probabilistic programming system considers the target posterior expectation as a first-class citi- zen. They aim to build an efficient estimation pipeline when target functions are previously known. In their implemen- tation, they currently use Annealed Importance Sampling as the choice of inference scheme. Our proposed methodol- ogy can join their suite of inference options when the target functions are more amenable to variational formulation.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.
Izmailov et al. (2021a) point out that the posteriors in Bayesian neural network can be bad at generalizing un- der specific dataset shifts. They uncover pathologies in the BNN posteriors that lead to poor generalization and present techniques that can possibly mitigate these. Different from them, we focus on understanding the problem of inaccurate PPD, estimation and how to improve estimation without changing the properties of the posterior.
References
Agrawal, A., Sheldon, D. R.,, and Domke, J. Advances in black-box VI: normalizing flows, importance weighting, and optimization. In NeurlPS, 2020.
Ambrogioni, L.,, Lin, K., Fertig, E., Vikram, S., Hinne,

---

Understanding the difficulties of posterior predictive estimation
M., Moore, D., and Gerven, M. Automatic structured variational inference. In AISTATS, 2021a.
Finke, A. and Thiery, A. H. On importance- toencoders. arXiv preprint arXiv: (90dv::/977.0444440 10477, 2019.
Ambrogioni, L., Silvestri, G., and van Gerven, M. Au- tomatic variational inference with cascading flows. In ICML, 2021b.
Geffner, T. and Domke, J. On the difficulty of unbiased alpha divergence minimization. In ICML, 2021.
Harper, F. M. and Konstan, J. A. The movielens datasets: History and context. Acm transactions on interactive intelligent systems (tiis), 5(4): 2015.
Bingham, E., Chen, J. P., Jankowiak, M., Obermeyer, F,, Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P., Hors-
fall, P., and Goodman, N. D. Pyro: Deep universal proba-
bilistic programming. JMLR, 2019.
Immer, A., Korzepa, M., and Bauer, M. Improving predic- tions of bayesian neural nets via local lincarization. In AISTATS, 2021.
Boyd, S. and Vandenberghe, L. Convex optimization. Cambridge university press. 2004.
Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman- S., and Zhang, Q. JAX: composable transformations of Nython+NumPy programs, 2018. URL http::/gitee com/ haopa//oool////ooooooooooooeeeee jax.
Izmailov, P., Nicholson, P., Lotfi, S., and Wilson, A. G. Dangers of bayesian model averaging under covariate shift. In NeurlPS, 2021a.
Izmailov, P., Vikram, S., Hoffman, M. D., and Wilson, A. G. G. What are bayesian neural network posteriors really like? In ICML, 2021b.
Burda, Y., Grosse, R., and Salakhutdinov, R. Importance weighted autoencoders. ICLR, 2016.
Kim, K., Oh, J., Gardner, J., Dieng, A. B., and Kim, H. Markov chain score ascent: A unifying framework of vari- ational inference with markovian gradients. In NeurIPS, 2022.
Burroni, J., Domke, J., and Sheldon, D. Sample average approximation for black-box vi. In UAI, 2024.
Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, and Riddell, A. Stan: probabílistic programming language. Journal of statistical software, 2017.
Kim, K., Wu, K., Oh, J., and Gardner, J. R. Practical and matching gradient variance bounds for black-bon variational bayesian inference. In ICML, 2023.
Cheng, Z., Yu, L., Xie, T., Zhang. S., and Zhang, C. Kernel semi-implice variational inference. ICML, 2024.
Kim, K., Ma, Y., and Gardner, J. Linear convergence of black-box variational inference: Should we stick the land- ing? In AISTATS, 2024a.
Costa, L. F. An introduction to multisets, arXiv preprint arXiv:21 10. 12902, 2021.
Dieng, A. B., Tran, D., Ranganath, R., Paisley, J., and Blei, D. Variational inference via upper bound minimization. In NewrlPS, 2017.
Kim, K., Oh, J., Wu, K., Ma, Y., and Gardner, J. On the convergence of black-box variational inference. In NeurlPS, 2024b.
Domke, J. Provable gradient variance guarantees for black- box variational inference. In NeurlPS, 2019.
Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. ICLR, 2015.
Domke, J. Provable smoothness guarantees for black-box variational inference. In ICML, pp. 2587-2596. PMLR, 2020.
Knoblauch, J., Jewson, J., and Damoulas, T. An optimization-centri view on bayes' rule: Reviewing and generalizing variational inference. JMLR, 2022.
Domke, and Sheldon, D. Importance weighting and variational inference. In NeurlPS, 2018.
Kompa, B., Snoek, J., and Beam, A. L. Second opinion needed: communicating uncertainty in medical machine learning. NPJ Digital Medicine, 2021.
Domke, J., Gower, R., and Garrigos, G. Provable conver- gence guarantees for black-box variational inference. In NeurlPS, 2024.
Kristiadi, A., Eschenhagen, R., and Hennig, P. Posterior refinement improves sample efficiency in bayesian neural networks. NeurlPS, 2022.
Filos, A., Farquhar, S., Gomez, A. N., Rudner, T. G., Ken- ton, Z. Smith, L., Alizadeh, M., De Kroon, A., and Gal, Y. systematic comparison of bayesian deep learning ro- bustness in diabetic retinopathy tasks. In /CML Workshop on Bayesian Deep Learning, 2019.
Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A., and Blei, D. M. Automatic differentiation variational infer- ence. JMLR, 2017.
10

---

Understanding the difficulties of posterior predictive estimation
Kušmierezyk, T., Sakaya, J., and Klami, A. Variational bayesian decision-making for continuous utilities. In NeurlPS, 2019.
C., Cyrino Oliveira, L.., De Baets, S., Dokumen-
tov, A., Ellison, J., Fiszeder, P., Franses, P. H., Frazier, D. T., Gilliland, M., Gönüi,, M. S., Goodwin, P., Grossi, L., Frushka-Cockayne, Y., Guidolin, M., Guidolin, M., Gunter, U., Guo, X., Guseo, R., Harvey, N., Hendry, D. F, Hollyman, R., Januschowski, T., Jeon, J., Jose, V. R. R., Kang, Y., Koehler, A. B., Kolassa, S., Kourentzes, N., Leva, S., Li, F., Litsiou, K., Makridakis, S., Martin, G. M.. Martinez, A. B., Meeran, S., Modis, T., Nikolopoulos, K., Önkal, D., Paccagnini, A., Panagiotelis, A., Panapakidis, I., Pavía, J. M., Pedio, M., Pedregal, D. J., Pinson, P., Ramos, P., Rapach, D. E., Reade, J. J., Rostami- B., Rubaszek, M., Sermpinis, G., Shang. H. L., Spiliotis, E., Syntetos, A. A., Talagala, D., Talagala, S., Tashman, L., Thomakos, D., Thorarinsdottir, T., Todini, E.. Trapero Arenas, J. R., Wang, X., Winkler, R. L.., Yusupova, A.. and Ziel, F. Forecasting: theory and practice. Interna- tional Journal of Forecasting, 2022.
Kušmierezyk, T., Sakaya, J., and Klami, A. Correcting predictions for approximate bayesian inference. AAAI, 2020.
Lacoste-Julien, S., Huszár, F., and Ghahramani, Z. Ap- proximate inference for the loss-calibrated bayesian. In AISTATS, 2011.
Liévin, V., Dittadi, A., Christensen, A., and Winther, O. Optimal variance control of the score-function gradient estimator for importance-weigheer bounds. NeurlPS, 2020.
Llorente, F., Martino, L., and Delgado, D. Target-aware bayesian inference via generalized thermodynamic inte- gration. Computational Statistics, 2023.
Phan, D., Pradhan, N., and Jankowiak, M. Composable effects for flexible and accelerated probabilistic program- ming in numpyro. In NeurIPS Program Transformations for Machine Learning Workshop, 2019.
Lopez, R., Boyeau, P., Yosef, N., Jordan, M., and Regier, J. Decision-making with auto-encoding variational bayes. In NewrlPS, 2020.
Lotfi, S., Izmailov, P., Benton, G., Goldblum, M., and Wil- son, A. G. Bayesian model selection, the marginal like- lihood, and generalization. In ICML, pp. 14223-14247. PMLR, 2022.
Rainforth, T., Kosiorek, A., Le, T. A., Maddison, C., Igl, M., Wood, F., and Tch, Y. W. Tighter Variational Bounds Not Necessarily Better. In ICML, 2018
Rainforth, T., Golinski, A., Wood, F,, and Zaidi, S. Target aware bayesian inference: How to beat optimal conven- tional estimators. JMLR, 2020.
Maddison, C. J., Lawson, J., Tucker, G., Heess, N., Norouzi, M., Mnih, A., Doucet, A., and Teh, Filtering varia- tional objectives. In NeurIPS, 2017.
Ranganath, R., Gerrish, S., and Blei, D. Black Box Varia- tional Inference. In AISTATS, 2014.
Martin, G. M., Frazier, D. T., Maneesoonthomn, W., Loaiza- Maya, R., Huber, F., Koop, G., Maheu, J., Nibbering, D., and Panagiotelis, A. Bayesian forecasting in economics and finance: A modern review. International Journal of Forecasting, 2024.
Reichelt, T., Golinski, A., Ong. L.,, and Rainforth, T. Ex- pectation programming: Adapting probabilistic program- ming systems estimate expectations efficiently. In UAI, 2022a.
Morais, M. and Pillow, J. W. Loss-calibrated expectation propagation for approximate bayesian pacissionmmkinn. arXiv preprint arXiv:2201.03118, 2022.
Reichelt, T., Ong, L., and Rainforth, T. Rethinking varia- tional inference for probabilistic programs with stochastic support. In NeurlPS, 2022b.
Neal, R. M. Annealed importance sampling. Statistics and computing, 2001.
Rezende, D. and Mohamed, S. Variational inference with normalizing flows. In ICML, pp. 1530-1538. PMLR, 2015.
Owen, A. B. Monte Carlo theory, methods and examples. 2013.
Roberts, J. and Tedrake, R. Signal-to-noise ratio analysis of policy gradient algorithms. NeurIPS, 2008.
Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., and Lakshminarayanan. B. Normalizing flows for probabilistic modeling and inference. JMLR, 2021.
Roeder, G., Wu, Y., and Duvenaud, D. Sticking the landing: Simple, lower-variance gradient estimators for variational inference. In NIPS, 2017.
Petropoulos, F., Apiletti, D., Assimakopoulos, V., Babai, M. Z., Barrow, D. K., Ben Taieb, S., Bergmeir, C.,, Bessa, R. J., Bijak, J., Boylan, J. E., Browell, J., Carnevale, C., Castle, L.,, Cirillo, P., Clements, M. P., Cordeiro,
Rudner, T. G., Key, O., Gal, Y., and Rainforth, On signal-to-noiss ratio issues in variational inference for deep gaussian processes. ICML., 2021.

---

Understanding the difficulties of posterior predictive estimation
Ruiz, F. J., Titsias, M. K., and Blei, D. M. Overdispersed black-box variational inference. In UAI,20 2016.
Wu, Y., Burda, Y., Salakhutdinov, R., and Grosse, On the quantitative analysis of decoder- generative models. ICLR, 2017.
Sendera, M., Kim, M., Mittal, S., Lemos, P., Scimeca, I.., Rector-Brooks, J., Adam, A., Bengio, Y., and Malkin, N. Improved off off-policy training of diffusion samplers. In NeurlPS, 2024.
Yu, and Zhang, C. Semi-implicit variational inference via score matching. ICLR, 2023.
Zimmermann, H., Lindsten, F., van de Meent, J.-W., and Naesseth, C. A. A variational perspective on generative flow networks. TMLR, 2023.
Stoyanov, V., Ropson, A., and Eisner, J. Empirical risk minimization of graphical model parameters given ap- proximate inference, decoding, and model structure. AISTATS, 2011.
Tucker, G., Lawson, D., Gu, S., and Maddison, C. J. Dou- bly reparameterized gradient estimators for monte carlo objectives. ICLR, 2019.
Tyralis, H. and Papacharalampous G. review of predictive uncertainty estimation with machine learning. Artificial Intelligence Review, 2022.
Vaart, A. W. v. d. Bayes Procedures, pp. 138-152. Cam- bridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 1998.
Vadera, M. P., Ghosh, S., Ng, K., and Marlin, B. M. Post- hoc loss-calibration for bayesian neural networks. In UAI, 2021.
Vehtari, A., Gelman, A., and Gabry, J. Practical bayesian model evaluation using leave-one-out cross- vooss-validation and waic. Statistics and Computing, 2016.
Vig. J., Sen, S., and Riedl, The tag genome: Encoding community knowledge support novel interaction. ACM Transactions Interactive Intelligent Systems (TiiS), 2 (3):1-44, 2012.
Wang, C., Sun, S., and Grosse, R. Beyond marginal uncer- tainty: How accurately can bayesian regression models estimate posterior predictive correlations? In Interna- tional Conference Artificial Intelligence and Statistics, pp. 2476-2484. PMLR, 2021.
Webb, S., Chen, P., Jankowiak, M., and Goodman, N. Improving automated variational inference with normal- izing flows. In JCML Workshop on Automated Machine Learning, 2019.
Wenzel, F., Roth, K., Veeling, B. S., Swikatkowski, J., Tran, L., Mandt, S., Snoek, J., Salimans, T., Jenatton, R., and Nowozin, S. How good is the bayes posterior deep neural networks really? ICML, 2020.
Williams, D. R., Rast, P., Pericchi, L. R., and Mulder, J. Comparing gaussian graphical models with the poste- rior predictive distribution and bayesian model selection. Psychological methods, 25(5):653, 2020.
12

---

Understanding the difficulties of posterior predictive estimation
A. Additional Comments
Here we offer some additional comments on the different aspects of our formulation that we could not cover in the main text.
A.1. True posterior predictive distribution vs. PPD,
In this paper, we are not interested in estimating the posterior predictive distribution under the true posterior ( (>D*)))p(=(()dd We are interested in estimating the posterior predictive distribution under the approximate posterior as defined in eq. (1). Accurate measurements of the PPD, are of independent interest and are used for forecasting, predictions, and inference method comparisons (see citations in discussion after eq. (1)).
A.2. Relationship of estimation error and variance of Rg
Let G be an unbiased estimator of the PPD, such that [[ pPD,- Then, using Jensen's inequality, we know that Ellog G] < log gpD making log G biased estimator of log PPD,- This bias the looseness in Jensen's inequality is directly related to the variance inversely related to the SNR of G. To see this, consider the K sample estimator Gx. A second-order Taylor series approximation in the vicinity of log PPD, (as done in the second-order delta method) gives the following about the bias or the looseness in Jensen's 's inequality:
1 VIG
1
log PPD, Elloo Gk] &
2K PPD
2K SNR(G)
These relationships hold when G is unbiased and apply to the estimators in eqs. (2) and (17). Therefore, we want the SNR to be high for lower estimation errors.
A.3. SNR vs variance of Rx
In this paper, we use signal-to to-ial--t-noieee ratio (SNR) of Rg defined as SNR (Rx) [[Rk]/VV[RR We study SNR because is equivalent to relative variance and bakes in the idea of how large the estimator variance is relative to the target quantity. The idea of relative variance becomes crucial when the target quantity is numerically small, as in the case of PPD, values. To make this precise, let's consider an example where the log PPD, = 100 (for reference, all estimates of log PPD, in tables In the absolute sense, the variance of this estimator is low; however, R carries more noise than signal. To why, note that PPD, exp(-100) 3.72 x 10-44 Intuitively, R varies on the scale of 10-10 (standard deviation) and will produce noisy approximations of the target value that is the order of 10-11, SNR naturally captures this intuition: SNR(R) = / /VD,/VV(RR) s 3.72 x 10-34 and flags the estimator as poor.
4 are lower than -100). Also, consider an unbiased estimator R with variance V(R) 10 20
A.4. Joint vs. Marginal PPD
The usage of "joint" PPD (as in eq. (1)) depends on both the model and the modeler's preference. The joint PPD (J (1D")=))p(z)d2) is natural metric whenever the test data are not IID (Williams et al., 2020) or one wants to measure correlations between new data points (Wang et al., 2021). It also shows up in the conditional log marginal likelihood metric, where a fraction observations used to evaluate the conditional likelihood over the remaining observations (Lotfi et al., 2022). The "margina" PPD (ʃ p(y")2)))) (2))))2)d)))) (z)d2) is a natural metric when one wants to measure the model's ability to predict single data point. Moreover, the dependence of SNR on the size of the test data is relatively weaker than the dependence on the mismatch and dimensionality (eqs. (6) and (7)).
B. Model Details
This section provides high-level overview of the models considered throughout the paper. For the first three models, we use synthetic data sampled from the model. For the hierarchical model, we use the MovieLens 25M dataset. For more details of the setup, see sections to M.
Exponential Family Models. We consider three exponential family models. A binomial model where p(ylz) is a binomial distribution with a known number of trials n and unknown success probability & E [0, 1], and p(z) is a beta distribution. An exponential model where p(y|=) ) is an exponential distribution with the unknown rate z E RR,, and p(z) a gamma
13

---

Understanding the difficulties of posterior predictive estimation
distribution. A normal model where p(y]=) is normal distribution with known variance o² and unknown mean E R, and p(z) is also normal distribution. See table 9 in section for full model details.
Table Summary of the data sets used for results tables and 2.
Model
T(D)
D
T(D*)
|D*|
Normal
10.08
100
4.96
100
210.85
Exp
7.00
100
39.37
100
11.74
Binomial
8.96
100
41.06
100
23.32
Linear Regression. We set p(y|=) = (wl== 1) and p(z) (= I), where E Rd and y E IR model from section 2.1. We use this model for the example in figs. and and experiments in section 5.2. For details on the scenarios in figs. and 8, see section K
Logistic Regression. The non-conjugate model with p(y|2) = (ly)sigmoid(=
(=Tx)) z)) and p(=) = \=200, J) where response E {0, 1,, latent variable IR4, and feature vector x E R4. We use this model in section 5.3 experiments. For approximate inference, we use full-rank Gaussian variational posterior. For more details, see section
Hierarchical Model. The hierarchical model using the MovieLens 25M-a dataset of 25 million movie ratings along with a set features for each movie (Vig et al., 2012). We randomly select 100 users after filtering those with more than 1000 ratings. We keep one-tenth the ratings for test data, and PCA the features to ten dimensions (see section for more details).
The task is model rating Vio 10, of user i for movie with given features Ziy- We use the model
p(0, v/) (O I) ((N))())))))) B(y, (osn
|sggmodd() xx,))
(20)
where and À together represent latent variables 2; are the global latent variables capturing preferences over users and A the local latent variables capturing preferences for user lr and Σ. are functions such that if 0 [0,.. 0z., m=() 0,
and (0) tri) tril(0z), where tril takes an unconstrained vector and outputs lower-triangula, positive definite matrix. n; is the number of ratings for user i. B is the Bernuulii distribution. We use this model for illustration in fig. and experiments in Section 5.4.
14

---

Understanding the difficulties of posterior predictive estimation
C. Proof for Theorem 2.1
For the proof sketch, see the theorem 2.1 in the main text. Here, we first provide the lemmas and corollaries needed for the proof and then provide the proof itself. Lemma C.1 rewrites the SNR of the naive Monte Carlo estimator from eq. (2) as function of moments of the estimators. Lemma C.3 then provides the general expression for the moments of the naive Monte Carlo estimator. Combining lemmas C.1 and C.3, gives the SNR in the form normalizing constants in eq. (4). Lemma C.4 gives the KI divergence between any two posterior distributions. Lemma C.5 then uses the result in lemma C.4 to give the expression for the sum of the KL divergences between three posterior distributions, where the first distribution is the same between the two KL terms (this is the general structure that appears in the eq. (3)). Corollary C.6 simplifies this sum KL expression for datasets with a particular structure. Fimally, we prove the theorem 2.1 by combining all the above results.
Lemma C.I. Let R& be the Monte Carlo estimator in (2). Then,
VK
[B[R]]
SNR (Rk)=
where log
(21)
exp (8)²
E RR
Proof. The proof follows naturally from the definition of SNR (Rx).
EE [[
SNR (Rx) VKSNR (Ri) VK-
(22)
V[R1
[R
KK
(23)
E
RH
E
K
60
(24)
3|A
VK
VK
(b)
(25)
/ exx (26)
v exp (8)²
-
(a) where (a) follows from the fact LHS and RHS of are equal for E[R,] and limit is the same at E[R]] = 0; and (b)
follows from the definition of
Definition C.2 noo-ooomalii function). Let D some dataset. Let p(D)2) be the likelihood and p(z) be the prior. p(D) {))p)) cxp V(D) where
Then, posterior distribution p(=|D)
V(D) = log ((D)=))(2)d..
(26)
Lemma C.3. Let p(D]z) be the likelihood and p(z) be the prior Let DD be some test data. Let p(D + D*1z) = p(D)=)p(DD1)) for any and D*. Let Rí be the Monte Carlo estimator for the PPD, under exact inference (eq. (2) with K and (D p(2(D)). Then.
exp V(D cD*)
[
(27)
exp V(D)
where e a non-negee non-negative integer and V is as definition C.2.
15

---

Understanding the difficulties of posterior predictive estimation
Proof. The proof is straightforward for 0. For 1, we have
E
R
(D( (b)
E [p(cD*)==)
(28)
f (2cD"\2)p(=(DDP)
(29)
[p((D=)))p(D) (s) =)p())))
(30)
exp V(D)
( fp(D + cD* 2))))))(2)d=
(31)
exp ((D)
@ exp V(D + cD*)
(32)
exp V(D)
where (a) follows from definition of (2), (b) and (c) follow from the i.i.d assumption on the datasets, and (d) follows from definition C.2. Note: we do not require points within a dataset to be i.i.d.
Lemma C.4. Les p(D)2), p(z), and p(z|D) be as in definition C.2. Let Da and Dr be the two multisets of data. Then,
KL (p(2) |l p(2)Db))
(33)
= E log (Dal=))
V(D.) + V(D,)
(34)
p(Db22)
(35)
Proof.
KL (p(2|D.) p(=(Da))
(36)
[(2] Do)
IE
log
(37)
p( = (Db)
p(D =)p(z) expl (Da))
EE
log
(38)
(2)
exp (Da)
p(Dul=)
E
log
-V(D)) V(D,)
(39)
p(D))z)
(40)
Lemma C.5. Let p(D)2), p(z), and p(z(D) be in definition C.2. Let D1, D2, and D3 be the three multisets of data. Then,
1 - (I Kl. (p(=)Da) p(=())) + E [oggp(DD) logp(D) log (22=))
p(=(D2))
(41)
(42)
V(Di) + V(22) V(Da).
(43)
Proof. Applying lemma C.4 to D3 and Dj gives
KL (p(2 Dt))
(44)
E [logp(Ds|z)
-
(45)
(46)
16

---

Understanding the difficulties of posterior predictive estimation
and applying it to D3 and D2 gives
KL (D(2//a) ((2)
(47)
E (log((Da)>) - log po V(Ds) + V(D2).
(48)
(49)
Now, multiplying the above two equations by 1 and adding them gives
2
KL pp=/Da)) Di)) KL (p(z 11 p(z(Dz))
(50)
1 E [ogg((D322 log (1 P1]2]]
1 V(Dz) + 2 (DD)
(51)
2
2
1
E (logp(Ds(z) - log pl (D2lz)] V(Ds) + 2 V(D2)
(52)
2
E logr(Dslzz)
log (D) log (D2l(D22=)
(53)
2
V(Di) + V(D2) V(Ds).
(54)
2
Corollary C.6. Let p(D)2), plz). and D(z(D) be in definition C.2. Let D1 = C.D D2 c.D + 20, and D3 c.D + D² be the three multisets of data where and o are non-negative integers. Then,
KL (p(z|Ds) II + K (I(taa)
(55)
V(Di)+ V(Dz).
(56)
2
Theorem C.7 (Repeated for convenience). Let (D(z)) be the likelihood and p(z) be the prior. Let D* be some test data. Let p(D + D*12) p(D]=)()))))))) >)p(D*(>) (2) for any D and D Let Rx (as in eq. (2)) be the Monte Carlo estimator for the PPD, under exact inference. Then, the signal-to-notes ratio of Rx is given by SNR (Rx) = /K/ /exp(5)2 where
8 N V(D) VVD++> V(D+D")
+ KL (p D*) p(2 2D*))
(57)
(58)
2
where V as definition C.2.
Proof.
(0) 1
[[
log E| 2
(59)
2
1
log E [12]]
log EE [Ra]
(60)
2
exp V(D + 2D*) - -log
exp V(D D*)
log
(61)
exp V(D)
exp V(D)
( V(D + 2D') + V(D) V(D+D')
(62)
2
(a) follows from lemma C.I, (b) follows from lemma C.3, and (c) follows from some simple algebraic manipulations.
Now, for the KL-divergence result, if we take the expression in corollary C.6, and plug D1 D and D3 D + 2D and
D3 = D + D*, then we get the same expression as eq. (58).
17

---

Understanding the difficulties of posterior predictive estimation
D. Discussion for informal result in eq. (5)
In this section, we provide detailed discussion for the informal result in eq. (5). First, we provide a lemma and a corollary that will be used in arriving the main result: Lemma D.I gives the KL divergence between any two Gaussians and Corollary D.2 simplifies the expression for the sum of the KL divergences between three Gaussians, where the first distribution is the same between the two KL terms (this structure appears in the SNR expressions in eq. (3)). Finally, restate the assumptions and calculate the result in eq. (5) using these results.
At the end of the section, we provide a note on the convergence of KL divergence and why it requires additional assumptions than just the convergence of the distributions.
Lemma D.1. Let N ((Ho, Eu) and (11, Σ1) be two Gaussian distributions of dimensionality d with Eo: E1 > 0. Then,
KL (N Mo) ||
=
2 20
1 2
-d
1
2
1
+ (11 - -10) _ (1z - 14)
2
1
+ In det Σn1 In |det Eol
(63)
2
Corollary D.2. Let Nuo, Eo). N(n, E1), and Ni2, En) be three Gaussian distributions of dimensionality d with Lo, E1, and Mg > 0. Then,
(L(N(HO: Σo) [|,, E1))
KKL 2o) [M2, Mn²)
tr ((2 Σ-
((-2)-' + Σイ Σn d
2
1
(pI -Mo) [
1 - -o) + 2 mO) T "E²' 1 Mo)
2
+ In |det Σn + 2 |det Σg| - In ddet Eol 2
For any dataset D, let âp be the maximum likelihood estimate and -Sp ! be the Hessian evaluated at the maximum likelihood estimate ? log ((D)e)), such that,
â = argmax log
and Sp¹ Vo (D)
(64)
The idea behind eq. (3) is to use the Berstein-von Misses theorem, also known as the Bayesian central limit theorem. While there are several variants of this theorem, Theorem 10.1 from Vaart (1998) suffices for our usecase. Note that Theorem 10.1 states the result in terms of the true parameters. For the extension to the maximum-likelihow estimate, see the discussion after Lemma 10.3 (Vaart, 1998).
Under eq. (3), we essentially assume that the posterior p(2|D) ( (z] (EMLE))) where I is the estimate of the Fisher information matrix evaluated at the maximum likelihood estimate. Since the estimate I I= (2MLE) ? (D))mM)))) the [D) cancels out and we get
p(=[D) R N (zláp, p2+|D+D' = A : and p D + 2D*) = N (zlep120= SD420*)
(65)
(66)
(67)
Now, assuming that the posteriors can be approximated by their corresponding Gaussian approximations under Bayesian CLT is not sufficient to guarantee that the KL. divergence between the original posteriors will also converge to the KL divergence between the Gaussian approximations. Convergence of KL divergences requires additional assumptions (to be precise, we need to additionally assume absolute continuity and uniform integrability, see section D.1 for discussion). For
18

---

Understanding the difficulties of posterior predictive estimation
simplicity, we side-step the convergence of KL-ddivergence discussion and directly analyze the expression for the sum of KL-diiergencess between the Bayesian CLT approximations. Using Corollary D.2, we can simplify follows.
(&D+DD (ED)
+ KKL (2D+D-,SD+D.) || M (
-tr ((25)
2 1 SD+D1 -d
1 + (2p 2D+D-) 1S (2D - 2D+D-) ) + (2D42D- - 2D+D-) S (2D42D- - 2D+D.)
2
+ In det Sol + In |dct SD+2D-1 - In |det SD+D-1.
(68)
The above expression is terse and hides the relationship of SNR to the dimensionality of the problem and the relative size of the datasets. To bring out these relationships, we make the simplifying assumption that the datasets D, D + 'D", and D + 2D" are similar enough such that the maximum likelihood estimate and the estimate of the Fisher Information matrix similar under the three datasets. Mathematically, we assume that
2p ~ 2DDD- R 2D42D
(69)
and empirical Fisher Information matrix (or the scaled Hessian) is given by
1 D 2D
D D
D
DDDD
SD420
(70)
Substituting from eqs. (69) and (70) into eq. (68), and simplifying in section D.2, we get
KL (N (ēD4D., ( M (E))
+ KI. (N (BD+D.,SD4D)) V (ED+20;
D D*
dlog
(71)
VIDI D 2D*
Finally, plugging the KL-divergences from eq. (71) into the definition of 8 in eq. (3), we get the result
[ D
dd log 1+ 1 |ID [D]
1
ô dog
(72)
V D 2D
V 2 D D
As we can see, the assumption about the similarity of the empirical Fisher Information matrix and the MLE is necessary to bring out the dependence of SNR on dimensionality and the relative size of datasets. This assumption over the datasets is relatively strong and the main reason we keep the above result informal. Note that the middle term in the above equation shows that is positive- the quantity inside the logarithm larger than one since [D + D*| is the arithmetic mean of [D] and [D 2D* which is always larger than the geometric mean / [D] [ + 2D* The right term clarifies that only the dimensionality and ratio of (D) and [D* matters.
D.1. Note on convergence of KT.-divereence
Although the Bayesian CLT or (Bernstein-von Mises theorem) provides convergence total variation between each posterior distribution and its normal approximation, this convergence alone is not sufficient to guarantee the convergence of the KL-divergence between two sequences of posterior distributions. KL divergence is highly sensitive to differences in how probability mass is allocated, especially in regions where one distribution assigns significant probability and the other does not. This sensitivity means that even with total variation convergence, the KI. divergence between the sequences can diverge or fail to converge to the KL divergence between the limits.
To ensure convergence of the KL divergence, additional assumptions of absolute continuity and uniform integrability are required. Absolute continuity guarantees that wherever one distribution assigns positive probability, the other does
19

---

Understanding the difficulties of posterior predictive estimation
1oo-preventing the KL divergence from becoming infinite due zero probabilities in the denominator of the density ratio. Uniform integrability of the log-density ratios is also required; it controls the contributions the KI. divergence from regions where the density ratios become extreme (either very large or very small). By ensuring that these extreme values do not disproportionately influence the KL divergence, it prevents divergence caused by small probabilities or heavy tails.
Collectively, these assumptions ensure that the differences between the distributions are well-hehaved for the KL. divergence to converge. We skip the formal proof as not the focus of this paper.
D.2. Note for the simplification from eq. (68) to eq. (71)
When the datasets D and 'D have the matching mean statistics, we have the relations in eqs. (69) and (70). Under eq. (69), the quadratic terms in eq. (68) are zero. We can simplify the term involving trace as follows:
((55))
+ 1 2550tz0 SD+D
1
1
2 ir(Spp
SD4D-) ) + 2
ur
2 tr
1
Sp
ID D
SD
1 tr
(1D+
2D*
5T¹
(
(
(D D"I
ST
-1
[D]
2
[D]
[D]
1
D
1
D + 2D*
[
tr
(S Sp)
tr (Sp Sp)
2
[
D*
2
[D]
D + D
1
D
1 DD + 2D'
d +
d
2
D D
2
D + D'
(
[D]
+ 2 D 2D*
d
D D
D D
D D
d,
where (a) follows from the relation in eq. (70) and (b) follows from the multiset notation (Costa, 2021).
Therefore, the first and the second term (d and -)) in eq. (68) cancel out, and the only remaining terms are the ones involving the logarithms the determinants the covariance matrices. These remaining terms can simplified as follows:
1
2 In det (Sp) + In det (SD+2D-) In (8D+D-)
2
n
1
[D]
[D]
In det (Sp) +
In det
SD
In det
SD
2
2
D + 2D-
D + D=
(d)
1
det (Sp) +
1
det (Sp)
d
log
DD]
In det (Sp) dlog
D [ + D+
2
2
2
D 2D
d
[D]
(D)
log
d log
2
D DD22D
2D-
DD
+ D*
@ d
log
D*
1
log
[D]
1
log [D +DD
2
DD1
dlog
[ D + 2D=
where (f) follows from eq. (70); (d) follows from log det (aX) dlog + log det X for any non-negative scalar a; this gives the final result in eq. (71); and (c) follows from simple algebraic manipulations.
20

---

Understanding the difficulties of posterior predictive estimation
E. Proof for Theorem 2.2
This section presents a proof for theorem 2.2 about the SNR the linear regression model. First, we formally define the Bayesian linear regression models in definition E.1. We then restate the assumptions in theorem 2.2 for convenience. Lemma E.2 presents the expression for the posterior distribution where we consider both training and test data under the assumption B1. Lemma E.3 presents the KL divergence between two posterior distributions with test data under the assumptions BI 1 and B2. Finally, we use the KL expression from lemmas E.2 and E.3 to prove theorem 2.2.
Definition E.1 (Bayesian Linear Regression Model). Consider the linear regression model with a Gaussian likelihood such that
p(yp = ((yp|Xpz,o'1))
(73)
where y/ E PR is the response vector, Xp E P is feature matrix, and 0² is the variance.
The conjugate prior a Gaussian distribution such that
p(z) (zluo, Eo)
(74)
where Ho is the mean and Eg is the covariance. Then, the posterior distribution is given by p(z|yp) = N (= Ep),
(75)
where
Ep= = (#xxxx+ + xo
and μD= =p=Ep Ep (50 -2Xov0 DS I5'HO)
(76)
Lemma E.2. Let p be the Bayesian linear regression model from definition E. Let BI and B2 hold. Let be non-negafive integer. Then,
= -N(ll==t DD+cD.),
(77)
where as |D] - DO,
1 (xxx)))) and PD+eD- + (XXXD) X yD
c
EDTED 7 c+1
(78)
Proof. We first massage the expressions for the covariance and mean posterior distribution such that we can use the Bl and B2.
ED
( XxXx
1
+ Xxxxx
(79)
(xxx8))
1 (1((XXX)) " (x ) -1
1
(80)
Based on the B1, we have
XD
XD
XD
Xp
DDtcD-XDtoDr =
= =
(81)
XD.
[D
Plugging eq. (81) into eq. (80) gives
D = -(1c+1)xx5X)) ( (((((())) --xxxxx
2:01) )5²))
(82)
21

---

Understanding the difficulties of posterior predictive estimation
similar massaging of expressions for the mean gives us
= =Ep
02 @XXDUD + 20¹1 1O)
(83)
Xxxxx + 1
1
DDDD + 2T²10 14
(84)
(xxx2) 1 1 (XxXD) - E )/'))²
1
0² XDUD Eg mu)
(85)
-e 1 + (XDxx) 2/²1 (Xxxx))
XoXD + 2/ 2 HD)
(86)
-े I (xDXp)
-1 2σ
(xxXp) - XTyD + (XDXD) - 20 mo)
(87)
Now, based on the B1, we have eq. (88).
XD
VD yD + A
XD
X DD+t- =
= xxc( 1)yD + cA)
(88)
[D]
y/D + +
Plugging eq. (81) and eq. (88) into eq. (87) gives
(
1
1 +
((c+ ((++1)x5XX) -1 E x51)
PD+cD*
(ce 1)X 1 XT ((e+ 1)p + cA) + 11)XTX) -2 -507 pm) c
1
(89)
Now, taking limits (as |D| + x: implies (XxXp) -1 D-1 0 from B2) on both sides of eq. (89) leads to simplification of terms inside the brackets and gives us
2n
lim EDTED. |D|
((c+1)XXXX)) ( ((c+1)xXXX) E ))
lim DD|+00
(90)
T xbxxx)
(91)
c+1
Similarly, taking the limits on both sides of eq. (87) leads to simplification of terms inside the brackets (under the obvious assumption that no is finite) and gives us
lim WD-ODD
ID +00
I+ I ((c+ ()xxx)) -1 I 'In))
1
= lim
|D]
0²
((c 1)XXXx))
-1 X ((ct 1lyD + +c) + ((c+ 1)xTxp) -1 "20010 2 mo)
(92)
(XFXx) xo c
c+1 A).
(93)
22

---

Understanding the difficulties of posterior predictive estimation
Lemma E.3. Let p be the Bayesian linear regression model from definition E.I. Let B1 and B2 hold. Let and 3 be two non-negative integers. Then, [D] DC,
KL
11 (D*)) - (ko + A Ma,
(94)
where Ka. is positive constant and Ma,s is a positive definite matrix such that
3++1 g +1
a+1 3+1
(3 a8) (a (8+ XD X
log
and
Ma.3
(95)
Proof. The result follows directly from plugging the mean and the covariance into the expression for KI. divergence between the two Gaussians, taking limits, and following simple algebraic manipulations. We know that
det Σ₂ 1 2 (1,))",,) E1 log
(u2 - 11) E (12 m1))
KL (N ED) (12, Σ2))
(96)
det Σ₁
Collecting the transpose terms, plugging in the covariance expressions from eq. (82) for the distributions p(= + + aDD) and p(z|D+ 3D*), and taking the limits (as [D] + oo implies (XFXD)
25 0 from B2) get
lim tr(Ez )1) 3 + d. PD
(97)
See section E.2.1 for details of the calculations. Collecting the determinant terms, plugging the covariance expressions from eq. (82) for the distributions p(=|D aDD)) and p(=|D + BDD), ) and taking the limits get
det Σn 100000 det Σ
logg a a+1
lim
log
d.
(98)
B 3+1
See section E.2.2 for details of the calculations.
Finally, plugging in the expressions for the covariance and the mean from eqs. (82) and (89) for p(zID + aD*) and p(z|D++ + BD*) in the quadratic term, and taking the limits we get
m1) = (B - (z)² 1 A Xo (xIXx) 'IA.. (a + 1)² (8+1)02
lim (M2 (1) DDi
(99)
See section E,2.3 for details of the calculations. Plugging the aforementioned limit results back into the KL--dvergence expression and identifying the correct coefficients gives the result.
Theorem E.4 (Repeated for convenience). Let p(z) (2lpoo Ea) and p(yp) 2) (WD) o21), where yD R are the responses, X E RiD[xd xd the features, Rdd are the weights, and o² is the known variance. Let be in eq. (3).
BI: The test features X consist of m copies of X and responses yD- consist of m copies of yo + A where A the mismatch vector of size (D).
B2: For sequence of increasingly large training datasets, posterior is dominated by data. Mathematically, as [D] + oc. (X,xx) Exp +
Assume Bl and B2. Then, lim 15 - 6pml = 0 and +00
1 + m
8tm log
/1 + 2m
m
X (XTXx) "XJA. "XJA.
+ 2g 2m + 3m +
(100)
23

---

Understanding the difficulties of posterior predictive estimation
Proof. 8 can be written average of the KL-divergences between the posterions p(=|D + D*) and p(=[D) and between the posteriors p(=|D D*) and p(=|D + 2D*). From the expressions of KI. divergences in lemma E.3, get
1
1
lim 8 [D] 000
2 Okm + A Mm. A) (km,2md + A 2m
(101)
1
1 A (Mmo M Mm,mm) A
(Rm,D + 2m) d + 2
(102)
2
Simplifying the expressions, we get
km, + kmm.2m
(103)
1 + m
+ m
1 + 2m
loog
+
log
+
(104)
1 m
+ 2m
1 + m
1 1 m
+ m
log
+ llog
(105)
+ 2m
+ m)²
llog
(106)
1 + 2m
+ m
2 2og
(107)
V1 2m
(108)
and
Mm,0 + Mm,mm
(109)
m²
1
XD
(x)
[X
X +
m
X (X -1 Xy
(110)
(m -(m+ 1)2 02
(m - 1)² (2m 1) 0²
m²
1
XD (x)
-1 X
1+ 2m +
(111)
(m + 1)2 02
m²
(mm + 1)
X (xfX,)
-1 XT
(112)
(m + 1)2 2m +1 σ²
2m3
X (XTXD) -1 X7
(113)
(m 1) (2m + 1) 0²
(114)
Plugging these back in eq. (102) we get the results.
E.1. Bounding the dúim in theorem 2.2
dum can be bounded by bounding three individual terms. First, log(1 + m)/ /1)/vi1+ and upper-bounded by d/4 log(m/2- + 1). Second, m²/ (2m² + 3m 1) is lower-bounded by 1/6 and upper- bounded by 1 1/2. Third, we have
1 2m is lower-bounded by d/4 log(m/2)
A XD A (X ' A
(115)
where U the left singular matrix of Xp containing singular left vectors. Then, from the propertics of the sint-sin vectors, [UT AIl terms power-bounded by 0 and upper- bounded by |||||. Combining these bounds, we get the bounds in eq. (7).
E.2. Note for the calculation of the KL terms in lemma E.3
We will consider the calculations for the three terms KI. divergence between two Gaussian distributions.
24

---

d. Understanding the difficulties of posterior predictive estimation
E. TRANSPOSE TERMS
The first term we focus on is the one involving transpose.
tr 2 En)
(116)
tr ED+aD.)
(117)
- -tr ((++))x3x))
20 1 ([ (()+))))) 2 -1= E 1 02 1 ((3 + Xp) - Σp
1
((a+1)x5xx))
1 (2 (la 1)xxXx) -
(118)
((++ 1)x
ta
1 42 ((a + Xb) -1 Σn )
(xxX))
(119)
a+1
tr ++
02 + - ((8 + 11xXfX1) -1 2
1
1+ + ((a 1)x - En
(120)
a+1
1+
ir
- 1 (18 ((++))XXX+) E
1
1 + ((a - 1)XxXx) D
(121)
a +1
(122)
On taking limits 'D] + oo (which is equivalent (XXXD) -1 En1 some simple manipulations gives
0 from B2) for the above expression and doing
lim tr(2z [n)
12 +00
B+1
8+1 a+1
+1 tr
02
=
(123)
E.2.2. DETERMINANT TERMS
The second term we consider is the one involving determinants.
log det Σe det Σ
(124)
log det Eg log det. Eg
(125)
log det ((3 1)X
((8+((1))XXx)) -
2:51'
(126)
1 (1 ((a 1)xxxp) "' xn!
+ log ((a + 1)XXXp)
(127)
ddoog + 1) - logdet det X log det
1 + ((B 518+(1+
1)xx)p)) -' 2000
(128)
+ dlog(a + 1) log det X + log det
( ((a
1
1)XxXX Σ
(129)
++ 1
logdet 1 1 ((8 + 1)x -1 "
-''''
dlog
(130)
+1
+ log det 1 (azl ((a (Ixxx)) -1 2
(131)
25

---

Understanding the difficulties of posterior predictive estimation
On taking limits for the above expression and doing some simple manipulations gives
det Σ2 = d dog
a+1
lim
log
(132)
100000 det Σ1
3+1
E.2.3. QUADRATIC TERMS
Noting that the terms involved are well-behaved functions of mean and covariance, we directly plug the limit expressions for mean and covariance from Lemma E.2.
lim (#2 - m11 Eng ( p1) [o] 1-
(133)
-(xxXx) -1 x (VD + (+1 - (XDXp) -1 X + ~)
B
1 xxxx)
(8+1)
Xb(1D ((xDX)) -
B
A) - (XxXx)) X + a 1
(134)
3+1
(xxxX)) -1 Gx
8
a+1 A
3+1
1
(8+1)
((X5X)) -1 X1
8
at 4)
(135)
B+1
() a)² 1
ND (xTXD) -T (xFX²) (xDXx)) -1 =XAA
(136)
(0+1)(8+11) 0²
(8-a)²
A X (xxX)) X1A
(137)
(a
26

---

Understanding the difficulties of posterior predictive estimation
F. Proof for Corollary 2.3
This section contains the proof for corollary 2.3. The overall idea of the proof is similar to the proof of theorem 2.1 in section C. Lemma F.1 provides the expression for the moments of the naive MC estimator for the conjugate exponential family models. Lemma F.2 provides the expression for the sum of KL divergences between three posterior distributions in the conjugate exponential family models, where the first distribution is the same for the two KL terms. Finally, prove corollary 2.3 using these two results.
Lemma F.1. Let the likelihood p[y]=) be exponential family (eq. (8) and prior p(z) = s(zlEo) be in the corresponding conjugate family (eq. (9)). Let D be multiset of training data, DD a multiset of test data, and let Rí be the Monte Carlo estimator for the PPD with exact inference (eq. (2) with 1). Let h(D") = II h(y). Then,
rED
E =h("") exp B(D + cD*
(138)
exp ((D)
where e is a non-negative integer and B is as in eq. (9).
Proof. Starting from the definition of Rí we have,
]
B = ()))
II p(ulz) VED=
(139)
II h(y) exp (T(v) d(z) A(=))
(140)
VED-
[
(M(D") exp (T(D') ((x) - (D*A(()))
(141)
(142)
where (a) follows from T(D*) = EyED T(y) and h(D") = IL h(y). Doing some basic manipulations, we get
E (M(D") exp
(143)
>(D')(() E exp @ M(D))
(144)
exp (sp(c"D")) dD*IA(2)) (2)
(145)
@ M(D') exp (CT(D") ((=) - <D A(=)) exp (T(D) T - (D)A(=))
(146)
exp ((())
MD'' exp (T +DD) T((2) ( (()))(((( dz
(147)
exp(B(fz))
@ h(D')C N (xp(B((D+.o)))
(148)
exp(B(Ep)
(149)
where (b) and (c) follow from the definition of s(2(ED) (eq. (9)) and the fact that the expectation is under the posterior; (d) follows from the multiset notation (Costa, 2021); (e) follows from the definition B in eq. (9).
Lemma F.2. a canonical exponential family p(a in) Mxx exp (T(x) T A(r)). the looseness of Jensen's equality applied to the nog-partition function A at points v, w, and stw
1
2 (A(v) + A(w)) A(u) = KL (p(xlu) p(z)))) + KL (p(z)u) p(z[w))
2
2
27

---

Understanding the difficulties of posterior predictive estimation
Proof. The KL-divergence between two canonical exponential family distributions with parameters and is given by
p(a/w) p(2(a)) p(z
(p(r//)) [[)(()) = E log
IE (T(x)) T T(x) A(w) + A(v))
(150)
p( (w)
(w v) T
[T(']] A(w) + A(e)
(151)
p(r/1))
(w u) T VA(w) - A(w) + A(e),
(152)
where (a) follows from the definition of the gradient of A.
Now, rearranging terms in eq. (152) gives an expression for log-partition function A at any point w in terms of the Log-partition function A at any other point v and the KL-divergence between the two distributions:
A(w) A(u) + (w v) VA(w) KL (p(x)w) p(x|o))
(153)
Replacing with u in eq. (153), gives
A(u) A(v) (u- w) v((u) KL (p(al) p()o)))
(154)
and replacing w with and with in eq. (153) gives
A(u) A(w) + (u - w) V - KL (p(x)) p(zp)))
(155)
On averaging eq. (154) and eq. (155) the VA(u) terms cancel out, and we get
A(u) 2 (A(v) + A(w))
] KL (p(z)/))))))))))))))))))))())( p((l)))
} KL (p(zlu) plzx)))
(156)
2
Finally, rearranging the terms, proves the result:
1 2 (A(b) + A(w)) A(u)
1 (p(rlu))) p(zl))) KL KL (p(x/w) || p(xlw)))-
(157)
2
Theorem F.3 (Repeated). Take model with a likelihood 2(y)z) in exponential family (eq. (8)) and prior p(2) s(2(Ea) in the corresponding conjugate family (eq. (9)). Let "D" be some test data. Let Rx be the Monte Carlo estimator for the PPD under exact inference (eq. (2)). Then, the signal-to-noise ratio is SNR(Rk) -VK/ V/eees eeepppp/ (0)2 for
KL (s(2 D + D²) s(2|D)) + (s s<
+
(158)
B (ED) + B(ED+2D=) B (ED+D-),
(159)
2
where for any dataset D, ED are the parameters that make the conjugate family s(2[ED) equal to the posterior density P(z(D) (eq. (10)), and as in eq. (9).
Proof. From Lemma C.1 get SNR (Rk)
VK for 8 log((E((R)) [B Using Lemma F.1, for c /exp(6)2-1
B(R
1
and c 2, we can simplify 8 as
E R
8 log
log R log [R.)
(160)
2
R
2
2
log h (D*) exp ((D 2D*)
- log hoggh(D") exp B(D + DD)
(161)
2
cxp B(D)
exp B(D)
b
log exp B(D + 2D*)
log
exp exp B(D + D
(162)
exp ((D)
exp B(D)
B(ED+2D.) + B(ED) B(ED+D+),
(163)
28

---

Understanding the difficulties of posterior predictive estimation
where (a) follows from Lemma F.1 for and 2, (b) follows from cancellations of log ((D)), and (c) follows from simple algebra.
Now, observe B in eq. (9) is the log-pan log-partition function of a canonical exponential family. Using Lemma F.2, and plugging ED, ED+D*, and w CD+2D- for conjugate prior family gives eq. (1
29

---

Understanding the difficulties of posterior predictive estimation
G. Proof for Theorem 3.1
This section contains the proof for theorem 3.1. The idea the proof similar to the proof of theorem 2.1 in section Definition G.1 defines the log constant and the posterior distribution for the augmented posterior with test data. Lemma G.2 provides the expression for the moments of the naive MC estimator under the approximate inference. Lemma G.3 provides the expression for the KL divergence between two augmented posteriors. Lemma G.4 provides the expression for the sum of KL divergences between three augmented posterior distributions, where the first distribution is the same for the two KL terms, Corollary G.5 simplifies the expression for the sum of KL divergences when the datasets have specific structure. Finally, we use these four results to prove theorem 3.1.
Definition G.1. Let (2] be the likelihood and p(2) be the prior distribution. Let qD(2) be the variational distribution. Let 2) be some test data. Then,
and (D(2|D") ((D))))g0(z))
ZD(D") = log p(D*)s)oo<zddz
(164)
(D*)))
Lemma G.2. Let 2() be the likelihood and p(z) be the prior distribution, Let (D(2) be the variational distribution. Let D' be some test data. Let p(D + 1) ((D)))(DD for any datasets and D Let Rx be the Monte Carlo estimator for the PPD, under upproximate inference (eq. (2) with K = 1). Then,
E R exp ZD(cD"),
(165)
where c is a non-negative integer.
Proof. The proof is straightforward for = 0 as (0) = log / (p(2)dz = 0. For c 2 1, we have
E [
E [p(D)
(166)
EE [p(cD*|=)]
(167)
p(cD"(=)gp(())))
(168)
exp (D(cD").
(169)
Lemma G.3. Let p(D)2), p(z), and (z(2) )) be as in definition G.I. Let Da and Dr be the three multisets of data. Then,
KL (L(((==(D.) (D(()D)) (logr(D. logr(Ds(=)) (D.) ZD(Di)
(170)
(171)
Proof.
KI. (go(= (D) y (=
(172)
((D.l=)op (x)
exp (2.)
E
log
(173)
p(Ds (+)
(DD)
p(Dalz)
exp ZD (Da) log exp (D()))
E
log
(174)
p(Dalz)
p(Da (2)
E
log
(D( + (D,)))
(175)
p(Dolz)
(176)
30

---

Understanding the difficulties of posterior predictive estimation
Lemma GA. Let p(])z), p(z), and (z(2) )) be as in definition G.I. Let D1. D2, and D3 be the three mullisets of data. Let
D1 D2, and D2 be the three multisets of data. Then,
1 KL (1D (2) |l (D(2> DD' 1)) + 1 zKL (gn(2 [10(2=()')))
(177)
2
E
log (Dg log p (Dil2) log (21z))
(178)
(Di) + (D2) Zd (Dz))
(179)
Proof. Applying the lemma G.3 to D3 and D gives
KL (gp(>|D3) || (D(2 1))
(180)
E (logp((a)z) - log (D)p( +
(181)
(182)
and applying it to D3 and D2 gives
KL (gp(2 (D(2lD"2))
(183)
E (Dgg((Dzl)) - log (esp(D22))) - (D())) + (D2).
(184)
(185)
Now, multiplying the above two equations by 1 and adding them gives
KL (QD Dz) (0(2(D")))
(186)
2
1
+ KL (gn(2)Ds)) 2))
(187)
2
E logr(Dalz) 1
log logg(D) (D)))) log logp(D2lz)
(188)
ZD(Di) + (D(22) (Dx))).
(189)
2
Corollary G.5. Let pD))2), p(2), and (2(2) be as in definition G.I. Let D1 = Ca D2 = c.D + 20,,D", and D3 c,D + D² be the three multisets of data where ca and o are non-negative integers. Then,
1
1 KL
KL (gD(>)Ds) ll go(2 1))
(190)
2
(D(Di) + (D2) (Da)
(191)
2
Theorem G.6 (Repeated for convenience). Let pl (D)z) be the likelihood and p(z) be the prior distribution. Let qo(() be the variational distribution. Let "D" be some test data. Let p(D + D*1z) = (PD)) for any datasets D and D* Let Rk be the Monte Carlo estimator for the PPD, under approximate inference (eq. (2) with = 1). Then, the signal-to-noise
ratio of Rx is given by SNR (Rx)= = /K/ eexp(0)2 - where
KL (gD(> go(()) KL (LD((2) || 11 (z|2DD))
(192)
2 _o(2D") -
-Zo(") (D")
(193)
where Zo and (>() are as definition G.I.
31

---

Understanding the difficulties of posterior predictive estimation
Proof.
1
log E[
8
(194)
2
E[R 2
1
log WE [4e]] log 13 [1
(195)
2
6 - - -DD(D")
(196)
Where (a) follows from lemma C.1 and (b) follows from lemma G.2. Lastly, plugging Di = and D2 2D* and D3 = D
into corollary G.5 and observing Zp (0) = gives the result.
32

---

Understanding the difficulties of posterior predictive estimation
H. Proof for Corollary 3.2
This section contains the proof for corollary 3.2. The proof follow the same structure as the proof for theorem 2.1 in section C. Lemma H.1 provides the expression for the moments of the naive MC estimator under the approximate inference in the conjugate exponential family models. We use these results to prove corollary 3.2.
Lemma H.1. Let the likelihood =(y]=) be as in eq. (8) and a prior p(z) s(z)Eo) be as in eq. (9). Let QD(2) s(2(n) be in the conjugate family (eq. (9)). Let D be some test data and let R1 be the Monte Carlo estimator for the PPD, under approximate inference (eq. (2) with K = 1). Then,
E Ri] = h(D")* exp [B (17 U(C')) - - (n)),
(197)
T(D)
e. is a non-negative integer, B is as in eq. (9) and (cD) = [D]
for any dataset D.
Proof. Starting from the definition of R2 we have,
E R 1 = E (1((112)) E
II p(ul=))
(198)
VED
-E
hly) exp (T(y) d(2)) -(())
(199)
VED-
(L(DD) exp (T(D") 6(z) - ID*
(200)
where (a) follows from T(D") = EyED T(y) and I(D") = II h(v). Doing some basic manipulations, we get
E (M(DD) exp (T(D*) o(=) ID'1A())))*
(201)
= -(DD"=E exp (CT(D") (b²)) - dD*)A(())]
(202)
c ((TDD))
T (())
"( exp
(203)
AA(=)
Ø(2)
1)
1(D")
exp
U(DDi)
-A(2)
(204)
@ MDD') / exp
2(2)
1)
U(CD") T
s(2(n)dz
(205)
A(=)
x(2)
((())
1)
f Texp
U(ED")
-A(2)
exp 17'
-A(=))
dz
h(DD))
(206)
exp ((n)
6(2)
exp (U(CD*)+ + n))
de
-A(z)
@ h(DD)C
(207)
exp B(r)
N(D') (
(208)
exp(B(n))
h(D")) e U((DD) -(n))
(209)
where (b) just collects the terms in the exponent into single vector; (c) defines U(cD) = [T()] |D]
for any dataset D;
(d) and (c) follows as expectation is under the variational distribution and the definition of conjugate family in eq. (9); (f)
follows from some simple algebra; (g) follows from the definition of B in eq. (9).
Theorem H.2. Take a model with a likelihood p(ylz) in an exponential family (eq. (8) and a prior p(z) = s(2(E0) in the corresponding conjugate family (eq. (9)). Let (z(2) = s(2|n) be an approximate distribution the corresponding
33

---

Understanding the difficulties of posterior predictive estimation
conjugate family (eq. (9) with parameters n, Let D' be multiset of test data and let Bi,g be the Monie Carlo estimator for the 2 (eq. (2) with K- 1). Then, the signal-to-noise ratio is SNR R1,g) Vexp(6y> for
1
KL
+ KL 2
+1 U((D))))
(210)
(m) + B B U
(211)
2
T(())
where B is as in eq. (9) and (cD)
for any dataset DD and non-negative integer
|D]
Proof. From Lemma C.I get SNR (Rk)=
Voxp(8)-1 for d 1 log [2 [R] [R)]))). Then
VR
1 log
log E [R]]] log IE [R]]]
(212)
[
(B(1) + U(2D')) - B(T)) - (B(n+ - B(i))
(213)
2
B(n + )(D))) + B(T) B(n+
(214)
2
where (a) follows from Lemma H.1 for 1 and and cancellations of log h(D") terms and (b) form simple algebraic manipulations.
Now, observe B eq. (9) is the log-partinn pe-partitiio function a canonical exponential family. Using Lemma F.2, and plugging v 7/, = ry (D*), and w = y + (2D*)) for conjugate prior family gives the eq. (15).
34

---

Understanding the difficulties of posterior predictive estimation
I. General experimental details
All our code is implemented JAX (Bradbury et al., 2018) and run on a single NVIDIA A100 GPU. In table 6, we provide the expressions for computation of different metrics from the results in tables to 3 and section 5.4.
Note BBVI. We rely using standard BBVI techniques for most of our experiments. The hope of BBVI is allow practitioners to not worry about designing special approximation families for each model p(D,=) =) (Ranganath et al., 2014; Kucukelbir et al., 2017; Agrawal et al., 2020; Ambrogioni et al., 202 Burroni et al., 2024). Instead, BBVI treats models as black boxes-only requiring access to V log pDD, x) to update the variational parameters using the stochastic gradients of a variational objective (for instance, IW-ELBO). Ongoing research in BBVI focuses on automating other algorithmic choices (Kucukelbir et al., 2017; Agrawal et al., 2020; Ambrogioni et al., 20;;a;b; Burroni et al., 2024). Such optimization schemes greatly improve the applicability of BBVI and come pre-implemented in popular probabilistic programming languages like Pyro (Bingham al., 2019), NumPyro (Phan al., 2019), and Stan (Carpenter et al., 2017). While we implement our own inference schemes for this paper, we expect the results to be similar if we use the aforementioned libraries.
Table 6: Summary of the expressions of metrics and their computations for the table Tables to 4. We report SNR (R) in terms of E [t] and V[[t] and report explicit form in tables and 8. We use S = 1000 for all our experiments. The results are then averaged over ten independent trials to generate mean and standard deviation numbers in tables 1 to 4
Expression Computation
Expression
Computation
Bloo Rx]
Aa ~ QD. 5 2.=1
log &
Σk D(D (2,,e)
SNR (Rk)
E|Rk]
V
π
B Blog RR]]
Pa,k ~ Tw: E
p (.. log x Ek-
SNR (I2
E[
V
(u(aa,))
Table 7: Mean of SNR for different estimators.
Expression
Computation
F[Rk]
Za,k ~ QD:
Ex
Ek-1
(DD
B[R
Ea,k ~ Two 1 15
EK
p a (Ps,))
3
K
k=1
Table 8: Variance of SNR for different estimators.
Expression Computation
[[Rx]
2g.k qD T Ex2l
1 Σx
K
((),)) - E[Rx]
V[RR]
2g,k Two T Σ
1 DK R
E
(u(F,
35

---

Understanding the difficulties of posterior predictive estimation
J. Exponential Family models: Additional Details
ô for Binomial
for Exp
5 for Normal
20
40
14
100.0
20,
10,00 0
30
12
T:D4
1.0
20
10.00
5
-20.0
10
100.0
--100.0~
50
100
0
100
50
100
[21]
[D*
D
Figure 9: contours. Setting exactly the same as Figure 5
SNR(RR) for Binomisl
SNB(R)) for Exp
SNR(B) for Normal
20
100000 20
10
10
10
30
TDD
110
10
10
-10
110
10
5
50
100
0
50
100
50
100
D*
[D*|
D
Figure 10: SNR contours. (Repeated for easier reference). Setting is the same as the Figure 5,
For each of the three models, we fix the number of training data points |D] 100 and number of test data points |D*| 100. Then, to sample the training data such that the mean statistics of the data T(D) = 10, sample from the likelihood distributions by carefully adjusting the parameters. This means, for normal we sample from M(10, 1); for Exp sample from (0.1); and for Binomial we sample from Binomial(100, 0.1).
Then, to sample the test data, we first select the region of high 8 from the Figure 5 and then roughly try to match the target mean statistics by carefully adjusting the parameters. For Normal, we sample from (5, 1) to target T(D*) =e 5; for Exp we sample from Exp(0.025) (0.02)) to target T(D*) x 40; and for Binomial we sample from B1nomial(100, 0.4) to target T(D*) 40. This strategy leads the numbers in table 5. Note, we only use one test and train setting for our experiments. The results reported in tables and 2 are averaged our ten independent estimations for a single data setting.
Table 9: For the three models: Normal, Exp, and Binomial, we identify the exponential family form from Section 2. For likelihood in eq. (8), we identify base measure hly), one-to-one parameter mapping ((z), and log-partition function A(z). Note, the sufficient statistics T(y) y for all models. For the conjugate prior in eq. (9), we identify the log partition function Bg)), where & (Er,E.) Ex)
Model
p(vlz)
h(y)
((2)
A(z)
= B(E)
x² exp(-
Normal N(vlz, o²)
V2
+
2
log 2no²
En
+
$2
Em
Exp
Exp(yl=)
2
log
[(6++1) log
gfa+1
Binomial
(yn(,,,,2)
log 1-2
n log(1 - 2)
log
T(Er 11) (15 [(ng, +2)
36

---

Understanding the difficulties of posterior predictive estimation
8 for Binomial
for Exp
5 for Normal
20
40
14
100.0
20.0
10,00 0
30
12
4.0
10
20
10
10.00
20.0
5
-20
10
100.0
-1000.0~
50
100
50
100
50
100
[D"]
ID*|
D
Figure 11: (Repeated for easier reference). 8 contours. Settings exactly the same as Figure 5.
6 for Binomial
8 for Exp
5 for Normal
20
20.0
40
10.0 0
15
30
1.
12
10
20
0
10
5
10.0
10
10.00
20.0
-20.0
50 [D1]
100
0
50
200
0
50
100
|D*|
D
Figure 12: ö contours. For each model, we first fix the training data set such that T(D) 10 (shown with red dotted line) and [D] 1000. For all the models, increasing the number of training data points results in lower ó for a given test data statistics when compared to Figure 11.
We learn a Gaussian variational approximation for each of the three models from Table 2. For the models with constrained latent variables (Exponential and Binomial), we transform 2 to an unconstrained space and then adjust the logarithm of the determinant of the Jacobian for correct density estimation (please, see (Kucukelbir et al., 2017, Section 2.3) for more details on such transformations). Our variational family has two unconstrained parameters: μ and To ensure positivity of standard deviation, we transform with the soft-plus function.
We consider two options to initialize and a: Laplace's 's approximation and standard Normal. To pick from the two options, we evaluate ELBO using 1000 samples and chose the option with higher ELBO value. For Laplace's approximation, we use JAX's BFGS optimizer (Bradbury et al., 2018) (for each model, BFGS took less than 50 estimations of logp(z, D)).
To learn the variational parameters, we optimize standard ELBO using ADAM (Kingma & Ba, 2015) with learning rate of 0.001 for 10, 000 iterations. For each iteration, we batch of 16 samples for estimating the DReG gradient (Tucker et al., 2019).
We learn a parameterized Gaussian proposal distribution for each of the three models from Tables I and 2. For the models with constrained latent variables (Exponential and Binomial), we transform & to an unconstrained space and then adjust the logarithm of the determinant of the Jacobian for correct density estimation (please, see (Kucukelbir et al., 2017, Section 2.3) for more details on such transformations). Our parameterized proposal distribution has two unconstrained parameters: f and a. To ensure positivity of standard deviation, we transform a with the soft-plus function.
We consider two options to initialize μ and a: Laplace's approximation and standard Normal. To pick from the two options, we evaluate IW-EIBO, using 1000 samples and chose the option with higher IW-ELBO,, value. For Laplace's approximation, we use JAX's BFGS optimizer (Bradbury et al., 2018) (for each model, BFGS took less than 50 estimations of log ((D >)) or pl D (2)).
To learn the proposal parameters, we optimize IW-ELBO using ADAM (Kingma & Ba, 2015) with learning rate of
37

---

Understanding the difficulties of posterior predictive estimation
0.001 for 1000 iterations. For each iteration, estimator for IW-ELBOMs uses M samples. We M 16 for all our experiments. Note, after counting the Laplace's approximation estimations, we use less than 20, 000 estimations of p(D (z)p(z|D) for learning the proposal.
a single sample of the DReG estimator. Note, single sample of DReG
Fatter tails. For the Binomial model, we observe that the estimates for pPD, are higher than the estimates for pPD,- This can be explained from two observations. First, the approximate posterior has fatter right tails than the true posterior, and second, the test data mean lies to the right of the training data mean (see table 5). This means that the approximate posterior places more mass in the region of test data and the p will be higher than ppD,- In fig. 13, we plot the densities for the exact posterior and the learned approximation @D- We also plot inset- zoomed-in version to highlight the fatter right tail of the approximate posterior. Remember, the variational approximation in the constrained space is obtained after transforming the unconstrained Gaussian variational approximation.
p
g(=
10
5
0
0.0
0.2
0.4
Figure 13: Fatter right tail of qp for Binomial model.
J.1. Empirical Validation for eq. (5)
We consider a model similar to the Normal model where likelihood p(ylz) is given by a multivariate normal N (ylz, Σ) with unknown mean E Rd and known variance Σ- = Lg. A multivariate Normal prior (=200, Ix) gives a conjugate model as section 2. For this model, we vary the number of latent dimensions E {1, 10, 100, 10000, 10000). For each d, we create a training data set with 1000 data points, and set test data DD D, that is, the mean statistics for training and test data sets match exactly. In fig. 14, we plot the ö from the approximation in eq. (5) (shown in blue dotted lines with crosses), and compare it against the from exact calculations in eq. (3) (shown in red solid lines with dots). The approximation is accurate for all d, and 8 scales linearly as predicted. This means for higher dimensional latent spaces, we can have extremely low SNR (Rt) even test data statisties match exactly the training data statistics.
10a
Exact
*>... Appros.
10"
101
109
10-1
10"
10'
10
Figure 14: 8 from approximation in eq. (5) (blue dotted line) accurate when compared to from exact expression in eq. (3) (red solid lines). Also, 5 scales linearly with d (eq. (5)).
J.2. Effect increasing the number of training data points
In fig. 12, we consider the effect increasing the number training data points from [D] 100 (fig. 11), |D] 1000 while holding the mean training statistics, T(D) 10, the same. As the number of training data points increases, 5 gets smaller for any given test setting. To understand why, note 8 as in eq. (3) involves two KL divergences: one between posteriors p(z + D)) and p(=[D), and the other between posteriors p( [D + D*) and p(=[D= 2D*). Intuitively, the number of training data points increases, we either require more test data or bigger mismatch between test data and training data for the two KL divergences to be large. Thus, for any given test data setting, we expect o to be smaller as the number of training data points increases.
38

---

Understanding the difficulties of posterior predictive estimation
K. Linear Regression: Additional Details
K.1. Experimental Details
We consider the exact inference settings and start with a baseline scenario where none of the three factors influencing SNR are too high. Thereafter, we independently increase the three factors: mismatch, the dimensionality of the latent space, and the size of the test data to create three additional scenarios. We the standard normal prior and likelihood with o²
Baseline. We set the number of training data points to 1000, the dimensionality of latent space d = 10, and the number of 1. We then forward sample training data set D and then generate the mismatched data Da by
mismatched copies adding mismatch vector A = 2 to the response vector ypp
More mismatch. We keep the training data same as in the baseline scenario and increase the mismatch vector to = 10.
More test data. We keep the training data same as in the baseline scenario and increase the number of mismatched copies m 10.
More dimensions. We keep the number of training data points, the number of mismatched copies, and the mismatch vector same in the baseline scenario and increase the dimensionality of the latent space to d = 100. We forward sample the training data set D and then generate the mismatched data Da by adding a mismatch vector A = 2 to the response vector yp.
Figure reports the results from estimating PPD, using simple MC estimator Rx from eq. (2) for K 100, 10 10% The error bands the five and ninety-five percentile intervals based on 1000 independent estimations.
For LIS, learn full-rank Gaussian proposal distribution by optimizing the IW-ELBO from eq. (18) with M 16 using the DReG estimator and ADAM optimizer with learning rate of 0.001 for 1000 iterations. We consider different initialization techniques for the variational parameters: Laplace's ' approximation and standard Normal, and pick the one that provides higher initial ELBO. For each optimization step. we use 8 copies to average the IW-ELBO gradient. For LIS, we learn the proposal once, and do 1, 000 independent estimations to estimate the error bands.
39

---

Understanding the difficulties of posterior predictive estimation
L. Logistic Regression: Additional Details
Baseline. We set the number of training data points to 1000, the dimensionality of latent space d = 10, and the number of mismatched copies m 1. We forward sample a training data set D and then generate the mismatched data Da by flipping the first A 0.1 fraction of the response vector yD-
More mismatch. We keep the training data same as in the baseline scenario and increase the mismatch fraction to A = 1.0.
More test data. We keep the training data same as in the baseline scenario and increase the number of mismatched copies to m 10.
More dimensions. We keep the number of training data points, the number of mismatched copies, and the mismatch fraction same as in the baseline scenario and increase the dimensionality of the latent space to d 100. We forward sample the training data set D and then generate the mismatched data DA by flipping the first A 0.1 fraction of the response vector yD-
We learn full-rank Gaussian variational approximation by optimizing the standard ELBO objective using the ADAM optimizer with learning rate of 0.001 for 1000 iterations. We consider different initialization techniques for the variational parameters: Laplace's 's approximation and standard Normal, and pick the one that provides higher initial ELBO. For each optimization step, we use 16 independent copies to average the ELBO gradient.
For LIS, we learn full-rank Gaussian proposal distribution by optimizing the IW-ELBO from eq. (18) with M 16 using the ADAM optimizer with learning rate of 0.001 for 1000 iterations. We consider different initialization techniques for the variational parameters: Laplace's ' approximation and standard Normal, and pick the one that provides higher initial ELBO, We use eight independent copies to average the gradient of IW-ELBO.
40

---

Understanding the difficulties of posterior predictive estimation
M. Hierarchical Model: Additional Details
We use MovieLens25M (Harper Konstan, 2015), a dataset of 25 million ratings for over 60,000 movies, rated by more than 160,000 users, and use movie features made of tag relevance scores collected by Vig al. (2012)).
Movielens25M originally uses 5 point ratings system. To get binary ratings, we map ratings greater than 3 points 1 and less than and equal to 3 to O. We pre-process the data to drop users with more than 1,000 ratings-leaving around 20M ratings. Also, we PCA the movie features to reduce their dimensionality to 10. We used a train-test split such that, for each user, one-tenth of the ratings are in the test set. This gives us se 18M ratings for training (and R 2M ratings for testing). Our of these we randomly select 100 users for experiments.
For Gaussian VI, we use a full-rank Gaussian. We optimize standard ELBO using ADAM for 1000 iterations with step-size of 0.001. For each optimization step, we use 16 copies to average the gradient.
For flow VI, use real-NVP flow with 10 coupling layers for all our experiments. We define each coupling layer to be comprised of two transitions, where a single transition corresponds to affine transformation of one part of the latent variables. For example, if the input variable for the t layer is 2(h), then first transition is defined as
(K)
21d = "I:d
2d+1:D 2d+1:D exp (b)(=())) (g(((1)))
(k)
(215)
where, for the function 8 and I, super-script a denotes first transition and sub-script & denotes the 1tth layer. For the next transition, the 2d-1:D part is kept unchanged and židd is affine transformed similarly to obtain the layer output 2(++11 (this time using s/ (2441:D) and tf (=)+ This is also referred to as the alternating first half binary mask. Both, scale(s) and translation(t) functions of single transition are parameterized by the same fully connected neural network(FNN), More specifically, for first transition in above example, single FNN takes =[x) input and outputs both sf (= and ig (())) Thus, the skeleton of the FNN, in terms of the size of the layers, is [d. H, H, 2(D d)] where, H denotes the size of the two hidden layers (H=32 for all our experiments).
(k)
The hidden layers of FNN use a leaky rectified linear unit with slope 0.01, while the output layer uses hyperbolic tangent for s and remains linear for t. We initialize the parameters of the neural networks from normal distribution [0, 0.001 2). This choice approximates standard normal initialization. We optimize standard ELBO with sticking the landing (STL) (Roeder et al., 2017) gradient using ADAM for 1000 iterations with step-size of 0.001. For each optimization step, we use 16 copies average the gradient.
To learn the proposal distribution for the learned IS estimator, we use real-NVP flow with architecture described above. We initialize it with parameters from the variational distribution. For the Gaussian VI, we fix the base distribution for the flow to the variational distribution. For flow VI, we use the same architecture for the proposal distribution and simply initialize using the parameters of the variational distribution. We optimize IW-ELBO with DReG estimator using ADAM for 100 iterations with siep of 0.001. For each optimization step, we use copies to average the gradient.
Table 10: Wall-elock times for trial of the MovieLens experiments (times were stable across trials with minor variations).
Method
No. of samples (K)
Time for Simple MC (s)
Time for IS (s)
Time for learning the proposal (s)
K 10x
8.82
17.11
33.76
Flow VI
K 106
4487.01
6781.45
28.06
K- = 103
7.15
20.13
34.29
Gaussian VI
K- 106
4094.95
6611 1.28
33.98
Running time. Our method has some overhead to optimize the IW-ELBO (independent of the number of samples K) and some overhead to perform importance weighting on each sample. When the SNR is low, naive Monte Carlo sampling might be faster. However, when a large number of samples is needed, the overhead to optimize the IW-ELBO can be negligible, and the per-sample overhead is small compared the benefits. For instance, table 10 plots the time for MovieLens experiments. If you look at Table 4, switching from naive Monte Carlo from 100 10 samples (with 1, 000 the computational cost) improves SNR less than switching to importance sampling with 10x samples (with around 1.5x x the cost).


Reviews:
[{"review_id": "oQYjkoEjNV", "reviewer": "Reviewer_CNiq", "summary": "This paper provides a theoretical investigation of Monte Carlo estimation of posterior predictive distributions. The signal-to-noise ratio (SNR) of Monte Carlo estimation is shown to decrease under three conditions: increasing mismatch between training and test, increasing dimensionality of latent space, and increasing size of test data. The contributions consist of (1) a theoretical analysis of the SNR in several settings (exact linear regression, exact conjugate exponential family, and approximate inference), and (2) a learned importance sampling procedure designed to improve the quality of the Monte Carlo estimate.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their feedback. We will address the minor concerns and offer some specific comments below. \n\n>Why is SNR a good metric for analyzing the quality of Monte Carlo estimation of the PPD?\n\nWe study SNR because it is equivalent to relative variance and bakes in the idea of how large the estimator variance is relative to the target quantity. The idea of relative variance becomes crucial when the target quantity is numerically small, as in the case of $\\textrm{PPD}_q$ values. To make this precise, let’s consider an example where the $\\log \\textrm{PPD}_q= -100$ (for reference, all estimates of $\\log \\textrm{PPD}_q$ in Tables 1, 2, 3, and 4 are lower than $-100$). Also, consider an unbiased estimator $R$ with variance $\\mathbb V(R)=10^{-20}$. In the absolute sense, the variance of this estimator is low; however, $R$ carries more noise than signal. To see why, note that $\\textrm{PPD}_q=\\exp(-100) \\approx 3.72\\times10^{-44}.$ Intuitively, $R$ varies on the scale of $10^{-10}$ (standard deviation) and will produce noisy approximations of the target value that is the order of $10^{-44}$. SNR naturally captures this intuition: $\\textrm{SNR}(R) = \\textrm{PPD}_q/\\sqrt{\\mathbb V [R]} \\approx 3.72\\times10^{-34}$ and flags the estimator as poor.\n\nSeveral works before us have studied the SNR of unbiased estimators of numerically small quantities (as cited in lines 59-61). For the works cited in the main paper, the focus is to study the behavior of gradient estimators since the value of the gradients can become very small as the optimization proceeds. In our extended search, we found two more works that study the SNR of gradient estimators for VI [1, 2] and one more work that studies SNR for the policy gradients in RL [3] (the definition of the SNR in [3] is different from our ratio of moments definition; however, the idea to study the SNR is similarly motivated). In hindsight, we agree this discussion should be more prominent. We will add these citations and the above discussion to the already existing discussion in lines 59-61.\n\n[1] Liévin, Valentin, et al. \"Optimal variance control of the score-function gradient estimator for importance-weighted bounds.\" NeurIPS, 2020.\n\n[2] Rudner, Tim GJ, et al. \"On signal-to-noise ratio issues in variational inference for deep Gaussian processes.\" ICML, 2021.\n\n[3] Roberts, John, and Russ Tedrake. \"Signal-to-noise ratio analysis of policy gradient algorithms.\" NeurIPS, 2008.\n\n>How do the theoretical results here influence the way that practitioners estimate PPDs or the way they select approximate posterior distributions?\n\nWe believe our theoretical results are of independent interest and also provide useful insights for practitioners. First, our work serves as a strong cautionary note. We hope that practitioners monitor (and report) the SNR value of their PPD estimates alongside the $\\log \\textrm{PPD}$ value. Second, whenever the SNR values of the naive estimator are low, our work shows that this does not necessarily reflect on the accuracy of the approximate posterior but simply on our ability to accurately estimate PPD. In cases of low SNR, we suggest practitioners use approaches like LIS to improve the reliability of the estimates and obtain a clearer picture of the relative performance of the different approximate posterior methods. \n\nWe briefly touch on the above discussion in the conclusion paragraph in Section 6 and will expand on this in the update. \n\n> Related work is contained only in the appendix, whereas model details are placed in the main paper after the conclusion.\n\nDue to space limitations, we moved the related work section to the appendix. In the camera-ready version, one more page is allowed, and we plan to use it to bring more of the discussion on the related works into the main paper (note that a small paragraph is present in Section 6). If the reviewers believe that the paper will benefit from moving *all* of the related work section to the main paper, then we propose moving the Model Details section to the appendix to create more space. \n\n> The related work section goes back only a decade or so, which makes it difficult to assess the significance of the results in this paper.\n\nAny research review can fall short of being exhaustive. We believe we made an honest attempt to cover the relevant literature. We cite research on PPD evaluation in Section B, and we will further expand on our motivation to study SNR (as mentioned in this rebuttal). If we are missing some relevant work, we will be more than happy to include it. However, we respectfully disagree that it is difficult to judge our work. We will be obliged if the reviewer can point us to specific literature that we have missed, and without which it becomes difficult to assess the significance of our analysis.\n\nOverall, we thank the reviewer for their time and hope that they will reconsider their position in light of the explanations and other reviews.", "claims_and_evidence": "The theoretical claims are supported by proofs in Appendices C - H and Tables 1 - 4 show the improvements due to the learned importance sampling.", "methods_and_evaluation_criteria": "The experimental setup seems reasonable, consisting of exponential family, linear regression, logistic regression, and a hierarchical model on the MovieLens data. Evaluation criteria include lower bounds on the log posterior predictive density and the signal-to-noise ratio, which make sense for this work.", "theoretical_claims": "I did not examine the proofs in detail. The results seem to intuitively correspond to expressions from information geometry.", "experimental_designs_or_analyses": "The experiments seem to be reasonable but I did not verify them in detail.", "relation_to_broader_scientific_literature": "The results should be better connected to the literature of analyzing the variance of Monte Carlo estimators. In this paper, the signal-to-noise ratio is the object of concern, and to my knowledge this has not been analyzed. The related work section goes back only a decade or so which makes it difficult to assess the significance of the results in this paper.", "essential_references_not_discussed": "References discussing e.g. Monte Carlo Variance and variance reduction techniques.", "other_strengths_and_weaknesses": "The structure of the paper is nonstandard which makes it more difficult to read. There is no background section building up to the results and justifying the use of SNR to analyze PPD estimation. Related work is contained only in the appendix whereas model details are placed in the main paper after the conclusion. More connections to the literature justifying the use of SNR and arguing for the significance/impact of these results are needed.", "other_comments_or_suggestions": "p. 2, l. 57: ratio", "questions_for_authors": "1. Why is SNR a good metric for analyzing the quality of Monte Carlo estimation of the PPD?\n2. How do the theoretical results here influence the way that practitioners estimate PPDs or the way they select approximate posterior distributions?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "STsQwOgBCc", "reviewer": "Reviewer_LsR8", "summary": "This paper provides a theoretical framework explaining the severe signal-to-noise ratio (SNR) degradation observed in naive posterior predictive distribution (PPD) estimators. \nIt rigorously demonstrate that, even with exact inference, SNR diminishes with increasing: (1) training-test data mismatch, (2) latent space dimensionality, and (3) test data size relative to training data. \nThese theoretical findings are empirically validated through numerical experiments on both synthetic and real-world datasets.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your thoughtful and constructive review. We appreciate your recognition of the theoretical contributions and the thoroughness of our analysis. Your feedback is invaluable and greatly encouraging. \n\nFor the camera-ready version, we will plan to use the extra space to add more commentary on the takeaways from sections 2.2 and 3 (as pointed out in [response to the reviewer hMnG](https://openreview.net/forum?id=TzfGuKazvf&noteId=ywqVSRwl8H)). We will also use the extra space to move more of the related works discussion to the main text (as indicated in [response to the reviewer CNiq](https://openreview.net/forum?id=TzfGuKazvf&noteId=hM1fsV2z27)). In fact, we are open to moving the model details to the appendix and, instead, moving all of the related works section into the main text. Please let us know if you have a strong preference for one or the other.", "claims_and_evidence": "The paper's claims are rigorously supported by comprehensive references and robust mathematical proofs, built upon well-justified assumptions", "methods_and_evaluation_criteria": "This paper analyzes scenarios involving exact posterior distributions, represented as $Q_{D}(z) = P(Z|D)$, as well as cases where  $P(Z|D)$ is approximated using variational inference (VI) or Laplace approximation. This comprehensive analysis provides a complete theoretical framework. \n\nThe theory's reliability is further strengthened by its validation on both synthetic and real-world datasets.", "theoretical_claims": "The proofs presented herein are, to the best of my ability, accurate.", "experimental_designs_or_analyses": "The numerical experiments provide clear empirical validation for all theoretical claims presented. A natural and compelling extension would be to investigate whether these theoretical findings hold within the context of Bayesian neural networks, a research area of considerable interest to the contemporary machine learning community.", "relation_to_broader_scientific_literature": "This paper identifies three scenarios leading to rapid signal-to-noise ratio (SNR) decay in posterior predictive distributions (PPDs) and proposes adaptive importance sampling as a mitigation strategy, providing valuable insights into PPD failure modes for the machine learning community.", "essential_references_not_discussed": "N/A.", "other_strengths_and_weaknesses": "N/A.", "other_comments_or_suggestions": "While the appendix provides a comprehensive review of related work, integrating key discussions into the main text would enhance understanding for readers unfamiliar with the field's cutting-edge developments.", "questions_for_authors": "N/A.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kmSLpwL5Yq", "reviewer": "Reviewer_qCWu", "summary": "This paper addresses the issue of unreliable posterior predictive density (PPD) estimates when using a simple Monte Carlo (MC) approach, highlighting the previously under-recognized issue of low signal-to-noise ratio (SNR). The authors theoretically analyze and empirically demonstrate that the SNR for posterior predictive estimation significantly decreases under three conditions: increasing mismatch between training and test data, increasing dimensionality of the latent variable space, and an increasing relative size of the test dataset. They propose an adaptive importance sampling method—learned IS—based on maximizing a variational proxy for SNR, to mitigate this issue, showing substantial improvement in predictive estimation accuracy.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your detailed and encouraging review. We appreciate your recognition of our theoretical analysis and the empirical validation of our proposed LIS method, as well as your insights into potential computational challenges. Your constructive feedback is truly invaluable.", "claims_and_evidence": "The main claim in this work is the signal-to-noise ratio (SNR) of the simple MC estimator can sometimes be extremely low, leading to unreliable estimates. To demonstrate this, the authors provide a theoretical analysis in section 2. They provide an analytical expression for SNR in Theorem 2.1 and they demonstrate it empirically on several models: linear regression, logistic regression, and hierarchical models.", "methods_and_evaluation_criteria": "The methods proposed in this paper are appropriate for the problem at hand and follow obviously from the analysis presented. Specifically, the authors identify the central issue of low signal-to-noise ratio (SNR) in simple Monte Carlo estimations of posterior predictive densities (PPDs). To mitigate this, they propose an adaptive importance sampling method, named \"Learned Importance Sampling\" (LIS). In this approach, they optimize a variational proxy to the SNR by maximizing an importance-weighted evidence lower-bound (IW-ELBO), to learn a more efficient proposal distribution.\n\nThe methodological choice is well-justified, as direct optimization of the optimal proposal distribution (which would maximize SNR directly) is often intractable. Hence, optimizing a proxy measure like the IW-ELBO, which has clear theoretical connections to minimizing estimator variance, is a reasonable and effective strategy (Section 4)​\n.\nThe authors also empirically demonstrate that this method significantly improves estimator accuracy across various scenarios, including exponential family models, linear regression, logistic regression, and hierarchical models, providing robust evidence that their methodological choice is practically effective and relevant for the stated problem.", "theoretical_claims": "I reviewed the correctness of the theoretical derivations provided, particularly Theorems 2.1 and 2.2. The proofs provided appear to be sound and accurately support the claims regarding the decay of SNR under specified conditions.", "experimental_designs_or_analyses": "The soundness of the experimental designs was checked, particularly the experiments with linear regression, logistic regression, exponential family models, and hierarchical models. The setups appear valid, and the experimental results strongly support the theoretical analyses. The use of illustrative examples clearly demonstrates the issues with low SNR and the effectiveness of the learned IS method.", "relation_to_broader_scientific_literature": "The paper's contributions lie within the broader literature of approximate inference and Monte Carlo estimation. It builds upon work on importance sampling and also related to prior research on variational inference, predictive uncertainty, and Bayesian model evaluation.", "essential_references_not_discussed": "The authors appear to have cited the relevant works adequately. No essential references seem missing.", "other_strengths_and_weaknesses": "Strengths:\n1. Clearly identifies and thoroughly analyzes a subtle but critical issue.\n2. Provides theoretical grounding and extensive empirical validation.\n3. Introduces a simple yet effective methodological improvement.\n\nWeaknesses:\nComputational complexity and scalability of the learned IS method, particularly in very high-dimensional or large-scale applications, may be a concern but it is less important since the contributions are primarily theoretical.", "other_comments_or_suggestions": "No", "questions_for_authors": "No", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QFB3qtsn2S", "reviewer": "Reviewer_hMnG", "summary": "The authors study the problem of Monte Carlo estimation of the density of the posterior predictive distribution for (approximate) Bayesian inference. They show that simple Monte Carlo estimation can have a low signal to noise ratio (SNR) if the training data and test data are substantially different, the dimension of the problem is high or the test set is much larger than the training set. They propose an importance sampling approach based on learning a sampling distribution using an Importance Weighted Evidence lower bound, and provide empirical support for this method at least partially mitigating the observed problems with low signal to noise ratio when using simple Monte Carlo estimation to estimate the (predictive) posterior density.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their detailed comments. We will address the typos, equation formatting, and other minor suggestions as is and offer comments to other concerns below. \n\n>... Bayesian CLT approximations ... provide additional detail ...\n\nWhile there are several versions of Bernstein-von Mises, we believe Theorem 10.1 from Vaart (1998) suffices for our use case. Note that Theorem 10.1 states the result in terms of the true parameters. For the extension to the maximum likelihood estimate, see the discussion between Lemma 10.3 and Lemma 10.4.\n\nIn short, we assume that the posterior $p(z \\vert \\mathcal D) \\approx \\mathcal N(z_{\\textrm{MLE}}, \\frac{1}{|\\mathcal{D}|} I^{-1}(z_{\\textrm{MLE}})),$ where $I$  is the estimate of the Fisher information matrix evaluated at the maximum likelihood estimate; such that,  $I^{-1} (z_{\\textrm{MLE}}) = -(\\frac{1}{|\\mathcal D|} \\nabla^2_z \\log p(\\mathcal D \\vert z))^{-1}$. Since the size of the dataset cancels out, we get the expressions in equations 65, 66, and 67. In the revision, we will add the discussion on Bernstein-von Mises and the above clarification after Corollary D.2.\n\n[1] Vaart, A. W. van der. \"Bayes Procedures.\" _Asymptotic Statistics_. Cambridge: Cambridge University Press, 1998. \n\n>Assumptions on the data-generating process for ...  equation 70.\n\nIn spirit, this assumption is not that different from assuming the data-generating distributions are the same. The quantity in eq. 70 ($\\frac{1}{|\\mathcal D|}S_{\\mathcal D}^{-1}$) is the estimate of the Fisher information matrix (FIM) under different datasets. We essentially assume that the train and test data sets are similar enough such that estimates (MLE and FIM) under these datasets will be the same. As the reviewer notices, it is tempting to only assume that the train and test data distributions are similar. However, even after assuming the same distribution, we will require additional assumptions for a result in terms of estimates under two given datasets. This assumption about the estimates is precisely the reason we (very explicitly) keep the result informal. We will add this discussion after eq. 70, and if the reviewers suggest, we will reword the assumption in terms of FIM.\n\n>What is being defined in definition C.2? \n\nThere is no specific assumption on the form of the posterior. We are essentially writing Bayes' rule. However, we define $V$ explicitly as we use it as a function of data $\\mathcal D$ in later expressions. We made an organizational choice to aid the reader, but we can move it to an inline definition if preferred. \n\n>The discussion of approximate inference ... but didn't provide insight ...\n\nThe main insight is that for approximate posteriors that closely resemble the true posterior, we expect the relationships from Section 2 to hold as is. However, for arbitrary distributions,  it is hard to make a precise statement beyond \"SNR depends on how much the *posterior* $q_\\mathcal{D}(z|\\mathcal{D}^∗)$ varies from the *prior* $q_\\mathcal{D}(z)$.\" Intuitively, the training data only enters these equations through the conditioning, or, equivalently, from the approximate posterior $q_\\mathcal D(z)$. Since arbitrary approximate posteriors may have arbitrary relationships with training data, it is not immediately obvious how to make a more precise statements. We hope that future research can explore this further. We will add this discussion at the end of Section 3 and in the limitations paragraph in Section 6.\n\n>What is the main (practical) insight I am supposed to take from section 2.2?\n\nThe analysis in Section 2.2 is complementary to the main claims of the paper. From the earlier parts in Section 2, we know that when the training and test datasets match and are large enough, the posteriors in eq. 3 are similar, and the SNR is high. The same is implied by the two conditions at the end of section 2.2. When $\\xi_\\mathcal D$ is large, the training dataset will be large (note that $\\xi_\\mathcal D$ includes both the statistics and the number of data points, see eq. 10). When the points $\\xi_\\mathcal D$ and $\\xi_{\\mathcal D + \\mathcal D^*}$ lie close to the ray emanating from the origin, the datasets have to be similar. Combined together, these two conditions also support the view that when the datasets are large and similar, we expect the SNR to be high. We favored the illustration in Figure 6 in place of words to derive the main takeaway. However, we will add some of this explicit discussion to the main text.\n\nOverall, we thank the reviewer for these insightful comments. We strongly believe that the concerns raised by them will make the paper stronger. We look forward to addressing any remaining concerns. Otherwise, we feel confident in our abilities to add the requested explanatory text, and hope the reviewer will consider a strong recommendation for our work.", "claims_and_evidence": "*Claim: Bias in estimating the log predictive density can lead to unreliable comparisons between methods.*\n\nFigure 1 illustrates the bias in estimating the log predictive density. While bias is clearly illustrated, it would be more compelling if the authors provided an example where the bias leads to issues in comparison between the methods. In the current example, it seems like we would reach the same conclusion, regardless of how many samples are used in estimation. And it seems reasonable to believe this is the correct conclusion (that flow VI is leading to higher predictive density). Also the standard error bars shown are disjoint for all numbers of samples, so it really seems comparison in this case is reliable, even if estimation is poor.\n\n*Claim: SNR is low when training and test data are 'mismatched'*\nFigures 2 and 8 show simulations supporting this claim. The linear regression case also provides support for this claim.\n\n*Claim: SNR decreases as the dimension of the latent space increases*\nFigure 2 and figure 8 (right) shows a simulation supporting this claim. Equation 5 provides a heuristic argument for this claim. I found the sketch of the heuristic argument in the main text incomplete. In particular, it wasn't clear to me what was meant by \"the corresponding Bayesian CLT approximations\". The link to section D provided didn't link to the appendix. When I went to section D of the appendix, I expected to see a clear statement of what is intended by A1 (i.e. what Baysian CLT is being used in this heuristic, along with a citation). But I couldn't find this. Please provide additional detail or point me to where it is already in the text. I strongly suspect some assumptions on the data-generating process for both test and train locations are needed for equation 70 (in particular, something like that all the data is sampled iid from some distribution). I also find assumption A2 difficult to understand at an intuitive level. Could you give concrete examples where it holds/does not hold? Currently think the argument is so heuristic that it does not provide significant insight beyond the simulations.  The linear regression case provides some support for this claim in a particular concrete (although somewhat limited) setting.\n\n*Claim: SNR is low when the test set is large compared to the test set*\nFigure 2 and 8 show simulations supporting this claim. The heuristic support for this claim follows the same argument as for the claim about dimension, and so my earlier comments apply. The linear regression case provides some support for this claim in a particular concrete setting.\n\n*Claim: Adaptive importance sampling can increase SNR relative to simple Monte Carlo*\n\nTables 1,2 and 3 as well as figure 8 support this claim.", "methods_and_evaluation_criteria": "The authors evaluate standard Monte Carlo and the proposed importance sampling in terms of expected log predictive density and estimated signal to noise ratio. While there is not a ground truth expected log predictive density available in many of the experiments, this metric is still useful as 1.) it is commonly used in practice for model comparison between methods, and so the quality of its estimation is of interest. 2.) the direction of bias in estimating this quantity is known, and so higher values must represent lower amounts of bias. The signal to noise ratio also seems a useful metric to report since it is the focus of most of the discussion of the paper.", "theoretical_claims": "I looked over the appendices for the proof of theorems 2.1 and 2.2. I the argument structure makes sense and I did not see errors. Overall the theoretical claims seem believable, although in certain places heuristic arguments made are imprecise (see details in other boxes).", "experimental_designs_or_analyses": "The simulation studies are reasonable illustrations of the arguments made in the paper. I looked over the model details provided in section 7, and did not have significant concerns.", "relation_to_broader_scientific_literature": "To the best of my knowledge, there is not a significant amount of literature studying the quality of Monte Carlo estimation of the (approximate) posterior based on independent samples. The proposed method in the paper builds on ideas from importance weighted variational inference.", "essential_references_not_discussed": "I did not see significant gaps in discussion of related literature (Appendix B). However, I am not familiar enough with several areas related to this paper in order to feel confident nothing has been missed.", "other_strengths_and_weaknesses": "Weakness: Some sections of the text did not have clear takeaways. For example, the discussion of approximate inference (section 3) gave formulas for the SNR, but didn't provide insight into when I should expect the SNR to be high or low with approximate inference. In the discussion for exponential families, there is a condition discussed in terms of statistics lying on a particular ray, but I didn't see how this can be related back to quantities related to the analysis (the available data, the prior and likelihood, the approximate inference performed).", "other_comments_or_suggestions": "The proof sketch for Theorem 2.1 repeatedly refers to steps as *simple*. I don't think this is useful. Either the reader will understand how this step is done and agree it is simple, or they will not in which case telling them it is simple is not helpful (and discouraging).\n\nappendix A.3, typo: “margina” \n\nPunctuation is missing in many equations in the appendix (e.g. eqn 34, 39, 88, 91, 107, 113). Eqns 40, 108 shouldn't be there as it is empty (likely results from an additional new line command).", "questions_for_authors": "What is the main (practical) insight I am supposed to take from section 2.2? I didn't see how the description around \"looseness\" in Jensen's inequality supported the main claims of the paper, nor how to turn it into practical insights about when I would expect the SNR to be low/high.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The paper received positive feedback for its theoretical and empirical analysis with two reviewers recommending acceptance and two leaning toward a weak decision (1 x accept, 1 x reject). While all reviewers acknowledge the novelty and rigor of the theoretical contributions, concerns were raised about novelty, clarity, and the limited discussion of related literature in the main text. The rebuttal addressed several concerns, but some issues remained (in particular the take-away of Section 2.2). All reviewers agree that the paper raises interesting points.
Decision: Accept (poster)