A Heavy-Tailed Algebra for Probabilistic Programming
Authors: Feynman T. Liang
Liam Hodgkinson
Michael W. Mahoney
Venue: NeurIPS 2023 poster (2023)
Keywords: probabilistic programming
static analysis
heavy tails
monte carlo
mcmc
variational inference
TL;DR: Develops an algebra for analyzing tails and operationalizes it in a PPL compiler to automatically inform the design of density approximator with provably calibrated tails.

Abstract: Despite the successes of probabilistic models based on passing noise through neural networks, recent work has identified that such methods often fail to capture tail behavior accurately---unless the tails of the base distribution are appropriately calibrated.  To overcome this deficiency, we propose a systematic approach for analyzing the tails of random variables, and we illustrate how this approach can be used during the static analysis (before drawing samples) pass of a probabilistic programming language (PPL) compiler.  To characterize how the tails change under various operations, we develop an algebra which acts on a three-parameter family of tail asymptotics and which is based on the generalized Gamma distribution.  Our algebraic operations are closed under addition and multiplication; they are capable of distinguishing sub-Gaussians with differing scales; and they handle ratios sufficiently well to reproduce the tails of most important statistical distributions directly from their definitions.  Our empirical results confirm that inference algorithms that leverage our heavy-tailed algebra attain superior performance across a number of density modeling and variational inference (VI) tasks.

A Heavy-Tailed Algebra for Probabilistic Programming
Feynman Liang
Liam Hodgkinson
Department of Statistics University of California, Berkeley feynman@berkeley. edd edu
School of Mathematics and Statistics
University of Melbourne, Australia
Ihodgkinsonuunimelll edu.
Michael W. Mahoney
ICSI, LBNL, and Department of Statistics University of California, Berkeley
mmahoney@stat. berkeley....c.. .dd edu
Abstract
Despite the successes of probabilistic models based on passing noise through neural networks, recent work has identified that such methods often fail to capture tail behavior accurately-unkee
the tails of the base distribution are appropriately calibrated. To overcome this deficiency, we propose systematic approach for analyzing the tails of random variables, and we illustrate how this approach can be used during the static analysis (before drawing samples) pass of probabilistic programming language (PPL) compiler. To characterize how the tails change under family
various operations, we develop an algebra which acts tail asymptotics and which is based on the generalized Gamma distribution. Our algebraic operations closed under addition and multiplication; they capable distinguishing aub-Gaussians with differing scales; and they handle ratios suffi- ciently well to reproduce the tails of most important statistical distributions directly from their definitions. Our empirical results confirm that inference algorithms that leverage heavy-tailed algebra attain superior performance across number of density modeling and variational inference (VI) tasks.
Introduction
Within the context of modern probabilistic programming languages (PPLs), recent developments
in functional programming [51]. programming languages [3), and deep variational inference (VI)
[4] combine to facilitate efficient probabilistic modelling and inference. But despite the broadening
appeal of probabilistic programming, common pitfalls such as mismatched distribution supports [32]
and non-integrable expectations [53] 55 60 remain uncomfortably commonplace and remarkably
challenging to address. In particular, heavy-tailed distributions arise in a wide range of statistical
applications and are known present substantial technical challenges 55 60) Recent innovations
aiming to improve PPLs have automated verification of distribution constraints [32), tamed noisy
gradient estimates [16] well unruly density ratios (53 (5 and approximated high-dimensional
distributions with non-trivial bulks [39]. To address the issue of heavy-tailed targets, approaches
which initialize with Gaussian tails have been proposed 25 [3]. However, these methods
typically require the use of optimization and/or sampling strategies to estimate the tails of the target
distribution. Such strategies are often unstable, fail allow for sufficiently wide array of possible tail bchaviours.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Analyze Target
Calibrate Tails
Refine Bulk
HatNoncal
10"
randoe.variable
def ZG:
recur (11
10-2
Sample
Sarrds
GGA
the www.a..ariillle
KG recor (1
R₂ 2
10-
chas Functional
10
def =2
Samples
x0
Quer
10-2
(1)
(2)
(3)
(4)
Figure 1: Our heavy-tailed algebra ensures that the tails of density estimators and variational approximations are calibrated those of the target distribution. Here, generative model expressed in a PPL (1) is analyzed using the GGA without drawing any samples (2) to compute the tail parameters of the target. representative distribution with calibrated tails is chosen for the initial approximation (3), and a learnable tail-invariant Lipschitz pushforward (see bottom of Table T and Theorem 2 is optimized (4) correct the bulk approximation.
Motivated by this, we introduce the first procedure for static analysis of probabilistic program that automates analysis of target distributions' tails. In addition, we show how tail metadata obtained from this procedure can be leveraged by PPL compilers to generate inference algorithms which mitigate number of pathologies. For example, importance sampling estimators can exhibit infinite variance the tail of the approximating density is lighter than the target; most prominent black-bor -bx VI methods are incapable of changing their tail behaviour from an initial proposal distribution [25 [ and
Monte-Carlo Markov Chain (MCMC) algorithms may also lose ergodicity when the tail of the target density falls outside of a particular family 42). All of these issues could avoided if the tail of the target is known before runtime.
To classify tail asymptotics, we propose three-parameter family of distributions which is closed under most typical operations. This family is based on the generalized Gamma distribution (Equation (2), and it interpolates between established asymptoties on uub-Gaussian random variables [31 and regularly varying random variables [35]. Algebraic operations on random variables can then be lifted computations on the tail parameters. This results in a heavy-tailed algebra that we designate as the generalized Gamma algebra (GGA). Through analyzing operations like X + Y, X², and X/Y at the level of densities (e.g., additive convolution px @ py), ), the tail parameters of target density can be estimated from the parameters of any input distributions using Table T
Operationalizing our GGA, we propose tail inferential static analysis strategy analogous to tradi- tional type inference. GGA tail metadata can be used to diagnose and address tail-related problems in downstream tasks, such as employing preemptively detect unbounded expectations. Here, we consider density estimation and VI where we use the GGA-GGGGpai comppted tail of the target density to calibrate our density approximation. When composed with learnable Lipschitz pushforward map (Section 3.2) the resulting combination is a flexible density approximator with tails provably calibrated to match those of the target.
seemannian-maiiflle methods 17 to sample heavy tails or
Contributions. Here are our main contributions.
We propose the generalized Gamma algebra (GGA) an example of heavy-tailed algebra for probability distributions. This extends prior work on classifying tail asymptotics, while including both Lub-Gaussian/ / sub-exponentials [31] as well as power-law Pareto-based tail indices [[]] Composing operations outlined Table one can compute the GGA tail class for downstream random variables of interest.
We implement GGA as abstract interpretation during the static analysis phase of a PPL compiler. This unlocks the ability to leverage GGA metadata in order better tailor MCMC and algorithms produced by PPL..
Finally, we demonstrate that density estimators which combine our GGA tails with neural networks (autoregressive normalizing flows [39] and neural spline flows (5] simultaneously achieves calibrated tails without sacrificing good bulk approximation.

---

2 The Generalized Gamma Algebra
First, we formulate our heavy-tailed algebra of random variables that is closed under most standard elementary operations (addition, multiplication, powers). The central class of random variables under consideration are those with tails of the form in Definition
Definition 1. A random variable X is said to have generalized Gamma tail if the Lebesgue density of [X] satisfies
Pl (x) xee
as I - oc,
(1)
for some c > 0, R. > 0 and p E R. Denote the set of all such random variables by g.
Consider the following equivalence relation G: x Y if and only if < PiX] (x)p/)())))( +o0 for all sufficiently large z. The resulting equivalence classes can be represented by their corresponding parameters V, p. Hence, we denote the class of random variables X satisfying Equation ① by (v, o, p). In the special case where = 0, for fixed v < -1, each class (1,, a, 0) for > 0 is equivalent, and denoted by R1l representing regularly varying tails. Our algebra operates on these equivalence classes of G, characterizing the change in tail behaviour under various operations. To incorporate tails which lie outside of G. we let R2 incorporate super-heavy tails, which denote random variables with tails heavier than any random variable in G. All operations remain consistent with this notation. Likewise, we let L denote super-light tails, which are treated in our algebra as class where p +oo (effectively constants).
Equation ① and the name of the algebra are derived from the generalized Gamma distribution. Definition 2. Let E IR, σ > O, and E R [R} be such that (v (+++1) 1)/p > 0. A non-negative random variable X generalized Gamma distributed with parameters v, o, P if has Lebesgue density
(x) Co.o, axx
0,
(2)
where Cu,a,p (v- /Tll///////////////// 1)/p) is the normalizing constant.
The importance of the generalized Gamma form arises due to a combination of two factors:
(i) The majority of interesting continuous univariate distributions with infinite support satisfy Equation (T) including Gaussians (v = O, م 2), caammm//cppoeettalccii--qppuu (v > -1,
p == 1), Weibull/Frechet (p is the log-normal distribution (see Example X in Appendix
1), and Student 1'/auchy/Pareto (R,). A notable exception
(ii) The set G is known to be closed under additive convolution, positive powers, and Lipschitz functions. We prove it is closed under multiplicative convolution well. This covers the majority of elementary operations on independent random variables. Reciprocals, exponentials and logarithms comprise the only exceptions, however, will introduce few "tricks" to handle these cases as well.
The full list of operations in GGA is compiled in Table and described in detail in Appendix A GGA classes for common probability distributions are provided in Appendix B All operations in the GGA can be proven to exhibit identical behaviour with their corresponding operations on random variables. with the sole exception of reciprocals (marked by t), where additional assumptions are required. The asymptotics for operations marked with an asterisk are novel this work. For further details, refer to Appendix(M A
Repeated applications. Provided independence holds, composition of operations in the GGA remain consistent unless one applies Lipschitz functions, logarithms, or exponentials. If one of these operations is applied, the tail becomes an upper bound, which remains consistent under addition, multiplication, and powers, but not reciprocals. Given that we are working with fixed class of tails, such behavior is inevitable, and it possible to perform sequence of operations for which the tail no longer becomes accurate.
Posterior distributions. A primary application of PPLs to perform Bayesian inference. To cover this use case, it is necessary to prescribe a procedure to deal with posterior distributions. Consider a setup where a collection of random variables XI,..., X,, are dependent on corresponding latent random elements Zi,... Z, well parameter 0 through functions f, by X; f( Z.). For
3

---

(1), 1 (22,02;;2)
Ordering
lim simsupp++s 222
+00.
(11,01,P1) $ (12, P2)
(02 (1 1.
if # P1.
Addition
ip P2
(v (1 22, -- + t )
ip = P2
Powers
(V,0,p)
( 1, {) for
1-V-2200 0,
if (v + 1)/p > 0 and p otherwise
Reciprocalsi
(v,0,p)-1
R2
Scalar Multiplication
d(,0,p)
= [v,
P1
+
+
o,
if P1, < 0
va
J. 1
if P1, > 0
(1, P2)
Pa
Rypl
if pı < 0, > 0 if fı 0, 0
Multiplication
where μ ==
1p.1
1 + lpzl
|pitles
P1p2l
o o= u() (pil) (oelpa))PPpa)
(n) + 12 P1) (N) 122 01 02.
if < P2
Product of Densities*
(1, (12,02, P2)
if p1 = P2 otherwise.
(v) + v2, 02, p2)
RR+1 if 2 R1
Exponentials#
exp(v.a.p)
otherwise.
[00|ol] 1,1)
v
Logarithms#
log(u,o,p)
=cc
otherwise.
Functions (L-Lipschitz)
S(X
L max X1. Xn)
Table 1: The Generalized Gamma Algebra. Operations on random variables (e.g., X1 + X2) are
viewed as actions density functions (e.g., convolution (MI, 01, P1) @ (12, P2)) and the tail parameters the result are analyzed and reported. In this table, * denotes novel results, and denotes that additional assumptions required.
simplicity, assume that each fi, f..--1....... f..1 where each fij is an elementary operation in Table To estimate the tail behaviour of conditioned on X, we propose an elementary approach involving inverses. For each operation fij. if fij is power, reciprocal, or multiplication operation, let Ríj be given according to the following:
Powers:
fay(2) = 3
Rg (1 11-,1,,0)
Reciprocals: f(((x)= =1
Rg = (2,1,0)
and otherwise, let R.j = 1.
Multiplication:
fu(l,y) = xy.
Rg = (1,

---

Letting f/" (ax,) 2) denote the inverse of fi in the first argument, we show in Appendix A
n.k Rij) 6w/=1
Zi)) 0|x = x= & (a, (i) &
&
&
where и denotes the prior for 0. Since the inverse of a composition of operations is composition of inverses, the tail of oā (x, Z;) can be determined by backpropagating through the computation graph for X; and sequentially applying inverse operations. Consequently, the tail behaviour of the posterior distribution for one parameter can be obtained using single backward pass. Posterior distributions for multiple parameters involve repeating this procedure one parameter at time, with other parameters fixed.
3
Implementation
3.1 Compile-time static analysis
Algorithm 1: GGA tails static analysis
To illustrate an implementation of GGA for static analysis, we sketch the operation of the PPL com- piler high-level and defer the code in Sup- plementary Materials for details. A probabilistic program first inspected using Python's 's built-in ast module and transformed to static single assign- ment (SSA) form 43). Next, standard compiler optimizations (e.g., dead code elimination, con- stant propagation) are applied and an execution of the optimized program is traced 4 58) and accu- mulated in directed acyclic graph representation. A treadth-first type checking pass, as seen in Algo- rithm [ completes in linear time, and GGA results may be applied to implement computeGGA () us- ing the following steps:
pass
Data: Abstract syntax tree for PPL
program
Result: GGA parameter estimates for all
random variables
frontier +- [rv Parents(rv) = 0];
tails (
while frontier * 0 do
next + frontier.
tails[next] < computeGGAAnext.op.
next.parent);
frontier <- frontier next.children();
end return tails
node has no parents, then is atomic distribution and its tail parameters are known (Table
Otherwise, the node is operation taking its potentially stochastic inputs (parents) to its output. Consult T for the output GGA tails.
3.2 Representative distributions
For each (v, a, p). we make a carefully defined choice of p on such that X ~ p. then X = (v, a, p). This way, any random variable f(X), where f is 1-Lipschitz, will exhibit the correct tail, and so approximations of this form may be used for VI or density estimation. Let X = (v,a,P) and 0 < 6 & denote a small parameter such that tails e are deemed to be "very heavy" (we chose 0.1). Inspired by ideas from implicit renewal theory [5] candidate distributions are as follows.
(p S 0) If p < 0, then px (z) c11 One such density the Student distribution, in this case, with |v| - 1 degrees of freedom if v < -1 (generate X ~ StudentT( [v] 1)).
(p > e) For moderately sized p > 0, the symmetrization of the generalized Gamma density 2b.
(p e) X (v. o. p) where is small, then will exhibit much heavier tails, and the generalized Gamma distribution in Case I will become challenging to sample from. In these cases, expect that the tail of should be well represented by power law. The generalized Gamma density (Equation
(2) satisfics E.X" = = (eppp(++++)) = ( for r > 0. Let a > be such that x 2. By Markov's 's inequality, the tail of X satisfies P(X > x) < 2 Therefore, we can represent tails of this form by the Student distribution with a + degrees of freedom (generate X ~ Studen(T(a)).
5

---

3.3 Bulk correction by Lipschitz mapping
While a representative distribution will exhibit the desired tails, the target distribution's bulk may be very different from generalized Gamma and result in poor distributional approximation. To address this, we propose splicing together the tails from a generalized Gamma with a flexible density approximation for the bulk. While many combinations Lipschitz operation in the GGA (Theorem 2 and post-compose neural spline flows [[5] (which identity functions outside of a bounded interval hence 1-Lipschitz) after properly initialized generalized Gamma distributions. Optimizing the parameters of the flow results in good bulk approximation while simultaneously preserving the tail correctness guarantees attained by the GGA.
possible, in this work we rely the
2.0
10
candidate
w/ flow target
Density
1.5
Density
10
1.0
10-3
0.5
10-5
0
5
10
10
10
10"
Figure 2: The candidate distribution chosen by the GGA calibrates tails to the target, but with incorrect bulk. A Lipschitz normalizing flow corrects the bulk (i) without changing the tail behaviour, as seen by the parallel tail asymptotics (black dashed lines) in (ii).
Example 1. Let A Igekk y E IR with F,, Aij iid
(- 1, 1). The distribution of T Ay = Eng Ayyy is convolution of normal-powert 21 and lacks a closed form expression. Using the
GGA (Table one can compute its tail parameters to be ( - 1, 3 2: 3) The candidate given by the GGA representative distribution (Section 3.2) is gamma distribution with correct tail behaviour, but a poor approximation otherwise. A learnable Lipschitz bijection is optimized to correct the bulk approximation (Figure ()). From the Lipschitz property, the slope of the tail asymptotics in log-log scale remains the same before and after applying the flow correction (Figure ii)): the tails guaranleed to remain calibrated.
Example 2. Consider Σt= X? where X, ~ (iudentT((). While we are not aware of a closed-form expression for the density, this example within the scope of our GGA. Empirical results illustrate that our method (Figure 3i)) accurately models both the bulk and the tail whereas Gaussian-biie based Lipschitz flows (Figure (i) inappropriately impose tails which decay too rapidly.
4 Experiments
We now demonstrate that GGA- density estimation yields improvements in tail estimation across several metrics. Our experiments consider normalizing flows initialized from (i) the parametric family defined Section 3.2 against (ii) a normal distribution (status quo). To further contrast the individual effect of using a GGA base distribution over standard normals against more expressive pushforward maps [5), we also report ablation results where normalizing flows are replaced by affine transforms, as originally proposed in K0 All experiments are repeated for 100 trials, trained to convergence using the Adam optimizer with manually tuned learning rate. Additional details are available in Appendix D All target distributions in this section are expressed as generative PPL programs: Cauchy using a reciprocal normal; Chi2 (chi-squared) using a sum squared normals; IG (Inverse Gamma) using reciprocal exponential; normal using a sum of normals; and StudentT using a normal and Cauchy ratio. Doing so tasks the static analyzer to infer the target's 's tails and makes the analysis non-trivial.
Our results the following tables share a consistent narrative where a GGA base distribution rarely hurts and can significantly help with heavy tailed targets. Standard evaluation metrics such as negative cross-entropy. ELBO, or importance- weighted autoencoder bounds (6] do not evaluate the quality of tail approximations. Instead, we consider diagnostics which do: namely, an estimated tail exponent

---

(i) GGA
(ii) Gaussian
10000
10000
7500
7500
Flow
Flow
5000
5000
2500
2500
2000
4000 6000 8000 10000
2000 4000 6000 8000 10000 Target
Target
Figure 3: Q-Q plots of density approximations of a heavy-tailed target E X? where X; ~ Studen(T(i)) initialized by our GGA candidate (i) and the Gaussian distribution (ii). While the expressive modeling capability of flows enables good approximation of the distribution bulk, Lipschitz transformations of Gaussians inevitably impose miscalibrated squared exponential tails which are sufficiently heavy as evidenced in (ii).
á, and the Pareto & diagnostic [60]. Except for when targets are truly light tailed oc in Chi2 and normal), GGA- -based approximations are the only ones to reproduce appropriate GPD tail index in density estimation and achieve a passing & below 0.2 in VI. Less surprising is the result that adding a flow improved approximation metrics, as we expect the additional representation flexibility to be beneficial.
Density Estimation. Given samples (z. from target density p. we minimize Monte-Carlo
estimate of the cross entropy (p. q) E,, [log (XX)
N EN log q(x:). The results are shown in Table 2 and Table 3 along with power-law tail index estimates â [I]]. Closeness between the target Pareto tail index a I1 and its estimate á in ((x) suggest calibrated tails. Overall, we see that normal (resp., Cauchy) based flows fails to capture heavy (resp., light) tails, while GGA-based flows yield good tail approximations (lower NLL, closer to target) across all cases.
Table 2:
Mean and standard errors (100 trials) of tail parameters å (smaller for heavier tails) for various density estimators and targets.
Target
Cauchy (a = 2) Flow
GGA Flow
Normal (a = o0) Flow
Cauchy
2
2.1 (0.03)
2.1 (0.07)
7.7 (2.5)
IG
1.9 (0.03)
1.9 (0.092)
7.3 (1.7)
StudentT
3
2.0 (0.06)
3.3 (0.45)
7.7 (2.3)
Chi2
oc
2.1 (0.07)
5.2 (1.6)
6.8 (2.4)
Normal
DC 2.9 (0.6)
8.2 (4.0)
8.4 (3.5)
Table 3: Mean and standard errors of log-likelihoode Ep log (X) for various density estimators and targets. While larger values imply a better overall approximation (row max bolded), log- likelioel dominated by bulk approximation so these results show that our method (GGA Flow) does not sacrifice bulk approximation quality.
Target
a
Cauchy (a 2) Flow
GGA Flow
Normal = oc) Flow
Cauchy
-.53 (0.05)
-3.22 (0.06)
-1.2 10 (6 x 10
IG
3355 (0.08)
-3.26 (0.05)
-266 10 (6 x 102)
StudentT
3
-2.12 (0.03)
-2.75 (0.04)
-2.92 (0.47)
Chi2
-2.30 (0.05)
-2.03 (0.04)
-2.24 (0.04)
Normal
-1.53 (0.03)
-1.41 (0.02)
-1.42 (0.02)
Variational Inference. For VI, the bulk is corrected through the ELBO optimization objective
x(X) N ΣN log p(a.) ~. q. Since the density p must also be evaluated, for simplic- ity, experiments in Table use closed-form marginalized densities for targets. The overall trends
Eg log p(X) R
also show that GGA yields consistent improvements; the & diagnostic (6)) indicates VI succeeds (h 0.2) when GGA with appropriately matched tails used and fails (i 1) when Gaussian tails are erroneously imposed.
7

---

Table 4: Pareto diagnostic ((60) to assess goodness of fit for VI (mean across 100 trials, standard
deviation in parenthesis) on targets of varying tail index (smaller a = heavier tails). A value > 0.2 is interpreted potentially problematic so only values not exceeding it are bolded.
Target
Normal Affine
Normal Flow
GGA Affine
GGA Flow
Cauchy
=
0.62 (0.26)
0.22 (0.059)
0.68 (0.038)
0.091 (0.04)
IG
a
8.6 (1.8)
8.2 (2.3)
2.0 (0.4)
2.9 (0.71)
StudentT
a
1.2 (0.16)
1.0 (0.43)
(0.082)
1.3 (0.097)
Chi2
a = oC
0.57 (0.081)
0.61 (0.067)
0.0093 (0.0067)
0.089 (0.044)
Normal
a DC
0.53 (0.17)
0.21 (0.067)
0.4 (0.086)
0.2 (0.089)
Bayesian linear regression. As a practical example of VI applied to posterior distributions, we consider the setting of one-dimensional Bayesian linear regression (BLR) with conjugate priors, defined by the likelihood ,| a ~ 3, o²) with a Gaussian prior 3|0² ~ (N on the coefficients, and inverse-Gamma Gamma prior with parameters ag and by on the residual variance o². The posterior distribution for 3 conditioned on o² and X, y is Gaussian. However, conditional on the pair (X, y). 0² is Gemmse-Gamma distributed with parameters ao + 2 and bu + 2 y p Σp)), where μ = I-1 X X 3 for the least-squares estimator, and Σ = x + Since o² is positive, it is typical for PPL implementations normalising flow starting from a Gaussian initialization will inappropriately approximate the inverse Gamma distributed p(o² XX y) with log-normal tails. On the other hand, Lipschitz flows starting from a GGA reference distribution will exhibit the correct tails. We assess this discrepancy in Figure 1 under an affine transformation on four subsampled datasets: super (superconductor critical temperature prediction dataset (23] with n = 256 and d = 154); who (life expectancy data from the World Health Organisation
apply exponential transformation. Hence, Lipschitz
the year 2013 41 with n 130, d 18): air (air quality data [14
with n = 6941, d = 11); and blog (blog feedback prediction dataset ( with n = 1024, d = 280). In Figure i), the GGA-based method seems to perfectly to the targets, while in Figure 4 ji), the standard Gaussian approach fails to capture the tail behaviour.
sppan
pla)x, y)
Sivge
logia
1000
100
-150
-150
-200
-200
10-1
10P
10A
10a
10-1
10a
10a
10a
0
Figure 4: Estimated densities for the posterior distribution of 0² in Bayesian linear regression under optimised exponential affine transformations from (i) GGA reference, and (ii) Gaussian reference.
Invariant distribution of SGD. For inputs X and labels Y from a dataset D, the least squares estimator for linear regression satisfies 3 = ming 1Ex.y=D(( - 3 To solve for this estimator, one can apply stochastic gradient descent (SGD) sampling over independent Xx, Ys D obtain the sequence of iterations
B&+1 = (I- -T- +Y,
for a step size 5 > O. For large 6, the iterates 3y typically exhibit heavy-tailed fluctuations 24 In this regard, this sequence of iterates has been used as simple model for more general stochastic optimization dynamics [22] 24). In particular, generalization performance has been tied to the heaviness the tails in the iterates [47. Here, we use our algebra to predict the tail behaviour in a simple one-dimensional setting where X ~ ~0, a²) and Y ~ (0, 1). From classical theory [ it known that X converges in distribution to a power law with tail exponent a > satisfying E _ 1. In Figure 5 we plot the density of the representative for 310 obtained using
8

---

algebra against kernel density estimate using 10th samples when E {0.4, 0.5, 0.6) and 6 = 2. In all cases, the density obtained from the algebra provides surprisingly close fit.
0 = 0.3
0= 0.4
d= 0.5
Deasiyy
Deasity
Density
10
10
10
10"
10"
10"
-10
10
Figure 5: Kernel density estimate of iterates of SGD (blue) vs. GGA predicted tail behaviour (orange)
Related Work
Heavy tails and probabilistic machine learning. For studying heavy tails, methods based subex- ponential distributions [8] and generalized Pareto distributions (GPD) (or equivalently, regularly varying distributions
44 have received significant attention historically. For example, [35 presents closure theorems for regularly varying distributions which are special cases of Proposition [ and Theorem Heavy tails often have a profound impact on probabilistic machine learning methods: in particular, the observation that density ratios pla))
q(x) tend to heavy tailed has resulted in new methods
for smoothing importance sampling [33 adaptively modifying divergences (55), and diagnosing VI
through the Pareto diagnostic [60]. These works are complementary to our paper, and our reported
results include & diagnostics for VI and á tail index estimates based on GPD.
Whereas those respective works lump the tail parameter in as another variational parameter and may be more generally applicable, the GGA may be applied before samples drawn and leads perfectly calibrated tails when applicable.
Our work considers heavy-tailed targets p(z) which is the same setting as 25 33
Probabilistic programming. PPLs can be broadly characterized by the inference algorithms they support, such as: Gibbs sampling over Bayes nets [13 48), stochastic control flow [19 58), deep
stochastic VI 4 12 Hamiltonian Monte-Carlo 8 59 Our implementation target beanmachine [50 is a declarative PPL selected due to availability of a PPL compiler and support for static analysis plugins. Similar to
4 [6]. it uses PyTorch 40 for GPU tensors and automatic differentiation.
Synthesizing an approximating distribution during PPL. compilation (Section 3 is also performed
in the Stan language by 30 and normalizing flow extensions in [57] We compare directly against.
these related density approximators in Section 4
Static analysis. There is long history of formal methods and probabilistic programming in the literature 226 22), with much of the research [0O concerned with defining formal semantics and [54] (see 3 for a recent review). Static analysis uses the abstract syntax tree (AST) representation of a program in order to compute invariants (e.g., the return type of a function, the number of classes implementing trait) without executing the underlying program. has traditionally been applied in the context of formalizing semantics [29]. and has been used to verify probabilistic programs by ensuring termination, bounding random values values 44 As dynamic analysis in a PPL. is less reliable due to non-determinism, static analysis techniques for PPLs become essential. As recent examples,
establishing invariants
4 to verify distribution supports and avoid - IIf log probabilities. More relevant to our work are applications of static analysis to improve inference. (38) and [22 both employ static analysis inform choice of inference method. However, both works do not account for heavy tails whereas the primary goal of BGA-based analysis is to ensure tails are properly modelled.
(32) proposes static analyzer for the Pyro PPL.
Conclusion
In this work, we have proposed a novel systematic approach for conducting tail inferential static PPL. analysis. We have done this by defining heavy-tailed algebra, and by implementing three- parameter generalized Gamma algebra into PPI. compiler. Initial results are promising, showing that improved inference with simpler approximation families is possible when combined with tail

---

metadata. While already useful, the generalized Gamma algebra and its implementation currently has some notable limitations:
The most significant omission to the algebra is classification of log-normal tails. Addition may be treated using [20], but multiplication with log-normal tails remains elusive.
Since the algebra assumes independence. handling of dependencies between defined random variables must be conducted externally. This can be addressed using symbolic package to decompose complex expressions into operations on independent random variables.
Scale coefficients for conditional distributions may often be inexact, as exact marginalization in general is NP-hard [28]. Treatment of disintegration using symbolic manipulations is significant open problem, with some basic developments 45
Compile-time static analysis is only applicable to fixed model structures. Control flow, open-
universe models [36]. and PPI.s to support them 4 are an important future research direction.
References
[1] Søren Asmussen and Hansjorg Albrecher. Ruin probabilities, volume 14. World scientific, 2010.
[2] Søren Asmussen, Enkelejd Hashorva, Patrick J Laub, and Thomas Taimre. Tail asymptotics of light-tailed Weibull-like sums. Probability and Mathematical Statistics, 37(2):235-25) 256, 2017.
[3] Ryan Bernstein. Static analysis for probabilistic programs. arXiv preprint 1909.0599 2019.
[4] Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman. Pyro: Deep universal probabilistic programming. The Journal of Machine Learning Research, 20(1):973- 978, 2019.
[5] Dariusz Buraczewski, Ewa Damek, Thomas Mikosch, al. Stochastic models with power-law tails. Springer Ser. Oper. Res. Financ. Eng. Springer, Cham, 10:978-3, 2016.
[6] Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. Importance weighted autoencoders. arXiv preprint arXiv: 1509.00519, 2015.
[7] Krisztian Buza. Feedback prediction for blogs. Data analysis, machine learning and knowledge discovery, pages 152-152. Springer, 2013.
[8] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 76(1), 2017.
[9] Kenta Cho and Bart Jacobs. Disintegration and Bayesian inversion via string diagrams. Mathe- matical Structures in Computer Science, 29(7):938-99 2019.
[10] Guillaume Claret, Sriram K Rajamani, Aditya V Nori, Andrew D Gordon, and Johannes Borgström. Bayesian inference using data flow analysis. In Proceedings of the 2013 9th Joint Meeting Foundations of Software Engineering, pages 92-102, 2013.
[11] Aaron Clauset, Cosma Rohilla Shalizi, and Mark EJ Newman. Power-law distributions in empirical data. SIAM review, (14666611 2009.
[12] Marco
Feras Saad, Alexander K Lew, and Vikash K Mansinghka. Gen: a general-purposs probabilistic programming system with programmable inference. In Proceed- ings of the 40th acm sigplan conference programming language design and implementation. pages 221 2-36.3, 2019.
[13] Perry de Valpine, Daniel Turek, Christopher J Paciorek, Clifford Anderson- Bergoman, Dun- can Temple Lang, and Rastislav Bodik. Programming with models: writing statistical algorithms for general model structures with nimble. Journal of Computational and Graphical Statistics, 26(2):403-413 2017.
10

---

[14] Saverio De Vito, Ettore Massera, Marco Piga, Luca Martinotto, and Girolamo Di Francia. On field calibration of electronic nose for benzene estimation an urban pollution monitoring scenario. Sensors and Actuators B: Chemical, 129():750--577, 2008.
[15] Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances in neural information processing systems, 32, 2019.
[16] SM Eslami, Nicolas Heess, Theophane Weber, Yuval Tassa, David Szepesvari, Geoffrey Hinton, et al. Attend, infer, repeat: Fast scene understanding with generative models. Advances in Neural Information Processing Systems, 29, 2016.
[17] Mark Girolami and Ben Calderhead. Riemann manifold Langovin and Hamiltonian monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 23(2):123-214, 2011.
[18] Charles M Goldie and Claudia Klüppelleegg. Subexponential distributions. A practical guide to heavy tails: statistical techniques and applications, pages 435-459, 1998.
[19] Noah Goodman, Vikash Mansinghka, Daniel M Roy, Keith Bonawitz, and Joshua B Tenenbaum. Church: language for generative models. arXiv preprint urXiv:1206.322 1206.3255, 2012.
[20] Archil Gulisashvili and Peter Tankov. Tail behavior of sums and differences of log-normal random variables. Bernoulli, 22(1):444-493, 2016.
[21] Rameshwar D Gupta and Ramesh C Gupta. Analyzing skewed data by power normal model. Test, 17(1):197220, 2008.
[22] Mert Gurbuzbalaban, Umut Simsekli, and Lingjiong Zhu. The heavy-tail phenomenon in SGD,
In International Conference on Machine Learning, pages 3964-3975. PMLR, 2021.
[23] Kam Hamidich. data-driven statistical model for predicting the critical temperature of
superconductor. Computational Materials Science, 154:346-355 354, 2018.
[24] Liam Hodgkinson and Michael W. Mahoney. Multiplicative noise and heavy tails in stochastic optimization. In International Conference on Machine Learning, pages 4262-4274. PMLR, 2021.
[25] Priyank Jaini, Ivan Kobyzev, Yaoliang Yu, and Marcus Brubaker. Tails of Lipschitz triangular flows. In International Conference on Machine Learning, pages 4673-4681. PMLR, 2020.
[26] Claire Jones and Gordon D Plotkin. A probabilistic powerdomain of evaluations. In Proceedings. Fourth Anmual Symposium on Logic in Computer Science, pages 186-187. IEEE Computer Society, 1989.
[27] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 2014.
[28] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
[29] Dexter Kozen. Semantics of probabilistic programs. In 20th Annual Symposium on Foundations of Computer Science (sfes 1979), pages 101-114. IEEE, 1979.
[30] Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M Blei. Automatic differentiation variational inference. Journal of machine learning research, 2017.
[31] Michel Ledoux. The concentration of measure phenomenon. American Mathematical Soc., 2001.
[32] Wonyeol Lee, Hangyeol Yu, Xavier Rival, and Hongseok Yang. Towards verified stochastic variational inference for probabilistic programs. Proceedings of the ACM on Programming Languages, 4(POPL):1-33 2019.
[33] Feynman Liang, Liam Hodgkinson, and Michael W. Mahoney. Fat-tailed variational inference with anisotropic tail adaptive flows, 2022.

---

[34] Arakaparampil M Mathai, Ram Kishore Saxena, and Hans Haubold. The H H-fumction: theory and applications. Springer Science & Business Media, 2009.
[35] Mikosch. Regular variation subexponentiality and their applications in probability theory,
1999.
[36] Brian Milch and Stuart Russell. Extending Bayesian networks to the open-universe case. Heuristics, Probability and Causality: A Tribute to Judea Pearl. College Publications, 2010.
[37] Nair, A. Wierman, and B. Zwart. The Fundamentals of Heavy Tails: Properties, Emergence, and Estimation. Cambridge University Press, 2022.
[38] Aditya Nori, Chung-Kil Hur, Sriram Rajamani, and Selva Samuel. R2: An efficient MCMC sampler for probabilistic programs. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 28, 2014.
[39] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1 2021.
[40] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
[41] Kumar Rajarshi. Life expectancy (who), 2018.
[42] Gareth O Roberts and Richard L Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, pages 341-363, 1996.
[43] Barry K Rosen, Mark N Wegman, and F Kenneth Zadeck. Global value numbers and redundant computations. In Proceedings of the 15th ACM SIGPLAN-SIGACT symposium on Principles of programming languages, pages 12-27, 1988.
[44] Sriram Sankaranarayanan, Aleksandar Chakarov, and Sumit Gulwani. Static analysis for probabilistic programs: inferring whole program properties from finitely many paths. In Proceedings of the 34th ACM SIGPLAN conference on Programming language design and implementation, pages 447-458. 2013.
[45] Chung-chier Shan and Norman Ramsey, Exact Bayesian inference by symbolic disintegration. In Proceedings of the 44th ACM SIGPLAN Symposium on Principles of Programming Lan- guages, pages 130-144, 2017
[46] N. Siddharth, Brooks Paige, Jan-W Willem van de Meent, Alban Desmaison, Noah D. Goodman, Pushmeet Kohli, Frank Wood, and Philip Torr. Learning disentangled representations with semi-supervised deep generative models. In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems 30, pages 5927-5937. Curran Associates, Inc., 2017.
[47] Umut Simsekli, Levent Sagun, and Mert Gurbuzbalaban. A tail-index analysis of stochastic gradient noise in deep neural networks. In International Conference on Machine Learning, pages 5827-5837. PMLR, 2019.
[48] David Spiegelhalter, Andrew Thomas, Nicky Best, and Wally Gilks. BUGS 0.5: Bayesian inference using Gibbs sampling manual (version ii). MRC Biostatisties Unit, Institute of Public Health, Cambridge, UK, pages 1-59, 1996.
[49] Nader Tajvidi. Confidence intervals and accuracy estimation for heavy-tailed generalized Pareto distributions. Extremes, 6(2): 11 11-123, 2003.
[50] Nazanin Tehrani, Nimar S Arora, Yucen Lily Li, Kinjal Divesh Shah, David Noursi, Michael Tingley, Narjes Torabi, Eric Lippert, Erik Meijer, et al. Bean machine: A declarative probabilistic programming language for efficient programmable inference. International Conference on Probabilistic Graphical Models. PMLR, 2020.
12

---

[51] David Tolpin, Man-Willem van de Meent, Hongseok Yang, and Frank Wood. Design and implementation of probabilistic programming language anglican. In Proceedings of the 28th Symposium on the Implementation and Application of Functional progrumming Languages, pages 1-12, 2016.
[52] Dustin Tran, Matthew Hoffman, Dave Moore, Christopher Suter, Srinivas Vasudevan, and Alexey Radul. Simple, distributed, and accelerated probabilisticc programming. Advances in Neural Information Processing Systems, 31, 2018.
[53] Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry, Pareto smoothed importance sampling. arXiv preprint arXiv::507, 026400266, 2015.
[54] Di Wang, Jan Hoffmann, and Thomas Reps. Pmaf: an algebraic framework for static analysis of probabilistic programs. ACM SIGPLAN Notices, 53(4):513 528, 2018.
[55] Dilin Wang, Hao Liu, and Qiang Liu. Variational inference with tail-adaptive f-divergence. Advances Neural Information Processing Systems, 31, 2018.
[56] George Neville Watson. treatise on the theory of Bessel functions. Cambridge university press, 1995.
[57] Stefan Webb, Jonathan Chen, Matrin Jankowiak, and Noah Goodman. Improving automated variational inference with normalizing flows. In ICML Workshop on Automated Machine Learning, 2019.
[58] David Wingate, Andreas Stuhlmüller, and Noah Goodman. Lightweight implementations of probabilistic programming languages via transformational compilation. In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statisties, pages 770-778.
JMLR Workshop and Conference Proceedings, 2011.
[59] Kai Xu, Hong Ge, Will Tebbutt, Mohamed Tarek, Martin Trapp, and Zoubin Ghahramani. Advancedhme. jl: A robust, modular and efficient implementation advanced hme algorithms. In Symposium on Advances in Approximate Bayesian Inference, pages 1-10. PMLR, 2020.
[60] Yuling Yao, Aki Vehtari, Daniel Simpson, and Andrew Gelman. Yes, but did it work?: Evaluat- ing variational inference. In International Conference on Machine Learning, pages 5581-5590. PMLR, 2018
13


Reviews:
[{"review_id": "aiTc3uuX6d", "reviewer": "Reviewer_9rGq", "summary": "The paper addresses the problem of density estimation of probabilistic models (expressed as probabilistic programs), with a focus on their tails. This is important for several Bayesian inference methods: importance sampling can exhibit infinite variance if the proposal has a lighter tail than the target, many black box variational inference methods cannot change their tail behavior and MCMC algorithms may lose ergodicity if the tail of the target is outside a particular class. This paper proposes a static analysis pass on probabilistic programs consisting of sampling from primitive distributions and arithmetic operations on the sampled variables (addition, multiplication, division, logarithm, exponential, and Lipschitz functions). The analysis is based on generalized Gamma distributions, which capture many tail behaviors of primitive distributions (e.g. Gaussian, Gamma, Weibull, Student-t, Pareto, but not the log-normal) and are closed under the above operations on random variables (or can be bounded if not represented exactly). As a consequence, a compiler pass can compute the parameters of a generalized Gamma distribution with the same (or at least providing a bound on the) tail behavior as the posterior distribution of the probabilistic program. Since the bulk of this distribution may be very different from the posterior, the paper combines it with neural spline flows: effectively, this computes the pushforward distribution under a Lipschitz function, resulting in a better bulk approximation while preserving the tails.\nThe experimental evaluation demonstrates that this approach to density estimation usually improves the tail estimation compared to status-quo initializations with a normal or Cauchy distribution, both for density estimation using normalizing flows and variational inference, where the target distributions are simple Cauchy, Inverse Gamma, Student-t, χ^2 and normal distributions. Similarly, the density estimation of the variance parameter of a Bayesian linear regression model is improved compared to the standard Gaussian approach. Finally, they observe that their method (without normalizing flows) provides a very close fit for the density estimation of the iterates of stochastic gradient descent applied to a least-squares linear regression.", "questions": "- Line 92 claims that closure under some operations is known. Could you provide references for this?\n- What exactly are the guarantees and assumptions for operations under which GGA is not closed? What happens if such operations are applied repeatedly? (cf. Weaknesses)\n- The operations in Table 1 only apply to equivalence classes $(\\nu, \\sigma, \\rho)$. What happens if the result of an operation is $\\mathcal L$ or $\\mathcal R_1$? Will those be treated as $(?, ?, \\infty)$ or $(-1, ?, 0)$ in subsequent operations? If so, what values are assumed for the question marks? Why don't the infinities lead to problems?\n- Section 3.3: could you elaborate on why the normalizing flows don't affect the tails (cf. Weaknesses)?\n- In Table 4, GGA Flow does not seem to work for IG and StudentT. Do you know what the reason for this is? It seems unexpected given all the claimed benefits.\n\nTypos:\n- Equation (1): should be $= cx^\\nu ...$ or $\\sim x^\\nu ...$\n- line 111: missing indices $f_{ij}$, $R_{ij}$.\n- line 131: shouldn't it say $\\rho \\le 0$ instead of $\\rho \\le -1$?\n- line 201: it is unclear where the \"conditional on\" clause ends, please rephrase\n", "limitations": "Limitations are discussed in the Conclusion. I believe a few items should be added to this list:\n- no control structure in the programming language (it is an expression language, cf. Weaknesses)\n- tails are not exact (only bounds) if certain operations are involved", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "The paper's main idea of using a static analysis to compute the tails of the posterior distribution of a probabilistic program is very interesting because tails are difficult to estimate. Such an ahead-of-time analysis with correctness guarantees is therefore very useful. In order to arrive at this result, the paper proves several new theoretical results on generalized Gamma distributions, which may be of independent interest. The experimental results show that in combination with normalizing flows, this often performs better than existing methods that assume a Gaussian base distribution. The presentation of the paper is generally clear and readable.\n\nOverall, this paper has an innovative idea with promising results, which are of interest to the probabilistic programming community.", "weaknesses": "While the paper's idea is very interesting, I have a few concerns about the soundness and presentation. I will update my rating if the author's answers address my concerns.\n\nGuarantees and assumptions: What exactly are the guarantees (mentioned in line 149) and assumptions of the static analysis? It would be good if this could be stated in one place, ideally as a self-contained theorem. Table 1 mentions that \"additional assumptions are required\" for several operations, but I was not able to find them in the paper. The assumption that the operands are independent is also very important and only really mentioned at the very end (line 270). It is also unclear what exactly the guarantees are for programs involving reciprocals, exponentials, and logarithms because the generalized Gamma distributions are not closed under these operations. In this case, will the computed tail be an upper or a lower bound? What happens if such operations are applied repeatedly? While there are theorems and proofs for each individual operation, I was not able to find a theorem clarifying the assumptions and guarantees about the composition of such operations, especially if the resulting tail cannot be represented exactly in the GGA.\n\nBulk correction: Section 3.3 (the bulk correction via normalizing flows) is informal and light on details. I assume the neural spine flows preserve the tails because they are 1-Lipschitz maps? This is not obvious to me, and I think it would be good if the paper could elaborate on why the normalizing flows don't affect the tails. Again, putting the statement in the form of a theorem with all necessary assumptions would be helpful.\n\nProbabilistic programming language: The grammar of the probabilistic programming language is only introduced by example, there is no comprehensive description, not even in the appendix. Furthermore, the name \"programming language\" sounds like there is support for control structure (branches, loops), which doesn't seem to be the case in this approach. Rather, the \"programming language\" seems to be a simple expression language with the grammar `E ::= sample(D) | op(E, ..., E)` where `E` is an expression, `op` is a supported operation, and `D` is a supported primitive distribution. It would be good if this could be clarified in the paper.\n\nIndependence assumption: a main assumption is that the operands to all operations are independent random variables. This seems to be difficult to guarantee in practice because posteriors of parameters are often correlated in probabilistic models. Does this mean that (for now), all operands of all operations need to be hand-checked for independence? Does this severely restrict the class of probabilistic models that can be handled by this approach? If so, this should be discussed. It is only briefly mentioned in line 270.\n\nThe following points are weaknesses, but less important than the previous points:\n\nExperimental evaluation: the experiments demonstrate benefits for very simple benchmarks (one-dimensional primitive distributions, one-dimensional Bayesian linear regression). If more realistic experiments (e.g. higher-dimensional examples) could be performed, this would strengthen the paper. However, I don't consider this a condition for publication because the results are already interesting.\n\nAbstract interpretation: the static analysis pass, as I understand it, performs abstract interpretation (a common technique in static analysis) where the abstract domain is given by the generalized Gamma algebra. It would be good to mention this keyword, as it helps put this method into context.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\nThank you for taking the time to provide a detailed review and giving us the opportunity to address your concerns. We will discuss each point not addressed in our overall author response sequentially.\n\n## Weaknesses\n\n**Guarantees and assumptions**:\nThank you for the excellent suggestion! To improve readability of Appendix A in the revision, we will book-end the results with the following main theorem. \n\n**Main Theorem**: Let $X_1,\\dots,X_n$ be independent random variables with generalized Gamma tails. \nThen\n- For any $i \\neq j$, $c \\neq 0$ and $\\beta > 0$, $X_i + X_j$, $X_i \\times X_j$, $c X_i$, and $X_i^\\beta$ also exhibit generalized Gamma tails as detailed in Table 1 (addition, multiplication, scalar multiplication, powers).\n- If $X_i$ and $X_j$ have densities $p_i$ and $p_j$, then the random variable with density proportional to $p_i p_j$ has generalized Gamma tails as in Table 1 (product of densities).\n- If $X_i$ is generalized Gamma distributed, then $\\frac{1}{X_i}$ has _exactly_ the generalized Gamma distribution described in Table 1 (reciprocals). \n- If $X_i$ has density $p$ continuous at zero with $p(0) > 0$, then $\\frac{1}{X_i}$ also has generalized Gamma tails described in Table 1 (reciprocals).\n- For any Lipschitz function $f$, the tail of $f(X_1,\\dots,X_n)$ is no heavier than the generalized Gamma tail detailed in Table 1 (Lipschitz functions). If $f$ is asymptotically linear then $f(X_1,\\dots,X_n)$ has _exactly_ the described tail. \n- $\\log X_i$ and $e^{X_i}$ have tails no heavier than the generalized Gamma tails in Table 1 (exponentials, logarithms). If $X_i$ is regularly varying, then $\\log X_i$ has  _exactly_ the described tail. If $X_i$ is exponentially distributed, then $\\exp X_i$ has _exactly_ the described tail. \n\nHopefully this clears up any confusion around the theoretical results.\n\n**Repeated application**: \nProvided independence holds, composition of operations in the GGA remain consistent unless one applies Lipschitz functions, logarithms, or exponentials. If one of these operations is applied, the tail becomes an upper bound, which remains consistent under addition, multiplication, and powers, but not reciprocals. Given that we are working with a fixed class of tails, such behavior is inevitable, and it is possible to perform a sequence of operations for which the tail no longer becomes accurate. Nevertheless, we have endeavored to create a system to balance theoretical consistency with practical value, and find that in most cases, the GGA algebra produces a tail that is either exact, or a very good estimate (the SGD example was particularly convincing to us). \n\n**Bulk correction**: Yes, neural spline flows are 1-Lipschitz linear functions in the tail. We implied this by writing \"identity functions outside of a bounded interval\" (line 147) but will revise to be more explicit for clarity.\n\n**Probabilistic programming languages**: We would like to point out the distinction between the host PPL and the symbolic expressions where the GGA is applicable. The GGA itself is not a PPL, but rather an ahead-of-time analysis method compatible with a PPL. Our host PPL (beanmachine) is based in Python and capable of branching and looping. At the end of a single execution trace, the expressions constituting a target random variable is analyzed by the GGA to determine its tail parameters.\n\nPerhaps more interesting, and indeed outside of the scope of the GGA we have developed, is the ability to directly analyze control flow structures. However, we note that analysis of a mixture random variable generated by branching control flows requires marginalization where even approximate inference is known to be NP-hard (line 274 and Koller and Friedman, 2009). Additional assumptions such as junction-tree structure or restricted distribution families may permit further analysis and we leave this for future work.\n\n**Experimental evaluation**: We have restricted our experiments to those where we know the solution exactly, as it is difficult to provide a good frame of reference otherwise.\n\n**Abstract interpretation**: Thank you for this keyword, indeed, this describes our approach very well. We will include it in the revised version.\n\n## Questions\n\nResponses to questions not already addressed by the \"Weaknesses\" section:\n\n - Line 92 references --- These results are cited in Appendix A, but we are happy to include these citations here as well. \n - Table 1 $\\mathcal{L}_1$ or $\\mathcal{R}_1$ results --- Again, we are taking some liberties to encode distributions that do not fit into the generalized Gamma tail framework. $\\mathcal{L}$ is treated as $(0,1,\\infty)$ and $\\mathcal{R}_1$ is treated as $(-1,1,0)$. Any operation which results in a super-light tail ($\\rho = \\infty$) becomes $\\mathcal{L}$ and any operation which results in a super-heavy tail ($\\rho \\leq 0$ and $\\nu \\geq -1$) becomes $\\mathcal{R}_1$. Each of the operations are consistent under these definitions, where $1/\\infty$ is treated as $0$. One deficiency here that we recently became aware of is that multiplication with a constant is not always equal to multiplication by $\\mathcal{L}$ --- this is the one operation where multiplication by $\\mathcal{L}$ can give different answers depending on how $\\nu,\\sigma$ is defined. This is unavoidable, but our procedure still offers an excellent heuristic even in this setting. In any case, we will include this discussion in the revised version. \n- Table 4 IG and StudentT --- The Pareto $\\hat{k}$ diagnostic (Yao, Vehtari, Simpson, Gelman, 2018) is based on a heuristic argument and as a test for divergence is known to yield false positives, which may be the  case in the two examples mentioned."}, {"review_id": "wRiLN7caSV", "reviewer": "Reviewer_E196", "summary": "During inference, we are often interested in the behavior of the tails of the distributions we are analyzing. Heavy or light tails may necessitate switching algorithms so that inference remains stable for example. This paper describes a calculus by which a probabilistic programming language may calculate the tails of distribution objects derived from operations performed on elemental distributions. They apply their method to perform some simple probabilistic programming applications.", "questions": "In line 134 you describe the representative distribution of a distribution with a tail with a small $\\rho$. You motivate your construction by matching moments, but could you describe more rigorously why this is a reasonable \n\nI also wonder why the authors don't consider incorporating discrete distributions such as Poisson and Geometric distributions, which can be treated as continuous in their setting, say by adding a uniform random variable, so a geometric would have a $(0, 0, \\log(p))$ tail for example. There are often cases where discrete data is treated as continuous, such as for RNAseq counts.", "limitations": "The authors describe limitations in section 6.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The authors apply their method to a huge number of distributions and their transformations. They summarize their rules helpfully in charts. It seems as though their method may also simply be added to other probabilistic programming languages. I would certainly appreciate being able to call the tail parameters of a distribution object in Pyro for example.", "weaknesses": "The authors in section 3 suggest a method for fitting a density that consists of fitting its tail using their calculus and then fitting the bulk using a flow model. This is a reasonable choice for their simple experiments and there are a number of other choices they could have made to fit the bulk. The exposition of the paper sets it up in figure 1 and section 3 as essential to their probabilistic programming method, but I don't believe that to be that case. It would clarify and strengthen the paper if authors clarified the relation of their flow method for estimating densities to the rest of their contributions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the positive assessment of our work! Yes, you are correct that the choice of spliced flow is non-essential. Our intention was to propose one such construction which was sufficiently flexible to capture the bulk while also respecting tail asymptotics computed by the GGA and we achieved this by proving Theorem 2 and operationalizing it in Section 3. Alternative constructions which preserve tail asymptotics may also perform well and we will clarify our contributions section to reflect this in the revision.\n\n> In line 134 you describe the representative distribution of a distribution with a tail with a small . You motivate your construction by matching moments, but could you describe more rigorously why this is a reasonable \n\nIt is known that very heavy subexponential distributions (small $\\rho$)---which often come from repeated multiplications---closely resemble a power law in practice, as shown in line 522. If we do not approximate this case by a power law, the representative generalized Gamma distribution often provides a very poor approximation to the bulk, inhibiting the final stage of bulk correction. The moment-matching approach is loosely derived from ideas in implicit renewal theory (Buraczewski, Dariusz, Damek, Mikosch 2016, which we are happy to cite in the revised version) and we found it to perform the best in practice. \n\n__Incorporating discretes__: Thank you for the suggestion! Yes, as you have already noted with the geometric distribution (although the class would be $(0,\\log p, 1)$), there is nothing to prohibit the use of the GGA for operating on discrete distributions by treating them as continuous through a suitable convolution as suggested."}, {"review_id": "9bkhG3sFwc", "reviewer": "Reviewer_mtnU", "summary": "The paper develops an algebra which acts on a three-parameter family of tail asymptotics based on the generalized Gamma distribution. The algebraic operations are closed under addition, multiplication, powers, and a full list is given in Table 1. With this algebra, tail calculation can be done automatically in probabilistic programming instead of analytically. The paper also proposed a method of splicing bulk and calibrated tails with experimental verification.", "questions": "1. In the case of multi-dimensional random variables and dependent tails, would you consider using an extreme-value copula?\n2. In the case of dependent tails and bulk to tail splicing/ calibration/ extrapolation, would you consider using an Archimax copula?", "limitations": "The limitations were clearly defined and addressed in the paper.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper develops a heavy-tailed algebra for probabilistic programming and also proposed a method of splicing bulk and calibrated tails with experimental verification.\n2. Example 4 in the appendix had a calculation for the product of two random variables. ", "weaknesses": "I did not notice any glaring weaknesses in the paper but I am also not a subject matter expert on probabilistic programming.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for reading and reviewing our paper. Since our work only focuses on univariate tails, we agree that copulas provide a promising direction which can be combined with other work applying normalizing flows to multi-variate heavy tails (ATAF paper) in order to improve multivariate heavy-tailed approximation. The choice of copula is indeed important here, so the extension will require further work. "}, {"review_id": "bHTspGFe2p", "reviewer": "Reviewer_bd7W", "summary": "The paper proposes a static analysis technique for probabilistic programming languages, which annotates random variables with metadata characterizing their tail behavior. In particular, generalized Gamma distributions are used for this purpose. It is shown that they are closed under a number of operations, including addition, multiplication, and, under some conditions, reciprocals. This gives rise to an algebra which is used to statically infer the tail behavior of random variables throughout the probabilistic program, including posterior distributions. At runtime, the random variables are then estimated by neural splice flows initialized to the inferred Gamma distributions, which guarantees the correct tail behavior while allowing flexibility for the bulk of the probability mass. It is shown that this approach yields the correct behavior in many instances in which conventional methods fail.", "questions": "- How does the system behave when the assumptions discussed above are not fully met, e.g., when there is an operation between dependent variables? Is it possible to make use of partial results for a given program?\n- The discussion on posterior distributions mostly focuses on two random variable connected by a chain of operations, but notes that other cases may also be covered. Is it really the case that any posterior distribution may be estimated this straightforwardly? Can the system as implemented cover complicated multivariate computation graphs?", "limitations": "The limitations have been well addressed, though their practical implications could be stated more directly (as discussed above).", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": " - The paper addresses an important and interesting problem in probabilistic programming.\n - The static analysis approach to analyzing tail behavior is novel to my knowledge.\n - The chosen approach is competently executed: The provided algebra covers a wide variety of distribution types and operations, many of which utilizing newly derived theoretical results.\n - Combining the inferred Gamma distribution parameters with neural splice flows is a clever way to maintain the correct tail behavior while allowing flexible estimation of the distribution's bulk.\n - Encouraging experimental results are provided for a number of tasks, including density estimation, variational inference, and Bayesian linear regression.\n - The paper is well written and easy to follow. Detailed derivations for the new theoretical results are provided in the appendix.\n", "weaknesses": "- The models used for the experimental results are somewhat small. While simple monovariate distributions are fine to get a sense for the behavior of the method, it would be good to include some larger programs to demonstrate the robustness of the method.\n- The paper mentions a number of limitations of the approach: It does not cover log-normal tails, operations between dependent variables, or dynamic model structures. Given these limitations, it would be especially useful to highlight some set of more complex models used in practice to which the method is nevertheless applicable.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the positive assessment of our work! \n\n> How does the system behave when the assumptions discussed above are not fully met, e.g., when there is an operation between dependent variables? Is it possible to make use of partial results for a given program?\n\n> The discussion on posterior distributions mostly focuses on two random variable connected by a chain of operations, but notes that other cases may also be covered. Is it really the case that any posterior distribution may be estimated this straightforwardly? Can the system as implemented cover complicated multivariate computation graphs?\n\nThis relates to questions raised by Reviewer 9rGq's --- in a nutshell, handling dependencies accurately will require symbolic treatment. \nYou are right that not all posteriors are covered, and it is actually NP-hard (line 274) to cover all possible posteriors. While the GGA is capable of analyzing some complicated multivaraite computation graphs, it is currently limited to computation graphs satisfying (or as the discussion around __independence assumptions__ shows,\ncan be equivalently rewritten to satisfy) the independence assumption for the input operands of every GGA operation.\nGenerally speaking, to exactly determine tail behavior at scale is likely to require further advancements on top of what we propose. \nHowever, we do agree that there is still partial utility when results are violated and will include this discussion in the revision. \nIn particular, since the GGA is quite simple, there is nothing preventing its usage at scale to provide an estimate of the tail behavior. To our knowledge, there is no alternative in the literature, and the GGA may reveal information about the target density (e.g. to diagnose issues) that would otherwise be opaque. \n\n> ...the method's evaluation is limited to unimodal distributions. The density estimation experiments have been done on standard distributions. I am unsure about whether this is a limitation of the evaluation or a limitation of the proposed method.\n\nWe believe evaluation against a multi-modal mixture target is reasonable provided the mixture is given explicitly (e.g. as a list of component random variables and their mixture probabilities) so that a PPL compiler can analyze it statically. \nDue to space constraints, we did not include an operation for a mixture of densities in the GGA although such a result is possible within the GGA. Indeed, denoting a mixture of densities by $\\cup$, this operation would be added to Table 1 in the form\n$$\n(\\nu_1,\\sigma_1,\\rho_1) \\cup (\\nu_2,\\sigma_2,\\rho_2) \\equiv \\max\\{(\\nu_1,\\sigma_1,\\rho_1),(\\nu_2,\\sigma_2,\\rho_2)\\}.\n$$\nAt this point, the framework is capable of naturally handling multi-modality through its\nnormalizing flow bulk adjustment.\nFor example, (Wehenkel \\& Louppe, 2019) or (Durkan et al., 2019) both provide multi-modal density\napproximations and could be used as the normalizing flow.\nSince multimodality is a bulk property (e.g. a Gaussian mixture with finitely many components continues to have Gaussian tails), it does not affect tail asymptotics\nand the GGA remains applicable.\nWith this in mind, multimodality is not really an assessment of the GGA, and so we focussed our attention on unimodal distributions with clearer comparisons. We are happy to include a multimodal example in the revised version, however. "}]

Author rebuttal:
We thank the reviewers for their time and valuable feedback. We appreciate that reviewers have recognized the novelty of our approach and its applications in ahead-of-time static analysis of probabilistic programming languages (PPL). Suggested minor changes and fixes have been incorporated into the working document. Below we respond to points shared by more than one reviewer and provide more specific responses to individual reviewers as comments.

**Independence assumption**: Independence of operands is evidently key to the algebra functioning correctly, as for example, the tail of $X^2$ differs from the tail of $X Y$ even if $X$ and $Y$ share the same tail. We will include a few words in Section 2 to be clear about this. To deal with arbitrary dependencies that arise due to composition requires a combination of symbolic methods with the GGA. For example, if $Y = X_1 + X_2$ and $Z = X_2 + X_3$, then even if $X_i$ are independent the GGA cannot be applied as-is to $Y + Z$. However, the equivalent symbolic representation $Y + Z = X_1 + 2 X_2 + X_3$ can be dealt with using the GGA. This is discussed in the conclusion (line 270), but we can include further discussion on this matter in the revision.

Meta-review: The authors give a systematic approach to analyzing the tail of random variables in the static analysis approach of a probabilistic programming language compiler, with accompanying empirical results. This is an important question as the tail behavior influences the choice/performance of MCMC and variational inference procedures. In response to reviewer questions, authors explicit stating guarantees for how tails of random variables transform under a catalogue of operations, which clarifies the theoretical grounding for their work. Limitations include the restriction to univariate distributions and independence assumptions, but the paper remains an excellent initial foray into the problem.

Reviewers gave uniformly high rating after discussion. I recommend acceptance.
Decision: Accept (poster)