Exact Bayesian Inference on Discrete Models via Probability Generating Functions: A Probabilistic Programming Approach
Authors: Fabian Zaiser
Andrzej S Murawski
Luke Ong
Venue: NeurIPS 2023 oral (2023)
Keywords: Bayesian statistics
probabliistic programming
exact inference
discrete models
probability generating functions

Abstract: We present an exact Bayesian inference method for discrete statistical models, which can find exact solutions to a large class of discrete inference problems, even with infinite support and continuous priors.
To express such models, we introduce a probabilistic programming language that supports discrete and continuous sampling, discrete observations, affine functions, (stochastic) branching, and conditioning on discrete events.
Our key tool is *probability generating functions*:
they provide a compact closed-form representation of distributions that are definable by programs, thus enabling the exact computation of posterior probabilities, expectation, variance, and higher moments.
Our inference method is provably correct and fully automated in a tool called *Genfer*, which uses automatic differentiation (specifically, Taylor polynomials), but does not require computer algebra.
Our experiments show that Genfer is often faster than the existing exact inference tools PSI, Dice, and Prodigy.
On a range of real-world inference problems that none of these exact tools can solve, Genfer's performance is competitive with approximate Monte Carlo methods, while avoiding approximation errors.

Exact Bayesian Inference on Discrete Models via Probability Generating Functions: A Probabilistic Programming Approach
Fabian Zaiser
Department of Computer Science
University of Oxford, UK
fabian, zaiserêca, ox ax. uk
Andrzej S. Murawski
C.-H. Luke Ong
Department of Computer Science University of Oxford, UK andrzej www.wwkkiieeae...mmmm dca. ackukk
School of Computer Science and Engineering Nanyang Technological University, Singapore
luke. ong@ntu. edu. ag
Abstract
We present an exact Bayesian inference method for discrete statistical models, which can find exact solutions large class of discrete inference problems, even with infinite support and continuous priors. To express such models, we introduce probabilistic pro- gramming language that supports discrete and continuous sampling, discrete observations, affine functions, (stochastic) branching, and conditioning on discrete events. Our key tool is probubility generating functions: they provide a compact closed-form representation of distributions that are definable by programs, thus enabling the exact computation of posterior probabilities, expectation, variance, and higher moments. Our inference method is provably correct and fully automated in tool called Genfer, which uses automatic differentiation (specifically, Taylor polynomials), but does not require computer algebra. Our experiments show that Genfer is often faster than the existing exact inference tools PSI, Dice, and Prodigy. On a range real-world inference problems that none of these exact tools solve, Genfer's performance is competitive with approximate Monte Carlo methods, while avoiding approximation errors.
1 Introduction
Bayesian statistics is a highly successful framework for reasoning under uncertainty that has found widespread use in variety of fields, such as Al/machine learning, medicine and healthcare, finance and risk management, social sciences, climate science, astrophysics, and many other disciplines. At its core is the idea of representing uncertainty probability, and updating prior belief's based on observed data via Bayes' law, to arrive at posterior beliefs. A key challenge in Bayesian statistics is computing this posterior distribution: analytical solutions usually impossible or intractable, which necessitates the use of approximate methods, such as Markov-chait Monte Carlo (MCMC) variational inference. In this work, we identify large class of discrete models for which exact inference in fact possible, in particular time series models of count data, such as autoregressive models, hidden Markov models, switchpoint models. We achieve this by leveraging probability generating functions (GFs) as a representation of distributions. The GF random variable X defined to be the function G(x) == E In probability theory, it is well-known tool study random variables and their distributions, related to the moment generating and the characteristic functions [14, Chapter 4]. In computer science, GF's have previously been used for the analysis of probabilistic programs [17, 4] and for exact inference on certain classes of graphical models [26, 27].
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Here we apply them uniformly in much more general context via probabilistic programming, enabling exact inference on more expressive Bayesian models.
We characterize the class of supported models with the help of a prohabilistic programming language. Probabilistic programming [25] has recently emerged as powerful tool in Bayesian inference. Probabilistic programming systems allow users to specify complex statistical models as programs in precise yet intuitive way, and automate the Bayesian inference task. This allows practitioners to focus on modelling, leaving the development of general-purpose inference algorithms to experts. Consequently, probabilistic programming systems such as Stan [3] enjoy increasing popularity among statisticians and data scientists. We describe a programming language, called SGCL (statistical guarded command language), that extends pGCL (probabilistic GCL) [17]. This language carefully designed to be simple yet expressive, and just restrictive enough to enable exact Bayesian inference on all programs that can be written in it.
Contributions We provide a new framework for exact inference on discrete Bayesian models.
(a) Our method is applicable to a large class of discrete models with infinite support, in particular time series models of count data, such as autoregressive models for population dynamics, Bayesian switchpoint models, mixture models, and hidden Markov models. To our knowledge, no exact inference method for all but the simplest population models was known before.
(b) The models are specified in probabilistic programming language (PPL), which provides flexibility in model definition, thus facilitating model construction. Our PPL supports stochastic branching, continuous and discrete priors, discrete observations, and conditioning events involving discrete variables.
(c) Every program written in the language can be translated to generating function that represents the posterior distribution in an automatic and provably correct way. From this generating function, one can extract posterior mean, variance, and higher moments, well as the posterior probability masses (for a discrete distribution).
(d) We have built an optimized lool, called Genfer ("GENerating Functions for inFERence"), that takes a probabilistic program as input and automatically computes the aforementioned set of descriptive statisties about the posterior distribution.
(c) We demonstrate that (I) on benchmarks with finite support, Genfer's performance is often better than existing exact inference tools, and (2) on a range of real-world examples that no existing exact tool supports, Genfer is competitive with approximate Monte Carlo methods, while achieving zero approximation error.
Computing the exact posterior distribution of probabilistic programs is intractable in general requires analytical solutions to integrals [9]. For this reason, existing systems either restrict the programming language to allow only tractable constructs (this is our approach) or cannot guarantee successful inference. In the former category are Dice [16], which only supports finite discrete distributions, and SPPL [22]. which supports some infinite-support distributions but requires finite discrete priors. Both are based on probabilistic circuits [5] or extensions thereof, which allow for efficient and exact inference. In the latter category are the systems PSI [9] and Hakaru [21], which do not impose syntactic restrictions. They rely on computer algebra techniques to find closed-form solution for the posterior. Such form need not exist in general and, even if it does, the system may fail to find it and the running time unpredictable and unscalable [22]. None of the case studies featured in our evaluation (Section 5.2) can be handled by the aforementioned tools.
Related Work
Probability generating functions are a useful tool in probability theory to study random variables with infinite support, e.g. in the context of branching processes [7. Chapter 12). In the context of Bayesian inference, probabilistic generating circuits (PGCs) leverage them boost the expressiveness of probabilistic circuits [28] and enable efficient inference [15]. However, PGCs cannot handle random variables with infinite support, which the focus of our approach. Generating functions have also been applied to probabilistic graphical models: Winner and Sheldon [26] find symbolic representation of the generating function for a Poisson autoregressive model and extract posterior probabilities from it. Subsequently, Winner et al. [27] extend this model to latent variable distributions other than the Poisson distribution, where symbolic manipulation is no longer tractable. Instead they evaluate generating functions and their derivatives using automatic differentiation, which enables exact inference for graphical models. Probabilistic programming is an elegant way of generalizing graphical models, allowing much richer representation of models [10, 25, 12]. Our contribution here is new framework for exact inference on Bayesian models via probabilistic programming.
2

---

In the context of discrete probabilistic programs without conditioning, Klinkenberg et al. [17] use generating functions to (manually) analyze loop invariants and determine termination probabilities. In follow-up work, Chen et al. [4] extend these techniques to automatically check (under certain restrictions) whether a looping program generates a specified distribution. Their analysis tool Prodigy recently gained the ability to perform Bayesian inference as well [18). It supports discrete distributions with infinite support (but no continuous priors) and is less scalable than our approach (see Section 5) since it relies computer algebra.
Limitations Exact posterior inference already PSPACE-hard for probabilistic programs involv- ing only finite discrete distributions [16, Section 6]. It follows that our method cannot always be performant and has to restrict the supported class of probabilistic programs. Indeed, our programming language forbids certain constructs, such as nonlinear transformations and observations of continuous variables, in order to preserve closed-form representation of the generating functions (Section 2). For the same reason, our method cannot compute probability density functions for continuous param- eters (but probability masses for discrete parameters and exact moments for all parameters are fine). Regarding performance, the running time of inference is polynomial in the numbers observed in the program and exponential in the number of program variables (Section 4). Despite these limitations, our approach shows that exact inference is possible in the first place, and our evaluation (Section 5) demonstrates support for many real- world models and efficient exact inference in practice.
Bayesian Probabilistic Programming
Probabilistic programming languages extend ordinary programming languages with two additional constructs: one for sampling from probability distributions and one for conditioning on observed values. We first discuss simple program that can be written in our language using simplified example based on population ecology (cf. [26]).
Suppose you're re a biologist trying to estimate the size of an animal population migrating into new habitat. The immigration of animals is often modeled using a Poisson distribution. You cannot count the animals exhaustively (otherwise we wouldn't It need estimation techniques), we assume that each individual is observed independently with a certain probability; in other words, the count is binomially distributed. For simplicity, we assume that the rate of the Poisson and the probability of the binomial distribution are known, say 20 and 0.1. (For more realistic population models, see Section 5.) As a generative model, this would be written X Poisson (20); Y ~ Binomial(X.. (X 0.1).
Suppose you observe Y = animals. The Bayesian inference problem to compute the posterior
distribution P(X- = x Y = 2) given by Bayes' rule as P(X = x Y = 2) P(X=x) P(Y P(Y=2|X X=x)
2( =2)
where P(X- x) is called the prior probability, P(Y 2 X x) the likelihood and IP(Y 2) the evidence or normalization constant.
Example 2.1. In our probabilistic programming language, this simplified population model would be expressed as:
X ~ Poisson Y Binomial(X, 0.1): observe Y = 2;
The syntax looks similar to generative model notation except that the observations are expressible as a command. Such a program denotes a joint probability distribution of its program variables, which are viewed as random variables. The program statements manipulate this joint distribution. After the first two sampling statements, the distribution has the probability mass function (PMF) pi (1( y) IP[X- x, Y y = Poisson(a; 20)- Binomiia (y;; x 0.1). Observing Y restricts this to the PMF (2 y) = P|X = x, Y = y = 2]. which equals IP = x. Y = 2 p1 (x. 2) if = 2 and 0 otherwise. Note that this is not a probability, but subprobability, distribution because the total mass is less than 1. So. as final step, we need to normalize, i.e. rescale the subprobability distribution back to probability distribution. This corresponds to the division by the evidence in P=(#.M)
Bayes' rule, yielding the PMF ps(2. y) = Σ Pa(a,,)) To obtain the posterior P[X = x Y = 2],
which is a distribution of the single variable X, not the joint of X and Y. we need marginalize p3 En P2(x.0)
and find E, (s y) Ex.y (2(p)) (1,,2,
(1.2) Fix=a.y=2 Pp FY=2
Y 2]. desired.
Programming constructs Next we describe our probabilistic programming language more for- mally. It is based the probabilistic guarded command language (pGCL) from [17] but augments
3

---

it with statistical features like conditioning on events and normalization, which is why we call it Statistical GCL (SGCL). Each program operates on a fixed set of variables X = {X,,.... x, tak- ing values in Rg A program consists of a list of statements Pi; ..;...; ; Pm. The simplest statement is skip, which does nothing and useful in conditionals (like pass in Python). Variables can be transformed using affine maps, e.g. X2 = 2X, + 7X3 + 2 (note that the coefficients must be non- negative to preserve nonnegativity of the variables). Programs can branch on the value of a (discrete) variable (e.g. if X E 31,3. 7} {...} else (...)), and sample new values for variables from distribu- tions (e.g. X ~ Binomial(10, (10,, 0.5) or Xx ~ Binomial(X,, 0.5)). The supported distributions are Bernoulli, Categorical, Binomial, Uniform (both discrete and continuous), NegBinomial, Geometric, Poisson, Exponential, and Gamma. They need to have constant parameters except for the compound distributions Binomial(X, p). NegBinomial(X, p). (oisson() ( X), and Bernoulli (X), where X can be a variable. One can observe events involving discrete variables (e.g. observe Xx E (3. 5)) and values from (discrete) distributions directly (e.g. observe 3 Binomial(X,, 0.5)). Note that observe m ~ D can be seen as a convenient abbreviation of Y ~ D; observe Y = m with a fresh variable Y. After observation, the variable distribution is usually not probability distribution anymore, but a subprobability distribution ("the numerator of Bayes' rule").
Syntax In summary, the syntax of programs P has the following BNF grammar:
::= skip P1: X2 == a1X1 + ................................ + an, + Xx E A (Pt) else {P2}
Xx ~ D Xx D(X;) observe Xs E A observe m D observe D(X,)
where P. P1 and P2 are subprograms; aj..... On C E R20: m E N: A a finite subset of the naturals; and D is a supported distribution. To reduce the need for temporary variables, the abbreviations t +=, and m ~ D {...) else else } are available (see Section 4). A fully formal description of the language constructs can be found in Appendix A.
Restrictions Our language imposes several syntactic restrictions to guarantee that the generating function of any program admits a closed form, which enables exact inference (see Section 3.1): (a) only affine functions are supported (e.g. no X2 or exp(X)), (b) only comparisons between variables and constants are supported (e.g. test for equality X Y), (c) only observations from discrete distributions on and only comparisons of such random variables are supported (e.g. no observe 1.5 ~ Normal(O, 1)) I)). (d) particular choice of distributions and their composites is supported, (c) loops or recursion are not supported. A discussion of possible language extensions and relaxations of these restrictions can be found in Appendix B.3.
3 Generating Functions
Consider Example 2.1. Even though the example is an elementary exercise in probability, it is challenging to compute the normalizing constant, because one needs to evaluate infinite sum: P[Y = Emen PLY = 2| = x) P[X = x] = =xren (5)0.12 0..g--. 20 20 turns out to be 2e (see Example 3.1), but is unclear how to arrive at this result in an automated way. If X had been continuous variable, would even have to evaluate an integral. We will present technique to compute such posteriors mechanically. relies on probability generating functions, whose definition includes an infinite sum integral, but which often admit a closed form, thus enabling the exact computation of posteriors.
z!
Definition Probability generating functions are well-known tool in probability theory to study random variables and their distributions, especially discrete ones. The probability generating function of random variable X is defined to be the function G(r)) == E[r-XX For discrete random variable supported on N, this can also be written as a power series G(z) = EmEn IP[X = n] x=. Since we often deal with subprobability distributions,
omit "probability" from the name and refer to G(x) = Ex~w[== [ where j is a subprobability measure, simply as generating function (GF). For continuous variables, it is often called factorial moment generating function, but we stick with the former name in all contexts. We will use the notation (f)) or (X)) for the GF of μ or X. Note that for discrete random variables supported on N, the GF is always defined on 1-, 1]. and for continuous ones at x = 1, but it need not be defined at other x.
Probability masses and moments Many common distributions admit closed form form (see Table 1). In such cases, GFs are compact representation of distribution, even if has infinite or continuous support like the Poisson or Exponential distributions, respectively. Crucially, one can
for their GF

---

Table 1: GFs for common distributions with constant (left) and random variable parameters (right)
Distribution D
gf(X ~ (x)
Distribution D(Y)
gf(X (Z)
Binomial(n, p)
(pz + p)"
Binomial(Y,
gf(() - + pa)
Geometric(p)
NegBinomial(Y,12
gf(Y) [-[1-pp=
1-11-ppz
(oisson(A) Exponential(A)
e((=11)
Poisson(A Y)
gf(Y) (
A-log
Bernoulli(Y)
1 (æ - 1) ())(
extract probability masses and moments of a distribution from its generating function. For discrete random variables X. PIx n.] is the n-th coefficient in the power series representation of G. so can be computed the Taylor coefficient O: + (()) (0) (hence the name "probability generating
function"). For a discrete or continuous random variable X, its expected value is E (1), and more generally, its n-th factorial moment is EX - 1) (X n + 1)] E
d
de
x=
((m)(1) (1) (hence the name "fuctorial moment generating function"). The raw and central moments can easily be computed from the factorial moments. For instance, the variance is V[X] = G" (1) (1)++"(1)-
G" (1)². We will exploit these properties of generating functions through automatic differentiation.
Multivariate case The definition of GFs is extended to multidimensional distributions in straightforward way: for random variables X = (X....., Xn), their GF is the function G(z)) ==
... -rx ....r We generally follow the convention of using uppercase letters for random and program variables, and lowercase letters for E[a*] where we write æ == (Ii,..., ...) and xx == of
the corresponding parameters of the generating function. Marginalization can also be expressed in i.c. to marginalize out X,, one simply substitutes 1 for 2n in G: G In-1) = G( In Fn- 1, 1). This allows us to compute probability masses and moments of random variable in a joint distribution: we marginalize out all the other variables and then use the previous properties of the derivatives.
terms of generating functions: to obtain the GF G of the joint distribution of (
3.1 Translating programs to generating functions
The standard way of describing the meaning of probabilistic programs is assigning (sub-Jprobability distributions them. An influential example Kozen's distribution transformer semantics [19], where each program statement transforms the joint distribution of all the variables X (Xi...., Xn). Since Kozen's language does not include observations, we present the full semantics in Appendix A. We call this the standard semantics of a probabilistic program and write [P]saa() for the transformation of y by the program P. As a last step, the subprobability distribution f has to be normalized, which we write normalize(1) (4) J This reduces the Bayesian inference problem for
a given program to computing its semantics, starting from the joint distribution Dirac(O,,) in which all variables are initialized to 0 with probability 1. While mathematically useful, the distribution transformer semantics is hardly amenable to computation as it involves integrals and infinite sums.
Instead, we shall compute the generating function of the posterior distribution represented by the probabilistic program. Then we can extract posterior probability masses and moments using automatic differentiation. Each statement in the programming language transforms the generating function of the distribution of program states, i.e. the joint distribution of the values of the variables X = (X,...., X.). Initially, we start with the constant function G = 1, which corresponds to all variables being initialized with since E E[a-| = 1.
The generating function semantics of program [P]gr describes how to transform G to the generating function [P] (G) for the distribution at the end. is defined in Table 2, where the update notation az(i - a) denotes (#1...., T4-1, a, T441s-..,
and 1,, means the n-tuple (1...., 1). The first five rules were already described in [4], so we only explain them briefly: skip leaves everything unchanged and Pi; P22 chains two statements by transforming with [P]]]et and then [Po] To explain linear
assignments, consider the case of only two variables: X, = 2X + 3X2 + 5. Then
(G) B[x +3X2 +3X,
5 22 = =)B(( (xfo2) X²] 2t G(u
For conditionals if X3 E A (P.) else (Pz). we split the generating function G into two parts: one where the condition is satisfied (G's e and its complement (G- GG. The former is transformed
5

---

Table Generating function semantics of programming constructs
Language construct P
skip
G(æ)
PiP
[PaIor ([
X& = a a +
rf G(æ') where 24. Fg and 2 == 2,21 for i
Xx A ( (12}
+ [Palee GxyEa)
where Gx,ea(x) = EiEA
of G( x2
Xx ~
G(x + (
X2 ((X,)
1+ G( for D (Binomial( (p), NegBinomiall p), Poisson(X- -)}
1,5 Ijj (zee)
G(æ(k - 11) 2z (3 - a. 1I) for D Bernoulli(
observe X3 E A Normalization
Gx = Lieg
normalize G) G(.
by the then-branch [Pi]gr. the latter by the else-branch [P2lgt. The computation of Gx is best understood by thinking of G as a power series where we keep only the terms where the exponent of rg is in A. Sampling X ~ D) from distribution with constant parameters works by first marginalizing out X3 and then multiplying by the generating function of D with parameter rY.-
The first new construct is sampling X ~ D(X,) from compound distributions (see Appendix B for a detailed explanation). Observing events X2 E À uses GxLE like conditionals, as explained above. Just like the subprobability distribution defined by program has to be normalized as a last step, we have normalize the generating function. The normalizing constant is calculated by marginalizing out all variables: G(I,.... 1). So we obtain the generating function representing the normalized G These intuitions can be
posterior distribution by rescaling with the inverse: normalize(G) == c(i.) made rigorous in the form of the following theorem, which is proven in Appendix B.2.
Theorem 3.1. The GF semantics is correct wrr. the standard semantics: for any SGCL program P and subprobability distribution il on R2o we have [ (g)) = stt((ssta(p))) In particular, correctly computes the GF of the posterior distribution of P as normalize([ ([PP (1)). Furthernore, there is some R > 1 such that [P]st(]) and [ormalize[[ (1)) are defined on far E RR" 1 < R] where Q₂ = R if the variable X, is supported on N and Q2 otherwise.
I; <
Novelty The semantics builds upon [17, 4]. To our knowledge, the GF semantics of the compound distributions Poisson(X (X,) and Bernoull((X,) is novel, and the former is required to support most models in Section 5. While the GF of observations has been considered in the context of a specific model [26], this has not been done the general context of a probabilistic programming language before. More generally, previous works involving GFs only considered discrete distributions, whereas we also allow sampling from continuous distributions. This is a major generalization and requires different proof techniques because the power series representation Eien IP[X i|æ*,, on which the proofs in [26, 17, 4] rely, is not valid for continuous distributions.
We can find the posterior distribution mechan- ically by applying the rules from the GF semanties. We start with the GF A(z, y) []10/"]== 1 corresponding to X and Y being initialized to O. Sampling X changes this to GF B ( =
Example 3.1 (GF translation). Consider Example 2.1.
A(1, (00,))c2((a-) (a-1)) = 1)) Sampling Y yields C(x, 8) = ((x(0.1y + 0.9), 1) = 22(y+9)-20
od Observing Y 2 yields D(I, 1 dyz Oxz 0) 2 18x---2020--- To normalize, divide by D(a.y) D(1,1) 2c -2, obtaining F(x, y) D(I.I) 2 since A(x, 1.
As described above, 10 individuals P(X- 10] 10t 10 E 1) = 991796451840e- and the expected value of the posterior E[X] 2 de (1, 1) = 20.
extract from this GF the posterior probability of, for example, exactly 0
4 Implementation & Optimizations
The main difficulty in implementing the GF semantics is the computation the partial derivatives. natural approach (as followed by [4, 18)) is manipulate symbolic representations of the generating

---

functions and to use computer algebra for the derivatives. However this usually scales badly, as demonstrated by Winner et al. [27]. because the size of the generating functions usually grows quickly with the data conditioned on. To see why, note that every observe Xx= = d statement in the program is translated to d-th partial derivative. Since probabilistic programs tend to contain many data points, it is common for the total order of derivatives to be in the hundreds. The size of the symbolic representation of a function can (and typically does) grow exponentially in the order of the derivative: the derivative of the product of two functions f g is the sum of two products fl gg + f-g'. so the representation doubles in size. Hence the running time would be 2224)) where d is the sum of all observed values, which is clearly unacceptable.
Instead, we exploit the fact that we do not need to generate the full representation of a GF, but merely to evaluate it and its derivatives. We implement our own automatic differentiation framework for this because existing ones are not designed for computing derivatives of order greater than, say, 4 or 5. In fact, it is more efficient to work with Taylor polynomials instead of higher derivatives directly. Winner et al. [27] already do this for the population model (with only one variable), and we extend this to our more general setting with multiple variables, requiring multivariate Taylor polynomials. In this approach, derivatives are the easy part as they can be read off the Taylor coefficients, but the composition of Taylor polynomials is the bottleneck. Winner et al. [27] use naive Old)) approach, which is fast enough for their single-variable use case.
Running time For n variables, naive composition of Taylor polynomials takes ((dan) time, where d is the sum of all observations in the program, i.c. the total order of differentiation, i.c. the degree of the polynomial. Note that this is polynomial in d. contrary to the symbolic approach, but exponential in the number of variables n. This is not as bad as it seems because in many cases, the number of variables can be kept to one or two, as opposed to the values of data points (such as the models from Section 5). In fact, we exploit the specific composition structure of generating functions to achieve O(d)) for = and (dd++3) for > 2 in the worst case, while often being faster in practice. Overall, our implementation takes (ss time in the worst case, where is the number of statements in the program, d is the sum of all observed values, and n is the number of program variables (see Appendix C.3).
Reducing the number of variables Given the exponential running time in n. it is crucial to reduce the number of variables when writing probabilistic programs. For one thing, program variables that are no longer needed can often be reused for a different purpose later. Furthermore, assignment and sampling can be combined with addition: Xx +=... stands for Xn+1 := ::...; X := X + +X,+1 and
Xx. +~ D for X,+1 X == X2 + Xn-1. The GF' for these statements can easily be computed without introducing the temporary variable X Similarly, in observe statements and if-conditions,
we use the shorthand ~ D with E NN for the event X.+1 m where X,+1 ~ D. Probabilistic programs typically contain many such observations, in particular from compound distributions. Hence it is worthwhile to optimize the generating function to avoid this extra variable Xn+1- Winner and Sheldon [26] can avoid an extra variable for a compound binomial distribution in the context of a specific model. We extend this to our more general setting with continuous variables, and also present optimized semantics for compound Poisson, negative binomial, and Bernoulli distributions. In fact, this optimization is essential to achieving good performance for many of the examples in Section 5. The optimized translation and its correctness proof can be found in Appendix C.2.
Implementation Our tool Genfer reads program file and outputs the posterior mean, variance, skewness, and kurtosis of specified variable. For discrete variables supported on N, it also computes the posterior probability masses up to configurable threshold. To have tight control over perfor- mance, especially for the multivariate Taylor polynomials, Genfer is written in Rust [20]. a safe sys- tems programming language. Our implementation is available on GitHub: ctthub.coo//zzaierrgeettt
Numerical issues Genfer can use several number formats for its computations: 64-bit floating point (the default and fastest), floating point with user-specified precision, and rational numbers (if no irrational numbers occur). To ensure that the floating-port point results are numerically stable, we also implemented interval arithmetic to bound the rounding errors. Initially, we found that programs with continuous distributions led to catastrophic cancellation errors, due to the logarithmic term in their GFs, whose Taylor expansion is badly behaved. We fixed this problem with a slightly modified representation of the GFs, avoiding the logarithms (details in Appendix C). For all the examples in Section 5.2, our results are accurate up to at least 5 significant digits.
7

---

Table 3: Comparison of inference times of tools for exact inference on PSI's benchmarks [9]
Tool
Genfer (FP)
Dice (FP)
Genfer (Q)
Dice (Q)
Prodigy
PSI
alarm (F)
0.0005s
0.0067s
0.0012s
0.0066s
0.011s
0.0053s
clickGrrah (C)
0.11s
unsupported
3.4s
unsupported
unsupported
46s
clinicalTriill (C)
150s
unsupported
1117s
unsupported
unsupported
timeout
clinicalTrrial2 (C)
0.0024s
unsupported
0.031s
unsupported
unsupported
0.46s
digitRecognition (F)
0.021s
0.83s
0.11s
2.7s
31s
146s
evidencel (F)
0.0002s
0.0057s
0.0003s
0.0056s
0.0030s
0.0016s
evidence2 (F)
0.0002s
0.0056s
0.0004s
0.0057s
0,0032s
0.0018s
grass (F)
0.0008s
0.0067s
0.0044s
0.0067s
0.019s
0.014s
murderMystery (F)
0.0002s
0.0055s
0.0003s
0.0057s
0.0028s
0.0021s
noisyOr (F)
0.0016s
0.0085s
0.019s
0.0088s
0.21s
0.055s
twoCoins (F)
0.0002s
0.0054s
0.0003s
0.0057s
0.0032s
0.0017s
5 Empirical Evaluation
5.1 Comparison with exact inference methods
We compare our tool Genfer with the following tools for exact Bayesian inference: Dice [16], which uses weighted model counting, PSI [9], which manipulates density functions using computer algehra, and Prodigy [18], which is based on generating functions like Genfer, but uses computer algebra instead of automatic differentiation. We evaluate them on the PSI benchmarks [9], excluding those that only PSI supports, e.g. due to observations from continuous distributions. Most benchmarks only use finite discrete distributions (labeled "F"), but three feature continuous priors (labeled "C").
We measured the wall-clock inference time for each tool, excluding startup time and input file parsing, and recorded the minimum from 5 consecutive runs with one-hour timeout." Dice and Genfer default to floating-point (FP) numbers, whereas Prodigy and PSI use rational numbers, which is slower but prevents rounding errors. For a fair comparison, we evaluated tools in rational mode and separately compared Dice with Genfer in FP mode. The results (Table 3) demonstrate Genfer's speed even on finite discrete models, despite our primary focus on models with infinite support.
5.2 Comparison with approximate inference methods
It is impossible to compare our approach with other exact inference methods realistic models with infinite support (Section 5.3): the scalable systems Dice [16] and SPPL [22] don't 't support such priors and the symbolic solvers PSI [9] and Prodigy [18] run out of memory time out after an hour.
Truncation As an alternative, we considered approximating the posterior by truncating discrete distributions with infinite support. This reduces the problem to finite discrete inference, which is more amenable to exact techniques. We decided against this approach because Winner and Sheldon [26] already demonstrated its inferiority to GF-based exact inference their graphical model. Moreover, it is harder to truncate general probabilistic programs, and even impossible for continuous priors.
Monte-Carla inference Hence, we compare our approach with Monte Carlo inference methods. Specifically, we choose the Anglican [24] probabilistic programming system because it offers the best built-in support for discrete models with many sfate-of-the-an art inference algorithms. Other popular systems are less suitable: Gen [6] specializes in programmable inference; Stan [3], Turing [8], and Pyro [1] mainly target continuous models; and WebPPL's [11] discrete inference algorithms are less extensive than Anglican's (e.g. no support for interacting particle MCMC).
Methodology The approximation error of a Monte Carlo inference algorithm depends on its settings' (e.g. the number of particles for SMC) and decreases with the number of samples the longer it is run. To ensure fair comparison, we use the following setup: for each inference problem, we run several inference algorithms with various configurations (settings and sampling budgets) and
'We were not able to run Hakaru [21], which uses the computer algebra system Maple as backend, but we do not expect it to produce better results than the other systems, given that PSI is generally more performant [9]. 2PSI was run with the experimental flag d-dp on finite discrete henchmarks for improved performance. "ee /tposs//poobprgg gi thub 1o/ hooanglia can/ weaaaiinfereet inference/ for a full list of Anglican's 's options.
8

---

0.35
ours
0.4
+
6.25
0.30
PMCNC
Eite (700)
0.3
0.3
IMMEMC
6.70
7Gtbb 2011
AMH
0.2-
RN-
c.35
SMC
SMC
6.30
0.1
6.05
0.00
0.0
3.5 Time (i
2.0
25
2 Time Iin k)
Time On k)
(a) Original model
(b) Modified arrival rate
(c) Two-type branching process
Figure 1:
Comparison of the population model and its modifications with approximate inference: mean and standard error of the computation time and TVD over 20 repeated runs.
report the approximation error and the elapsed time. To measure the quality of the approximation, report the total variation distance (TVD) between the exact solution and the approximate posterior distribution, as well as the approximation errors of the posterior mean μ., standard deviation o, skewness (third standardized moment) S, and kurtosis K. To ensure invariance of our error metrics under translation and scaling of the distribution, we compute the error of the mean as الرμا where f the true posterior mean, á its approximation, and a the true posterior standard deviation; the relative error of the standard deviation (because is invariant under translation but not scaling), and the absolute error of skewness and kurtosis (because they are invariant under both translation and scaling). To reduce noise, we average all these error measures and the computation times over 20 runs and report the standard error as error bars. We run several well-known inference algorithms Eetropolis-Hastingg (LMH),
implemented in Anglican: importance sampling (IS). Lightweight Random Walk Metropolis- Hastings (RMH), Sequential Monte Carlo (SMC), Partiele Gibbs (PGibbs), and interacting particle MCMC (IPMCMC). For comparability, in all experiments, each algorithm is run with two sampling budgets and, if possible, two different settings (one being the defaults) for a total of four configurations. The sampling budgets were 1000 or 10000, because significantly lower sample sizes gave unusable results and significantly higher sample sizes took much more time than exact method. We discard the first 20% of the samples, standard procedure called "burn-in".
Note that this setup is generous to the approximate methods because we only report the average time for one run of each configuration. However, in practice, one does not know the best configuration, so the algorithms need to run several times with different settings. By contrast, our method requires only one run because the result is exact.
5.3 Benchmarks with infinite-support distributions
Population ecology Our first benchmark comes from [26, 27] and models animal populations. We have seen a simplified version in Example 2.1. Here we model a population N₂ at time steps k O,..., m. At each time step, there is Poisson-distributee number of new arrivals, which are added to the binomially distributed number of survivors from the previous time step. Each individual
observed with fixed probability so the number of observed individuals is binomially distributed:
New ~ Poisson(Ax); Survivors: ~ Binomial(Nix11.))
N₂- = New + SurvivorsHi observe yk ~ Binomial(Nee,
where the model parameters λ E R, 8 E 1] are taken from [26]: the detection probability is set to 0.2 ([26) considers a range of values, but we pick one for space reasons); and the observed number ya E N of individuals in the population at time step l is simulated from the same ground truth as [26]. The goal is to infer the final number Nn of individuals in the population. We set the population size model parameter and hence the observed values which influence the running time to be 4 times larger than the largest in [26] (see details in Appendix D.3). The results (Fig. Ia) show that our method is superior to MCMC methods in both computation time and accuracy since it is exact.
Modifications While this model was already solved exactly in [26], our probabilistic programming approach makes it trivial to modify the model, since the whole inference is automated and one only needs to change a few lines of the program: (a) we can model the possibility natural disasters af- fecting the offspring rate with conditional: Disaster ~ (ernoulli(0.1)) if Disaster = 1 1 (News. ~ Poisson(X') else [New; ~ Poisson( ())). or (b) instead of a single population, we can model popu- lations of two kinds of individuals that interact, i.c. multitype branching process (see Fig. Ic). None

---

02
0.53
PM Nokss 24C
32x
cons
0003
60
32
0063
02
(a) Switchpoint: TVD error
(b) Mixture: TVD error
(c) Mixture: histogram
(d) HMM: TVD error
Figure 2: Plots for the switchpoint, mixture, and hidden Markov model.
4-1
un
1M4 SMC
(a) Error of the mean
(b) Error of the std. dev. (c) Error of the skewness
(d) Error of the kurtosis
Figure 3: Comparison of moments for the switchpoint model: approximate inference vs our method.
of these modifications can be handled by [26] or [27]. The results of the first modification (Fig. 1b) are very similar. The more complex second modification takes longer to solve exactly, but less time than approximate inference with 10000 samples, and achieves zero error.
Switchpoint model Our second benchmark Bayesian switchpoint analysis, which is about detecting a change in the frequency of certain events over time. We use the model from [23] with continuous priors and its 11l real-world data points about the frequency of coal-mining accidents. We compare both the moment errors (Fig. 3) and the TVD errors (Fig. 2a). In both cases, the approximations are less accurate and take longer than our exact method.
Mixture model We consider binary mixture model on the same data set, with equal mixture weights and a geometric prior for the rates: each data point is observed from a mixture of two Poisson distributions with different rates and the task is to infer these rates. Due to their multimodality, mixture models are notoriously hard for approximate inference methods, which confirmed in Fig. 2b. Even the runs with the lowest error cover only one of the two modes (ef. the sample histogram in Fig. 2c).
Hidden Markov model We a hidden Markov model based on [22, Section 2.2], but involving infinite (geometric) priors. It is a two-state system with known transition probabilities and the rate for the observed data depends on the hidden state. We run this model on 30 simulated data points. For this model as well, our method clearly outperforms approximate methods (Fig. 2d).
To our knowledge, our approach is the first to find exact solution these problems, except the very first problem without the modifications, which appeared in [26]. For brevity, we only presented the most important aspects of these benchmarks, relegating their encoding probabilistic programs to Appendix D.3. Code and reproduction instructions are provided in the supplementary material.
Conclusion
By leveraging generating functions, we have developed and proven correct a framework for exact Bayesian inference on discrete models, even with infinite support and continuous priors. We have demonstrated competitive performance on a range of models specified in an expressive probabilistic programming language, which our tool Genfer processes automatically.
Future work It is natural question how our method could be integrated with more general probabilistic programming system. For example, in a sampling-base based inference algorithm, one could imagine using generating functions to solve subprograms exactly if this is possible. More generally, it would be desirable to explore how the compositionality of the GF translation can be improved. As it stands, our GF translation describes the joint distribution of all program variables we never reason "locally" about a subset of the variables. A compositional approach would likely facilitate the application of the GF method functional probabilistic languages like Anglican.
10

---

Acknowledgments and Disclosure of Funding
We would like to thank Maria Craciun, Hugo Paquet, Tim Reichelt, and Dominik Wagner for providing valuable input on this work. We are especially grateful to Mathieu Huot for insightful discussions and very helpful feedback on a draft of this paper.
This research was supported by the Engineering and Physical Sciences Research Council (studentship 2285273. grant EP/T006579) and the National Research Foundation, Singapore, under its RSS Scheme (NRF-RSS2022-009) For the puppose Open Access, the authors have applied a CC BY public copyright license any Author Accepted Manuscript version arising from this submission.
References
[1] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, A. Szerlip, P. Horsfall, and N. D. Goodman. Pyro: Deep universal probabilistic programming. J. Mach. Learn. Res., 20:28:1 1-2:66, 2019, URL http::/jmll org/ http:/pppraa///////1 10-200//88433... html.
[2] R. P. Brent and H. T. Kung. Fast algorithms for manipulating formal power series. J. ACM, 25(4): 581-595, 1978. doi: 10.1 00/144/422292232200999 URL attps:://oiis /ttp::////oo.o///////////// .14 1145/322092. 322099.
[3] B. Carpenter, A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. Brubaker, J. Guo, P. Li, and A. Riddell. Stan A probabilistic programming language. Journal of Statistical Software, 76(1), 2017. ISSN 1548-7660. doi: 10.186377jss.w076..01 URL https: //tp:://wwww... osti htt.l//ooobbboob/io// 14302022
[4] M. Chen, J. Katoen, I.. Klinkenherg, and T. Winkler. Does a program yield the right distribution? verifying probabilistic programs via generating functions. In Shoham and Y. Vizel, editors, Computer Aided Verification 34th International Conference, CAV 2022, Haifa, Israel, August 7-10, 2022, Proceedings, Part I, volume 13371 of Lecture Notes in Computer Science, pages URL https://doi.org//
79-101. Springer, 2022. doi: 10. 1007 3-08---03-----00 -31--1111111155
[5] Y. Choi, A. Vergari, and Van den Broeck. Probabilistic circuits: A unifying framework for tractable probabilistic models, 2020. URL ProbCirc20, pdf. Accessed: 2023-10-22.
cs. ucla. edue adu/papers/
[6] M. F. Jusumano-Townee F. A. Saad, A. K. Lew, and V. K. Mansinghka. Gen: general- purpose probabilistic programming system with programmable inference. ACM SIGPLAN International Conference on Programming Language Design and Implementation, PLDI 2019.
ACM, 2019. doi: 10. 1221144/3312211333323
[7] W. Feller. An introduction to probability theory and its applications, volume I. Wiley Series in Probability and Statistics. John Wiley & Sons, Nashville, TN, 3 edition, jan 1968.
[8] H. Ge, K. Xu, and Z. Ghahramani. Turing: A language for flexible probabilistic infer- ence. In A. Storkey and F. Perez-Cruz, editors, Proceedings of the Twenty-Firs, Interna- tional Conference Machine Learning Research, pages 1682-1690. PMLR, 09-11 Apr 2018. URL https: //proceedings..o wtt.iiiinnnmm..ss//// press/ v84/
Artificial Intelligence and Statistics, volume 84 of Proceedings of
/84//e188// htt/ge18bb.... html...tm..
[9] T. Gehr, S. Misailovic, and M. T. Vechev. PSI: exact symbolic inference for probabilistic programs. In S. Chaudhuri and A. Farzan, editors, Computer Aided Verification 28th Inter- national Conference, CAV 2016, Toronto, ON, Canada, July 17-23, 2016. Proceedings, Part I. volume 9779 of Lecture Notes in Computer Science, pages 62 23-83.33 Springer, 2016. doi: 10. 1007/978-3-319-415288
URL ht //tp:///do...o http::/dooooooo////// wow......////.//////.../// /ww.://0000000------999999999990
[10] Z. Ghahramani. Probabilistic machine learning and artificial intelligence. Nat., 521(7553):452- 459, 2015. doi: 10.1038/nature 14541. URL
[11] N. D. Goodman and A. Stuhlmüller. The Design and Implementation of Probabilistic Program- ming Languages. http: /ttp://dippl.... org, 2014. Accessed: 2023-4-6.

---

[12] A. D. Gordon, T. A. Henzinger, A. V. Nori, and S. K. Rajamani. Probabilistic programming. In J. Herbsleb and M. B. Dwyer, editors, Proceedings of the on Future of Software Engineering. FOSE 2014, Hyderabad, India, May 31 June 7, 2014, pages 167-181. ACM, 2014. doi: 10.1 (411145/25938822559955
URL URL htsp:: httes:/ddooooooo http///..ooog//./// 11.000.1111 1145 2593882 2593900.
[13] A. Griewank and Walther. Evaluating derivatives principles and techniques of algorithmic differentiation, Second Edition. SIAM, 2008. ISBN 100--00-89871-659--1 doi: 10.
9780898717761. URL www.htpps///oomm.mmm http://0...13/1119999977777771.1
[14] A. Gut. Probability: Graduate Course. Springer Texts in Statistics. Springer, New York, second edition, 2013. ISBN 378-1-4614-4707-88 978-1-41777-1- doi: 10.1007/ 978- 778-1441144708-53 URL http:///.....000..00000077611161
[15] J. Harviainen, V. P. Ramaswamy, and M. Koivisto. On inference and learning with probabilistic generating circuits. In R. Evans and I. Shpitser, editors, Uncertainty in Artificial Intelligence, UAI 2023, July 31 4 August 2023, Pittsburgh, PA, USA, volume 216 of Proceedings of Machine Learning Research, pages 829-838 PMLR, 2023. URL https: httttps:/pprceediiiii mlr press/ v216 2216/harviainen22 html.
[16] S. Holtzen, G. V. den Broeck, and T. D. Millstein. Scaling exact inference for discrete probabilistic programs. Proc. ACM Program. Lang. (4000PSLA):10::1140000000 2020. doi: 10.1 10.1145/3428208. URL https: /tt..//diiiooo org/10. 10. 1145/144293 3428208
[17] L. Klinkenberg, K. Batz, B. L. Kaminski, J. Katoen, J. Moerman, and T. Winkler. Generating functions for probabilistic programs. In M. Fernández, editor, Logic-Based Based Program Synthesis and Transformation - 30th International Symposium, LOPSTR 2020, Bologna, Italy, September 7-9, 2020, Proceedings, volume 12561 Lecture Notes Computer Science, pages 231-248. Springer, 2020. doi: 10. (0000797----30 68446-4) 978- 378-3-030-1 68-46-4-44422
URL https: http::://ooo..ooooo 1007
[18] L. Klinkenberg, Blumenthal, M. Chen, and Katoen. Exact bayesian inference for loopy probabiliiticc programs. CoRR, 0txx/2377.07314 2023. doi: 10448550aaxxx.22377071444 URL
https://doi.orgg. www.o//....../////////////// 48550/ 2307 07314.
[19] D. Kozen. Semanties of probabilistic programs. J. Comput. Syst. Sci., 22(3):328-350, 1981. doi: 10.1016/0022-00008811888000 (9000366. URL ittps://doii org/10. 101100000002200000 90036-2.
[20] N. D. Matsakis and F. S. K. II. The rust language. In M. Feldman and S. Taft, editors, Proceedings of the 2014 ACM SIGAda annual conference High integrity language technology. HILT 2014, Portland, Oregon, USA, October 18-21, 2014, pages 103-104, ACM, 2014. doi: 00.1145/266317 26622663666 188.
URL tts: wttp:///ooooooo wwd.dooooorg.. .rg.ooggggoooooo /10. 1..111.. /111/1/6611111 2663188
[21] P. Narayanan, J. Carette, W. Romano, C. Shan, and R. Zinkov. Probabilistic inference by program transformation in hakaru (system description). In O. Kiselyov and A. King, editors, Functional and Logic Programminng 13th International Symposium, FLOPS 2016, Kochi, Japan, March 4-6, 2016, Proceedings, volume 9613 of Lecture Notes in Computer Science, URL https://doi. org/
pages 62-79. Springer, 2016. doi: 10. 10211007777722 -------96-9---------
[22] F. A. Saad, M. C. Rinard, and V. K. Mansinghka. SPPL: probabilistic programming with fast exact symbolic inference. S. N. Freund and E. Yahav, editors, PLDI '21: 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation, Virtual Event, Canada, June 20-25, 2021, pages 804 819. ACM, 2021. doi: 10.1145/34534833445490 URL attps://dois org org/ 10. 1145 01144///43483 3454078.
[23] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck. Probabilistic programming in python using pyme3. PeerJ Comput. Sci., 2:e55, 2016. doi: 10.7717/peerj-es.555 URL https:////ddd /doi .rrg/ 10. N0.7717/peerj-et
[24] D. Tolpin, J. van de Meent, and F. D. Wood. Probabilistic programming in anglican. In European Conference on Machine Learning and Knowledge Discovery in Databases, ECML PKDD 2015,
volume 9286 of LNCS. Springer, 2015, doi: 111007/978--331--2466
12

---

[25] J. van de Meent, B. Paige, H. Yang, and F. Wood. An introduction to probabilistic programming. CoRR. abs/1809. 2018. URL http:////arx http://awxii.orabra////aaaaas//e aba/a/a9..1.... 10756.
[26] K. Winner and D. Sheldon. Probabilistic inference with generating functions for poisson latent variable models. In D. D. Lee, M. Sugiyama, U. von Luxburg, I. Guyon, and R. Gar- nett, editors, Advances in Neural Information Processing Systems 29: Anmual Conference on Neural Information Processing Systems 2016. December 5-10, 2016, Barcelona, Spain, pages 2640-2648, 2016. URL attps://proceeddngg cc148666222677222008777444488883 www.:/eeeaaaa-aaaaateeeeaa Abstract html.
wew.iiinnnnnn.rr.eeeeeeee neurips cc/ wc/paper//ee/e/eecceee/eeeeeeeea /20///200/hha/a
[27] K. Winner, D. Sujono, and D. Sheldon. Exact inference for integer latent-variable models. In D. Precup and Y. W. Teh, editors, Proceedings of the 34th International Conference on Machine Learning, ICMI. 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine Learning Research, pages 3761-3770. PMLR, 2017. URL http: / //proceedings.no mll. prwss/ 070/winners 1 html.
[28] H. Zhang, B. Juba and G. V. den Broeck. Probabilistic generating circuits. Meila and T. Zhang. editors, Proceedings of the 38th International Conference on Machine Learning, ICMI. 2021, 18-24 July 2021, Virtual Event, volume 139 of Proceedings of Machine Learn- ing Research, pages 12447-12457. PMLR, 2021. URL http: http:///oooeeddeee ings mlr, alr.press. v139 #139/zhang2/ html.
13

---

A Details on the Probabilistic Programming Language SGCL
A.1 minimal grammar
minimal grammar of SGCL looks like this:
::= skip P1 X2 a1
an + [(P)) {P2}
| fail Xx ~ D Xx D(X,)
where 01, a2, E R20. p E [0, 1] are literals, and A is finite subset of N. The X₃ for k = l..... are variables that can take on values in Rzo- To explain, skip leaves everything unchanged, P1: P2 rums P1 and then P2. and X == a X + performs an affine transformation where a E R3, and Rz0 ensure that the support of X₃ is subset of Rzo). The if construct allows branching based on comparisons of variable X₁ with constant natural numbers, with the requirement that X₂ be discrete variable on and X D samples X₁ from a primitive distribution (with constant parameters)
D {Bernoullip)), Categorical(po. Pm Binomial(m, Uniform Unifomm NegBinomial(m.p2) (), Geometric(p), Poisson(A). Exponential(a)) Gamma(a, B), Uniform bj
1p E [0, 1], a. R20, 1,m E N, a, b E Rzo, a - b}
whereas X ~ DX,,) samples Xg from a compound distribution (with a parameter that is again variable)
D(X,) = (Binomial(X, pX, NegBinomial(X;; Poisson(A X2), Bernoulli(X,)
|pe 00, 1),j E (1......), E [2o))
A.2 Syntactic sugar
For convenience, we add some syntactic sugar, which is expanded follows:
X + a X +ew X (a(k lk a 1)) X +
+~ D ~ X ~ D; Xx += Xnt1 (Xn+1 fresh)
observe o if (skip} else {fail}
m ~ D (Pi) else [22] ~ X,+1 ~ D; if Xx+1 = m (Pi) else (P2} (Xa+1 fresh)
if Xxx (Pi) else [P2] if X2 {m) {Pi} else (P2}
-6 (P) else [P]] ~ if o (P2) else {Pt}
o A 62 (Pi) else [P2] if ci (if Φ2 (Pt else (P2) else (P2) else
if o1 V 62 (Pi) else { ~ if di [Pi] else (if Φ2 {P1} else [P2}}
A.3 Standard measure semanties
Throughout the text, we always assume that N carries the discrete d- algebra P(N) and IR the Borel r-algebra ((R). Products of these spaces (e.g. x IK) are endowed with the product -llgebrr. In the following, we write for (II,..., .r,) and use the update notation æ +
Notation and Conventions
Ii+1, xn,) In fact, we even write updates æli; j + di,..., . + ax) which stands for æ(i1 a.].... [ix + ax). Finally,
a) to denote ..,,,..,,,;.........
write 1, for the k-tuple (1, (L...., 1) and similarly for O₃- We write for the zero measure or the constant zero function. We also use Iverson brackets [o] for logical condition ø: [Ø] is 1 if o is satisfied and otherwise. For example, the indicator function 1, of a set S can be defined using Iverson brackets: 1s(x) = [a E S}. For measure y R", define (S) x E S X₂ E A)).
Standard transformer semantics The standard semantics for such a probabilistic program uses distribution transformers [19]. Since language includes conditioning. the distribution of program states will not always be probability distribution, but a subprobability distribution on RS for n variables, i.e. a function μ B(R - [o, 1] with ((R < 1. The semantics of a program P is thus measure transformer [P] Meas(R Meas (R We assume that the initial measure
14

---

Iskipl std (11) = 1
P1: 23 std (11) = [22] sta ([Pi
[ = X + sts (1) 14(X ER" x X + E ES)
lif X E A (P) else [P (4) = [PlIad + [P2lsta -x,eA)
where (X)EA(() = m((X S Xx AA)
[X ~ (1)) = / / X Y E S) dD) d [ ((x)))sa((4)()) / [X| Y
dx(x)
[tdi =
normalize(1)) = I dx(()
Figure 4: The standard measure semantics
[-]sed.
is the probability distribution where every variable X₁ is 0 with probability 1, i... X Dira The semantics is shown in Fig. 4 where X is vector of length
The first rules were essentially given by [19], with slightly different notation and presentation. Kozen defines some of the measures rectangles (S = S1 x ... x 5,) only, which uniquely extends to measure on all sets the product o-algebra. By contrast, we directly define the measure on all sets in the o-algebra using integrals.
We quickly explain the rules. The skip command does nothing, so the distribution f of the program is unchanged. Chaining commands P1; P2 has the effect of composing the corresponding trans- formations of μr. Affine transformations of program variables yield measure of S that is given by the original measure of the preimage of the event S under the affine map. For conditionals X E (P) else { we consider only outcomes satisfying X₃ E in the first branch and the complement in the second branch and then sum both possibilities. Sampling from distributions X
a bit more complicated. Note that yuS) = J[[ S du( X). To obtain the measure after sampling X2 D. we care about the new state X + Y] being in S, where Y has distribution D. This is exactly expressed by the double integral. The rule for compound distributions X₁ ~ D(X;) is essentially the same the one for distributions with constant parameters X D, except the dependency on X, needs to respected.
The two new rules (not derivable from Kozen's work) are for fail and the implicit normalization (ormalize(μ)) as final step the end of program. The rule for fail departs from Kozen's semanties in the sense that it yields a subprobability distribution (i.e. the total mass can be less than 1). Intuitively, fail "cuts off" or assigns the probability 0 certain branches in the program. This is called hard conditioning in the PPL literature. For example, if X1 = 5 (skip) else {fail) has the effect of conditioning on the event that X, 5 because all other scenarios are assigned probability
The normalization construct corresponds to the division by the evidence in Bayes' rule and scales the subprobability distribution of subprogram P back to probability distribution. The semantics 1 is Tap(Xj so the distribution f is scaled by the inverse of the total probability mass of thus normalizing it.
B Details on Generating Functions
An exhaustive list of the supported distributions and their generating functions can be found in Tables 4a and 4b.
15

---

(a) GFs for common distributions
Distribution
gf(X ~ (x)
Dirac(a)
x²
Bernoullí(p)
pa
Categorical(po...
Poo)
Σ P:t
(b) GFs for common compound distributions
Binomial(n, p)
(x +
p)
Uniform (a..b)
1+1
Distribution D(Y)
gf(X ~ (x)
[6 6--+1][11-2 -7
Binomial
gf(Y) pp px)
NegBinomial(r..p)
(1=1p) p)=
NegBinomial(Y,,)
8f(Y) 1-<1-pp)>
Geometric(p)
1-( p/ Ax
Poisson(A Y)
gf())
Poisson(A)
Bernoulli(Y)
+ (a - 1) (gf(Y))(1)
Exponential(A)
À À- log
Gamma(a; 3)
3-log 3
x +
Uniform(us (a, b]
(6-a) log
I =
(skiplsr(C) G
[P; Po (G) = [P.Ia((P)
[X = a x T
cla (G)(x) = = G(æ')
where x/. == xht and x. == x for i + k
if x E (PP) else (G)) = [Pi] ( (Gx,ea) + [Palerr Gx,ea)
and A C N finite where Gx,ed(x)= = Σ G([ + 0)) o
il
iEA
ifX, has support in N
[Xx (G) G(x(k - (D) [x ~ D(X,) (G) = G(x(k + 1,j - x; (f((()))
for D (Binomial(-,p)) Negiinmiall (-,p), Poisson(A.-))
[X ~ Bernoulli X, Mar (G) = Gz + 1)) + x 1) (,G(x + 1))
if X, has support in [0. 1]
[fail]] (G) = 0
G normalize(G) G(1,))
Figure The generating function (GF) semanties [-
16

---

B.1 Generating function semantics
The full generating function semanties can be found in Fig. 5. The first five rules in Fig. 5 were presented in [17,4] The GF for sampling from Binomial( j.P) p) was first given in [26]. Chen et al. [4] also implicitly describe the GF semanties of sampling from Binomial(X,.P)) p) and NegBinomial(X,..
by expressing it as summing X, iid samples from Bernoulli(p) and Geometric(p). respectively.
Let us diseuss each rule in detail. The skip command leaves the distribution and hence the generating function unchanged. Chaining commands Pi: P3 is again achieved by composing the individual transformers for each command. To obtain some intuition for the affine transformation case Xx := a X + C, suppose we only have 2 variables: X1 == 2X + 3X2 + The idea is that the transformed generating function is given by [X1 == 2X1 + 3X2 + (G) = E[u] 2X +3.x+ { (2f Xa] = .ं G(r xfw2).
The semantics of if X E (PI) else } uses GX, If G' is the generating function for mea- sure μ then GxæE is the generating function for the measure yxx, eXLe((S), := nu( {æ E S a E A}). To get intuition for this, remember that for discrete variables X2, we have G G 22) Ex,,x,
X2))) xx xxNxxx22 If we want obtain the generating function for X1 5, we want
to keep only terms where Ti has exponent 5. In other words, we need the Taylor coefficient g(0.)))) This should give some intuition for Gr,E The semanties of if X A (Pi) } else { transforms GxLE in the then-branch and the complement G - GXLE in the else branch, summing the two possibilities.
s!
For sampling, say X1 D, we first marginalize out the old X₃ by substituting for 21 in G: [| = E
Xi x X²] = G(1, Then we multiply with the GF of the new distribution D. The rules for compound distributions D(X;) are more involved, so we offer the following intuition.
For the compound distributions (nn) Binomial(n. (n. p), Neg p), Poisson(Amm) essential use of their property that (f()))) = (((()())())) So if G(z) == B is the GF of X and Y ~ D(X), then Ex Ex -xx Ey-xx))"""'xx""
we make
Ex
Exi= gf((()))())
Ex (ax x (gf( (D(I) = R((x (D)))) = (f(x..f(D(1))))) This not a fully formal argument, but it should give some intuition for the rule.
For D(p) = Bernoulli(p), the intuition is that if G(x) E|ax1 is the GF of X and if Y Bernoulli then (x, v) = E[x - ((1- - X)y* + = (y- - 1)Xx*] =
G(x) (y (rr).
The GF semantics of fail is the zero function because its distribution is zero everywhere. For normalization, note that substituting 1, in the generating function computes the marginal probabilities. similarly to the substitution of in the semantics of sampling. It scales the GF of the subprogram P by the inverse of the normalizing constant obtained by that substitution.
We hope these explanations give some intuition for the GF semantics, but of course, proof is required.
B.2 Correctness proof
For the correctness proof, we need to describe the set of æ E R² where the generating function is defined. For R > 0 and measure fs R3 define it to be the product of open intervals
T(R, 11) == (Q1, x ... x (Qn, R)
where Q = R the i-th component of μ is supported on N (i.e. p(S) = p( E | E NJ) for any measurable S) and Q = 0 otherwise. The case split is necessary because the generating functions are better- and thus have a larger domain, if the measure is supported on N. In fact, it is important that 2; 0 is possible in the N-case because otherwise we would not be able to extract probability masses from the GF.
Theorem B.1 (Correctness). The GF semantics is correct w.E.t. the standard semantics [-]ssta: for any SGCL program P and subprobability distribution μ on Ro0 we have (g)())) = g =()))))) In particular, it correctly computes the GF of the posterior distribution of P as ([
If (f() is defined on T(R,y) for some R > 1 then [P] is defined on T(R', [P]st())) for some R' > 1. In particular, the GFs [P]e((1) and thus ([Pml are defined on T(R, [P]ssa(p) for some R > 1.
17

---

Proof. As usual for such statements, we prove this by induction on the structure of the program. The program can take one of the following forms:
Skip. For the trivial program, we have [P] = G and [P]sta(()) = μs, so the claim is trivial.
Chaining. If P is P1; P2, we know by the inductive hypothesis that ( (((((ilssta and similarly for Ps. Taking this together, we find by the definitions of the semantics:
Pele
[Pa]et (e (g((())))) P [2 ([Plsaa [f; sel (p)).
[ ((ta(())))))))))))
The claim about the domain of [P] follows directly from the inductive hypotheses.
Affine assignments. If P is X2 = a X cc, we have (A)( = m1 E R" X X + +c] E 5}) and thus
gf ([[PP)bc(44))(2))
d(([Pbad(())))
f a
Ne-1 as tan
tnte k+1
.., Xe+1 ... Nx X,)
-k- uk
(0x-12 ai 17-1) (1 Ne..
(IR+ 2 peett (In ) yohn dp(X (X)
an I IK axtl = ((( s.... T IR = (P({{(((()))())
Fnakk
The claim that it defined on some T(R' [P]sst((()) holds if we choose R' =
min( / R, minf... 2VR) R) > because then | < VR( VR) s R for i + k and
1 < (R)R) I VR VRSA R.
Conditionals. If P Xx A(P) else (12), the GF semantics defines (gf())))x which is the same as gf ( by Lemma B.2. Using this fact, the linearity of gf(-), and the induction hypothesis, we find
gf([Pi]saa gf([ ((((P))aa((xx)))))))))))))))))))
gf gf mxLEa))
gf(ux,EA)) (f (gf)))( (gf(H))xNEA)
=
= (P](t(
Also by Lemma B.2, we know that (gf()))) x,eA (and thus (gf(m) - (gf(H)) defined on T(R. μ)), so by the induction hypothesis, both [Pi] (P)]g(((f(())))x. xxe4) and (Pllgr(ff( gf(x)ee)) are defined on some T(R). Hence the same holds for [P]))((f(()))
Sampling. If P is X D, we find:
gf ((((Pasa()))(()) = xx d((P)at(()))
/ (alkk + 11)x / 2x aD(Y)dd(())))))) du(x)
/ (1nk 11) du(x): / dD())
(=(k) + ())
(p))
Regarding the domain of [P] (g))) note that gf(D) is defined all of R for any finite distri- bution and for the Poisson distribution (cf. Table 4a). For D E (Geometric(p)) Neg p)). the GF gf(D) is defined on (- 1-p' I-p), so pick R² = min(R, r-1 > 1. The GF of
18

---

Exponential(A) A) is defined for log < A.i. 2 < exp(A), so we can pick R' := min e²) > 1. The GF of Gamma(o, (9 is defined for log x < B, f.e. x < exp(3). so we can pick R' == min R, e0) > 1. The GF of Uniform(a, b] defined on E. so H² == R works.
Compound distributions. If P Xx ~ D(X,), we use the fact that (f() = (=) for D(n) (Binomial(n. p), NegBinomial(n, p), Poisson(X n.), which can easily be checked by looking at their generating functions (cf. Table 4b).
/ (([P)ssd(4)))
/ dD(K)))
(ack 11)* x xx (D(x))(()dxXX) / (arkk 11) (a)(
/ (ook + () dp(X)
(aok 1,j (f()) dp(Xx (f()) + +,ا x, (f(())))
(f)))((f()))
Regarding the domain of (Pler(ff())), we have to ensure that [ag . (z() < R. For the (Inomial(1,p)) distribution, we choose R' = vTt > such that | || + paaa < vR(11 p + Pv R) S VR((1 - R) + RV)) < R. For the NegBinomial(1) distribution, we choose R' = min( v R -R-P = min( R, 1+ p(vR-1)
VR-PVR
because for any r E (-R', we have
p
S
PVR R
VR
1 (1 p)z,
R)
VR-p
P
VR.
and thus [ TグT-(1-10) p)ar
R. For the Poisson(A) distribution, we choose R' min( / R,
(og(R))
because then bry - exp(A(<x - 1))1 < exxxp(A
log(R) R.
2A
2A
For the Bernoulli distribution D(X,) Bernoulli(X) X,), we reason as follows:
gf (x))(()= =...
(ar +
/ (ar(k>> 1> 11) x))x(11 (1 + X,(** 1)) du(x)
[(alk - 11) x dx(x) (xx - 1) / (= + 11)* dp(x)
gf()) + 1]) + (Ix 1) / x. 1,(2 + 11) x du(x)
(f()) 1]) + 2イ( 1) 0, / (ook + 11) dp(x)
gf()) 1J) + ( - 1) (,ef(n)(elk++ 1]) ()))
Note that the interchange of integral and differentiation is allowed by Lemma B.3. It also implies that (ff is defined for zy E (0, R) by Lemma B.3. Fail. If P is fail, then [P] (1) is the zero measure and [[]))r(g(()) is the zero function, so the claim holds trivially. This GF is clearly defined everywhere.
19

---

Normalization. For normalization, we make use of the linearity of gf(-).
gf (normalize(p)) gf
(Tdu(x))
gf(H) J d)(())
gf(H)
f1 ... 1x (
gf(())
(I,)
normalize(gf(u))
Furthermore, if gf)) is defined on T(R, (t), then so is normalize(gf(1)) on T(R, normalize(j))) =
p).
This finishes the induction and proves the claim.
Lemma B.2. Let A C N be a finite set. Let be measure on IR" such that its k-th component supported N, ie u(S) 1(fx S X NJ) for all measurable S C R22- Define (XLEA(S): = p({x E S Xx E A)). Let == gf().. Then
0,C((((++)) 0J)
(()) = (
o
SEA
Furthermore, G is defined on T(R then is GX,EA on (x,e) I).
Proof. Since UX+EA EiEA MX where we write -X == Mx4e(1). we have:
gf(ux;=i)(x)) = x ddxx=i(X)
/
Xs xxx.... Xk-1
Xk+1
xxx [Xx = o d
.....II--
Tk Fh41
x 1 X1
Nu Xx11
Xn
0'x Xx
du(xx
il
Or.
=x==0/
The reason for the last step is that for the function f(x) x² we have
0 ad
ifl
(0)
T if
Since 0'-i is for > and for - i, we find that f() = = I] il. Above, we used this fact
with = X₁ and Tg. We continue:
1/
/
Xe-
O'r
Ixal
Xe+1
x
du(X)
Onk
X du(x))
a
drt
i!
21=0
0
x (x((
dr.
i!
.. G(zik + 0))
i!
desired. Note that the integral and differentiation operators can be interchanged by Lemma B.3. T
20

---

Lemma B.3. If the integral J xx du(X) ) for a measure on R2 is defined for all a E ((R,H), we have
Oi ... Dam x du(X)= / Or, .h..0. Dr xx dy(X)
and both sides are defined for all E M.)
Furthermore, the right-hand side has the form xxxxxxxxxxx x du(X) for a polynomial p in n vari- ables and w E No. with the property that pl (X) = O whenever X; < wi for some j. in which case we define p( X)
xxx xxxx
= O, even if x is undefined.
Proof. The proof is by induction on m. If m = then the statement is trivial. For the induction step 0, may assume that
Oig 0....... 0:m / x xx d)(x))
= / 0.2..... D. x du(x)) = 1--pp((X dj((X)
By Lemma B.4, we reason
0.,....... -mm / xx du(X) Or1 / p(x)x=* x- du(x)
/ (Xis wi, )e) X )ax x 20 ni +11 d)())(()
/ 0r, ... ... n. d)(x)
which establishes the induction goal with the polynomial p( X) == (Xi, - wi, (x() and i ==
wlin + Wi, + +1]. The property p p(X)= O whenever X3 < w for some j still holds due to the added factor (X - we,) taking care of the case X,, = - Wi 1.
Lemma B.4. Let w E N", be a polynomial in n variables, and p. a measure on R We define p(X)x* x
0 whenever p(X) 0, even if æ S 0 and X2 < wy for some j. Suppose pl p(X)x is defined y-almost everywhere and j-integrable for all æ E T(R, u). Then
0z: / p(x)xx du(x)-
/
Oz. p(x) x du(X)
/ (Xi w/ w))p(XXet x du(X)
holds and is defined for all p.
Proof. The proof is about verifying the conditions of the Leibniz integral rule. The nontrivial condition is the boundedness of the derivative by an integrable function. The proof splits into two parts, depending on whether the i-th component of μu is supported on N. If it is, then x; E (-R, Rt) by the definition of T(R. μ1) and can be nonpositive, but X; w. and p(X) A happens plalmost never. Otherwise, we have negative. Hence the two cases need slightly different treatment.
(0, R) and it is thus guaranteed to be positive, but X; w, may be
We first deal with the case that the i-th component not supported on N. Let 0 < R' < R < R.
We first prove ø2z Leibniz integral rule. For this purpose, we need bound the derivative:
a. f p(X)x* (x(X)
se, ((X) X) ) dp(() for any H; E (e, I2), using the
i dz
< 1 _ wxIl((X)axx- x;
pl
w
$ ((X, w; w))pp( X x wib
1 X; wi [pllp((X x
If X; > wi, we can choose M > 1 sufficiently large such that [X, - wgln X. < (X; wy) Rt < - whenever X, - wA > M. (Such an M can found because the exponen-
tial function y +
Rd
grows more quickly than the linear function y + y.) If X, - w, S M
21

---

then we also have IX; wi | x, x S .MM.. R"x:-";, So the derivative is thus bounded on the set {X X; > w;} by M. 1p X (a + R*]) x- which is integrable by assumption.
On the complement set {X 0 X w,), we find X, wy wy|x 1, x.=w: < |X, w. x
wie Hence the derivative is bounded by |p(X ](x e)) e)) which is again integrable by assumption. So the derivative is bounded by an integrable function and thus, by the Leibniz integral rule, interchanging differentiation and the integral is valid for all E I²) and thus for all I; E (0, R).
Next, consider the case that the i-th component of μ is supported on N. Let 0 < R' < R" < R. Note that the set {X X, < w, A p(X) o 0} has measure zero because otherwise the integrand p( X J)x - would not be defined for Ia < O due to the negative exponent X; wy. Hence, for X 2 w p( x) and can bound the partial derivative for | < Rt:
j--ammstt all X, we have p(X)=
p( (x) x)a= x x-
S sI(x; wx:p(( X |ac x w +1]
Or:
$ [X 2 w||XX - w/ [|pppp(()]))))) ||2]] x X- -ww++u++++
x X 2 w + + 1]] X; wi X X -- w[i-h
X wii + w. + 1] X; w.llp(( | Rd -
X Z w||Xx - ws wl X [|a R1 x
because the factor [X, w/| vanishes for X, wy and thus ensures that negative values for the exponent X - wwli + w, + 1] don't matter, that we can use the monotonicity of | X wiw in lz].
Similarly to the first case, for a sufficiently large M > max(1, w:). we have |X; w; W X'' -o < R"X,-";, whenever X; w; > M. If X; w; < M then [X 2 w||X; wal I M and we also have X 2 wi X; w/ W.'R'' < M R" The derivative is thus bounded by M |p(X) (()li + R" R"IXX- >- which integrable by assumption. By the Leibniz integral rule, interchanging differentiation and the integral is valid for all Ia E R², R²') and thus for all
r, E (-R, R).
L
Remark B.1. As an example of what can go wrong with derivatives zero if the measure is not supported on N. consider the Dirac measure μ = ()rac()). It has the generating function G(x) == (f)((()) vz defined for Rzo- Its derivative is G'(x)= 2/0 which is not defined for = 0.
B.3
Possible extensions to the probabilistic programming language
The syntax of our language guarantees that the generating function of the distribution of any proba- bilistic program admits a closed form. But are the syntactic restrictions necessary or are there other useful programming constructs that can be supported? We believe the following constructs preserve closed forms of the generating function and could thus be supported:
Additional distributions with constant parameters: any such distribution could be supported as long as its GF is defined on [0, 1 + e) for some > 0 and its derivatives can be evaluated.
Additional compound distributions: we spent a lot of time trying to find additional compound dis- tributions with closed-form GF since this would be the most useful for probabilistic models, but we were largely unsuccessful. The only such distributions we found are Binomial(m., X&),, which could just be written as sum of iid Bernoulli Xg)-variables and Gamma(Xx. .9) with shape parameter X3., which could be translated to GF with the same idea as for Binomial(X,, p). NegBinomial(X,, p), and Poisson(X Xx))
Modulo, subtraction, iid sums: The event Xx mod 2 = 0, and the statements X₁ max (X,- m, 0) (a form of subtraction that ensures nonnegativity). and Xx == sum_ (i__idd) X1) (i.e. X₃ is assigned the sum of X, independent random variables with distribution D) can also be supported as shown in [4, 18].
Nonlinear functions on variables with finite support: if a variable has finite support, arbitrary functions on it could be performed by exhaustively testing all values in its domain.
22

---

Soft conditioning: another supportable construct could be scorea x* for E [0, 11, which multiplies the likelihood of the current path by aNx,, We think score (Xx;) where q(t) == m Li=o 9; is polynomial with coefficients di E [0, 1] and with S on the domain of X could also be supportable using similar techniques as for observations from Bernoulli(X,,
None of these extensions scemed very useful for real-world models, in particular, they are not needed for our benchmarks. For this reason we did not formalize, implement, or prove them.
Support for loops is another big question. Bounded loops are easy enough to deal with by fully unrolling them. However exact inference for programs with unbounded loops seems extremely difficult, unless the inference algorithm is given additional information, e.g. probabilistic loop invariant as in [4] or at least a loop invariant template with a few holes to be filled with real constants [18]. However, finding such a loop invariant for nontrivial programs seems exceedingly difficult. Furthermore, the variables of looping programs may have infinite expected values, so the generating function may not be definable at the point 1. But in our semantics, evaluating is very important for marginalization. For these reasons, we consider exact inference in the presence of loops an interesting but very hard research problem.
C Details on implementation and optimizations
Overview Our tool takes SGCL program P with n variables input, and a variable X; whose posterior distribution is to be computed. It translates it to the GF [P]lar(1) according to the GF semantics and normalizes it: G := (rrmalize[[[]]ee (1)). This G' is the GF of the posterior distribution of the program, from which it extracts the posterior moments of X, and its probability masses (if X, is discrete).
Computation of moments First, it computes the first four (raw) posterior moments from the factorial moments, which are obtained by differentiation of G:
M1 E[X,] a,G(1,)
M2 E[X?] [X(X,)(X, - =f
M3 == EXX;]
+ - 1) [(X.)
M4 = E[X,|| = EXX,(X,, - 1)(, -04-2)(44 -X),x, - 3) 6E - 1)(X,
+7E[XX(X;- - - E[X]
=0+G(1,)) + 60, +
From the raw moments, it computes the first four (centered/standdrddzze)) moments, i.e. the expected value (n = E(X.). the variance (o² == V(X.)), the skewness (E[(X; / and the kurtosis
(X, /) )
= [[X = Mj
²² = v[X.] E[(X; 11)
E[X?| [x + 11² M2 112
B[(X; -4)31
|| 3x 3,2 -
Skew[X,]
x3
σ3
M3 3uM2 -21
g3
B((X; - 11)41 0¹
B|X,] - 41E[X 6nPE[X]] 41 + 4
Kurt[X.]
M4 411M3 + .4. M2 3/d¹
Computation of probability masses If X; is a discrete variable, the tool computes all the prob- ability masses P[X; = m) until the value of m such that the tail probabilities are guaranteed to be below certain threshold, which is set to P[X, 2 m] < 256 This is achieved by setting
23

---

m - / (X 2u4 where : B[X,]. Then we find
P 2 m| E | A 2 B = =
14)² -4) 256E[ 1 S 256
by Markov's inequality. Note that in practice, the tail probabilities are typically much smaller than usually in the order of 10-5, For all & E (O,....m),, the tool computes the posterior probability ki by marginalizing out all the other variables (substituting for them) and computing the Taylor coefficient at 2; 0:
F[X
1 k! 0r -(
P(X: kj
would be desirable to compute posterior densities for continuous distributions as well. In fact, there are mathematical ways of recovering the probability density function from a generating function via an inverse Laplace transform. However, this cannot automated in practice because requires solving integrals, which is intractable.
Implementation details Our tool Genfer implemented in Rust [20], a safe systems programming language. The main reasons were low-level control and performance: the operations on the Taylor polynomials to evaluate derivatives of the generating function need to be fast and are optimized to exploit the structure of GFfs arising from probabilistic programs (see Appendix C.3). C or C++ would have satisfied the performance criterion as well, but Rust's language features like memory safety, enums (tagged unions), and pattern matching made the implementation lot more pleasant, robust, and easier to maintain. The first author's 's experience with Rust was another contributing factor.
The coefficients of the Taylor polynomials are stored in a multidimensional array provided by the poating-point numbers and unbounded rational numbers are provided by the rug library", which is an interface to the GNU libraries GMP (for rationals) and MPFR (for floats). Our implementation is available on GitHub: github.com//raiier/gentte
ndarray library*. Arbitrary-precision
C.1 Automatic differentiation and Taylor polynomials
The main difficulty in implementing the GF semantics is the computation of the (partial) derivatives. This seems to be a unique situation where we want to compute d-th derivatives where d is in the order of hundreds. The reason is that the observe X₁; construct is translated into d-th (partial) derivative and if d is real-world observation, it can be large. We have not come across another application of automatic differentiation that required derivatives of total order of more than 10.
As a consequence, when we tried PyTorch, a mature machine learning library implementing automatic differentiation, the performance for derivatives of high order very poor. Therefore, decided to implement our own automatic differentiation framework. Our approach is to compute the Taylor expansion in all variables up some order of the generating function instead of symbolic expres- sion of the d-th derivative. This has the advantage of growing polynomially in d, not exponentially. Contrary [27], is advantageous use Taylor coefficients instead of the partial derivatives because the additional factorial factors are can easily lead to overflows.
Taylor polynomials More formally, we define the Taylor polynomial Taylor (G) of a function G: 1. + R at R of order the polynomial
Taylor ao (G) =
Σ
go
a
where we used multi-index notation: |ox| stands for on +... + On ae! for an!.... an, and d for 00'
...0 "..... The coefficients ca O G(w) are the Taylor coefficients of G at and are stored in multidimensional array in our implementation.
io/ 'tttps:/crrates co/crate//rrg
24

---

Operations The operations on generating functions (Appendix B) are implemented in terms of their effect on the Taylor expansions. Let G. H. I= R be two functions with Taylor expansions Taylor ur (G) = Ejaled Ia (z w) and Taylor to (H) Ejarga (a w) Addition F
G+ H of two generating functions is implemented by adding the coefficients: Taylor (G H) = La + ha)(2 w)a Scalar multiplication F' c G is implemented by multiplying the coefficients: Taylor (c-G) = Ejalsa( 9a (xx -w) Multiplication F = G H of two generating
functions is implemented by the Cauchy product:
Taylor wo (G H) Σ
Σ aa--aa-a
Dar Mag
(x -x)
Sd (X)
Division, exponentiation, logarithms, and powers can be implemented as well (see [13, Chapter 13] for details). Partial derivatives essentially correspond to shifting the index and a multiplication:
Taylor w ((0,)
d-1 (Ox.))=
Σ (as + 1)a
|elsd-1
The most complicated case composition/substiuttion, i.e. F(xx) = G(æ + H(x))). For this, ha ha(( w)s and we let Taylor (G) Ejlsa w²) Then we have (where "h.o.t." stands for "higher order terms")
we let w' : wik + ho] where Taylor (H) = Eiassa
G(æ(k + H(x))) = G x: k Elalca w)) + hoo.))
Σ Ja
x k ho Eisisisa ha(x w))
8 + thott ) a + h.o.t.
aa| Sd
Σ Ja (ax w) k+ Eislelsa (ale-ww)s))) ))))) a + h.o.t.
la Sd
because by definition x + hol -ol-wi w' (x w) [t + 0). This means that we can obtain Taylor (F) by substituting Σ Eiglaia (a w for ox w/ in Taylor ( (G).
Example C.1 (Calculations with Taylor polynomials). Here we show how the expectation for Example 2.1 can be computed using Taylor polynomials. We use the same naming as in Example 3.1.
Recall that to compute E[X], we need to evaluate 0, E(1, 1), so we compute the Taylor expansion of E at (1, 1) of order 1, i.c. Taylor
To compute E(z, y) = D(1,I)' we need D(x, y), so
Diz.y)
we need Taylor[1.1) (D). For D(x. y) = 2 O0 0), we need the Taylor expansion of C at (1, 0) of order 1 + 2 = 3. For Czz,v) v) (0x + 0.9). 1). we need the Taylor expansion of B at (1 0.9), 1) = (0.9. 1) of order 3. Since Bzx, y) exp(20(x 1[) = exp(20((x 0.9)
0.1)) = exp(20 - 0.9) - 2), this Taylor expansion is
4000
Taylor (B) = -2 2 20e (I 0.9) 200c (x 0.9)²
2 (r 0.9)
3
Since C(1 + (x: 0),0+ y) B((1 + (.: 1y000.1y+009)) B(0.9 + 0.9(z 1) + 0.1y 0. - 11y, 1), the Taylor expansion of Ct (1,0) is given by replacing (æ - 0.9) with 0.9(x - 1) 0.1y + 0.1 - 1)y and we get
Taylor (1.0) (C)
(C)) = e --2(1 2y 2y² 3 q/² + 18(x - 1) + +8(xz - 1)y + 40(x 1)y²
+ 162(x 1)2 360 1)2y + -x)3)
For = =00 0), we compute the second partial derivative wrt. y, substitute y = 0 and then multiply by y² = (1+ +++17 (y - yielding
1 3/00 (4 + 8y + 80(x 1) + higher order terms)l
D(z v)
(4 + 80(z 1) + higher order terms)(1 + (y - 1))²
- + 80(x - 1) + 8(y - 1) + higher order terms)
As consequence, we find that D11, 1) 4e -2, so Taylor = 1 20(x 1) + 2(y 1). So
[[X] (1, 20, as desired.
25

---

Memoizing intermediate results When computing the Taylor expansion of function G, it is important memoize the intermediate Taylor expansions of subexpressions of G they oceur more than once. For example, conditionals are translated to generating function that uses the previous generating function twice. Evaluating it repeatedly would double the running time for each conditional. With memoization, it is possible to handle programs with lots of branching (e.g. > 2000 paths in the mixture model, cf. Table 5) in a reasonable amount of time.
C.2 Optimizing Observations from Compound Distributions
Observations from compound distributions are generally the bottleneck for the running time because the construct observe d ~ D D(Xx) ) syntactic sugar for X,+1 D(xx): observe X = which introduces new variable, which is immediately discarded after the observation. This expansion has the semantics
[observe d (G)(() [[observe X = (] D(Xx)) Tn+1)
Od
d' Ox [xat1
In+11
containing an extra variable In+1, which worsens the running time of the Taylor approach signifi- cantly. We can find more efficient ways of computing this semantics.
Theorem C.1. Observing from compound binomial distributions can be implemented without intro- ducing new variable:
1 [observe ~ Binomial(X p) (nrr))" ( (1 p)ax))
Proof. A proof for the discrete setting was given in [27]. In our setting, which allows continuous distributions, we cannot use the usual argument about power series anymore. Instead, we reason as follows:
Jobserved d ~ Binomial(Xe:P))lf((
Od
[Xn+1 Binomial(Xa)ls xn
d! 0un+1 gd G(ælk ay (pen+1 + 1- p)])
In-1=0
d! On
I=+1=0
GaG(ask + ox (pxn+1 + p)l) (porx)"
d!
xm41=0
d! (pmx) ( Cg + (1 p))aa))
where follows by iterating the following argument:
G(ælk T& (pæn+1 + 11 (11pp]) p)])
0x
= axG + 2k (pen+1 + 1 - --]) O(o (pxn+1 + 1 p))
axnt1
OgG(allk + 2b (prn+1 + 11 p)]) prr
which holds by the chain rule.
Theorem C.2. Observing from compound Poisson distributions be implemented without intro- ducing a new variable:
[observe d /oape
where
(G)(())(z)) == Ax G())
26

---

Proof. Let F R² + R be smooth function. Then
F( + In e ((en+111)))
0znt1
O,F(alk + ay e 1)) 12 X
(xx(()) oy
Inductively, get
Glælkk + - (( D (z)k (C) - (
Oxd
Substituting In+1 H 0 in the final expression and dividing by d! yields:
[obseeve d Poisson (riison(AXx)lloe (IXx)]le(G))
[Xn+1 ~ Poisson Poisson (X))]er #n41)) Od
Intz)
d! dad
#o41=U
d G(zlk + a& (
d! du
Int1=0
Dd.
Theorem C.3. Observing from compound negalive binomial distributions can be implemented without introducing new variable:
d ofF(kkk p. p'z (1 - p)d La.i
[observe d ~ NegBinomial(Xx. Plls
Σ
where the numbers La, are known Lah numbers and defined recursively follows:
Lo.0 ==
Laji = 0
for i < or <i for S sisd+1
La+1, = (d + i)L + La.-1
Proof. Let R + be smooth function and y R + R given by y(z) = 1-(1-p)= Then
(xxx)
p)) p)v(x)² and
(1-(1-p)) p)a)
F&
F F(x + 1 - (1 p)
0zn+1
O (xlk + T&P ((2n+ 1)) p.s. (1 - p1))((2n44))
We claim that write
dd
p
F(alk 1 (1 p)
Dxd n+
F( + v(2n+1))l p'z (1 p)d Las v(zn+1)
0
The claim is proved by induction. First the case d 0:
F(zlk + 1- (1 p)rnt1
pæz
= F(ze(h: + 1 (1 pa
(1 p)" 2 Lo.0 u(In+1)"
D)
Next, the induction step:
8d+1
F(x[k +
1 (1 p
0z 2x
x
27

---

par Σ OkFlolk + par v(xn+1))) p'af. (1 - L) y(anti)dtt
0ont1
oft F(xz(k + per v((n+1))) pazy (ma+1) p1a1(1 p)_II Ld (0
1=0
0g F(x pa v((n+1))) p/ar p)² La ()y(((nt1) 1) y (vn+1)
Σ 0 F(x + par v(xn+1))) pl sit (1 - (1-pp)+(((((n)))) y(un++)² 1a
0. F( - par v(2n+1))) p'zk (1 - 1ltpp)=+((yxxxpppp y(En+1)² (d (2)1)((((2)))) ddti-1
1=0
o Fax pak v(xn+1))) p'at (1 Ledi 1 d+i+1
d+1
oe F(x pa v(=n+1))) plat (1 p) + Ta (x dditt1
+
4+1
Σ
(ark + peu s((n+1))) piaj (1 p)a+1((())))))))) ((d + id + 1da-1))((nt1)) dtiti )
1=0
F(x(k + paz v(zn+1))) plat (1 p) La (x-1,)((2111 dtit1
Substituting Fn+1 + 0 in the final expression and dividing by d! yields:
[observe d ~ (igBinommal(xxxP)eeee
1 gd
[Xa+1 N
knt1)
d! dr
lan mm
gd
Gl 1- (1 1p)xx+t lmm
2k p
d! Ord n+1
+ plack (1 (1-p) L)
Theorem C.4. Observing from compound Bernoulli distributions can be implemented without introducing new variable:
G(z) Gz0) 0,G(x)
d=0 d d otherwise
[observe ~ Bernoulli(Xx)] =
0
Proof. We argue as follows:
[observe d ~ Bernoull((Xx)]le((G]]
3d
[X [en
d! Ord
#a11=0
gd
(G((z) + - 1) - G(x))
d! dr.
28

---

If d 0. the right-hand side is G(z) G,)) If 1, the right-hand side (x), For
larger d 2 2, the dl-th derivative vanishes, which completes the proof.
C.3 Analysis of the running time
If the program makes observations (or compares variables with constants) di,..., .,.....,,,. .......,mm.... the running time of the program will depend on these constants, Indeed, each observation of (or comparison of a variable with) d, requires the computation of the d,-th partial derivative of the generating function. In our Taylor polynomial approach, this means that the highest order of Taylor polynomials computed will be d ;= > ii= dy. If the program has n variables, storing the coefficients of this polynomial alone requires (dd)) space. This can be reduced if the generating function is a low-degree polynomial, in which case we do not store the zero coefficients of the terms of higher order.
Naively, the time complexity of multiplying two such Taylor polynomials is O(d2)) and that of composing them is O(d"). However, can do better by exploiting the fact that the polynomials are sparse. In the generating function semantics, we never multiply two polynomials with n variables. One of them will have at most two variables. This brings the running time for multiplications down to (d Furthermore, the GF semantics has the property that when composing two generating functions, we can substitute variables one by one instead of simultancously, bringing the running time down to (n d.n+1). Even more, at most one variable is substituted and the substituted terms have at most two variables, which brings the running time down to (d
Theorem C.5. Our exact inference method can evaluate the k-th derivative of a generating function in O(s dr and O(sd)) ) for n = 1, where s is the number of statements, d (for data) is k plus the sum of all values that are observed compared against in the program, and n. is the number of variables of the program.
Proof. The above arguments already analyze the bottlenecks of the algorithm. But for completeness' sake, we consider all operations being performed on generating functions and their Taylor expansions (cf. Appendix C.1).
The maximum order of differentiation needed is d and since there are variables, the maximum num- ber of Taylor coefficients that need to be stored is Odd"). Addition of two such Taylor polynomials and scalar multiplications trivially take Odd)) time. Since in the GF semantics, a GF is only ever multiplied by one with at most two variables, the multiplications take 2 time. Division only happens with a divisor of degree at most 1, which can also be implemented in O(d"+1) time [13, Chapter 13]. Exponentiation exp(...) ), logarithms log(...). and powers (...) ym are only performed on polynomials with one variable, where they can be implemented in ((d)) time [13, Table 13.2).
Finally, the only substitutions that are required are ones where a term with at most two variables is substituted for a variable, i.e. of the form p(z=lk: + g(21, x2))) where p(x:) = Pa. Then
pk + g(21, x2)l) (o( + q(11, x2) (pi(x) + g(x1, (p2(x) + +...))
where (x) = E Easa= aak=00) This needs d multiplications of g(#1, in 2 variables and polynomial in variables, each of which takes +d)) time. In total, the substitution can be performed in 33) time. Finally, if there is only one variable, the substitution can be performed in O(d)) time because q can only have one variable.
In summary, the input SGCL program has s statements, and each of them is translated to a bounded number of operations on generating functions. These operations are performed using Taylor expan- sions of degree at most d, each of which takes at most 0(d++3) time (and at most (d²) time for 1). So overall, the running time is O(s d3+3) (and O(sd3) for n n==1)
Potential improvements of the running time Using various tricks, such as FFT, the running time of the composition of Taylor polynomials, i.c substitution of polynomials, can be improved to O(dlog 3n/22 [2]. This is still exponential in the number of variables n, but would improve our asymptotic running time for n < 6 program variables. tion of multiplication and composition based on the fast Fourier transform (FFT), but discarded this option because it led to much greater numerical instabilities. is possible that the running time could instead be improved with the first algorithm in [2], which relies on fast matrix multiplication instead
fact, we experimented with an implementa-
29

---

of FFT. In addition, it should often be possible to exploit independence structure in programs, so the GF can be factored into independent GFs for subsets of the variables. This not implemented yet because it was not relevant for our benchmarks.
Performance in practice Note that, for many programs, the running time will be better than the worst case. Indeed, substitutions are the most expensive operation, if the substitute expression is of high degree. For some models, like the population model, the substitutions take the form of Theorem C.I. i.e. Ga (1 p)ar.)) which is of degree 1. Such a substitution can be performed in (dd) time, because one simply needs to multiply each coefficient by a power of (1 - p)) This explains why the population model is so fast in the experimental section (Section 5). In other cases as well, we can often make use of the fact that the polynomials have low degree to reduce the space needed to store the coefficients and to speed up the operations.
C.4 Numerical issues
If our implementation is used in floating-poinn point mode, there is the possibility of rounding errors and numerical instabilities. which could invalidate the results. To enable users to detect this, we implemented an option our tool to use interval arithmetic in the computation: instead of a single (rounded) number, it keeps track of lower and upper bound on the correctly rounded result. If this interval is very wide, it indicates numerical instabilities.
Initially, we observed large numerical errors (caused by catastrophic cancellation) for probabilistic programs involving continuous distributions. The reason for this turned out to be the term log(x) that occurs in the GF of the continuous distributions. The Taylor coefficients of log(z) at points 2 1 are not bell-behaved::
1-1)"2-1
log(x) log(z)) +
1
because the coefficients grow exponentially for (due the term 2 ") and the alternating sign exacerbates this problem by increasing the likelihood of catastrophic cancellations.
To fix this problem, adopted slightly modified representation. Instead of computing the Taylor coefficients of the GF G directly, we compute those of the function H(x) == G(x') where xy == x; X; is discrete and x² == exp(:;;
X; is continuous, thus canceling the logarithm." 6 (The reason we don't use x' == exp(x) for all variables is that for discrete variables, we may need evaluate at x = 0, corresponding to z; = -o, which is impossible.) It is straightforward to adapt the GF semantics to this modified representation. Due to the technical nature of the adjustment (case analysis on whether an index corresponds to a discrete variable or not). we do not describe it in detail. With this modified representation, avoid the catastrophic cancellations and obtain numerically stable results for programs using continuous distributions as well.
D Details on the Empirical Evaluation
a laptop computer with Intel@ TM 15-8250U CPU @ 1.60GHz processor and 16 GiB of RAM, running Ubuntu 22.04.2 and took few hours to complete overall. The code to run the benchmarks and instructions on how to reproduce our experiments are available in the supplementary material.
All experiments
D.1 Comparison with exact inference
We manually patched the tools Dice, Prodigy, and PST measure and output the time taken exclusively for inference to ensure a fairer comparison and because the time tool takes to startup and parsing times are not very interesting from research perspective. We each tool 5 times in and took the minimum running time to reduce the effect of noise (garbage collection, background activity on
Note that this boils down to using the moment generating function (MGF) (x,))(xx.X.) for continuous variables X, and the probability generating function (PGF) B[a;"') for discrete variables X.. In mathematics, the MGF is more often used for continuous distributions and the PGF more commonly for discrete ones. It is interesting see that from numerical perspective, this preference is confirmed.
30

---

Table Summary of the benchmarks: n is the number of variables, is the number of observations, d is the sum of all the observed values, ie. the total order of derivatives that needs to be computed, the number of statements, and p the number of program paths. Benchmark
d
continuous priors?
population model (Fig. la)
4
254
13
no
population (modified, Fig. 1b)
4
254
21
16
population (two types, Fig. lc)
2
8
277
30
n0
switchpoint (Fig. 2a)
109
188
12433
111
yes
mixture model (Figs. 2b and 2c)
2
218
188
329
2
6
10
no
HMM (Fig. 2d)
60
152
2
10
n1
the computer). In rational mode, we ran Genfer with the option --rational and Dice with the option -mm--tyee For the "clinicalTrial" benchmark, Genfer required 400-bit floating-point numbers via --rrecision 400 because 64-bit numbers did not provide enough accuracy. For better performance, PSI was run with the flag dd on finite discrete benchmarks.
The benchmarks were taken from the PSI paper [9, Table 1], excluding "HIV"" "LinearRegressionn" "TrueSkill", "AddFun/max", "AddFun/summ,, and "LearningGaussin"" While Genfer supports continu- ous priors, the excluded benchmarks use continuous distributions other ways that are not supported (e.g. in observations). Note that Dice and Prodigy have no support for continuous distributions at all. These benchmarks were translated mostly manually to the tools' respective formats, and we checked that the outputs of the tools are consistent with each other.
D.2 Comparison with approximate inference
Experimental setup As explained in Section 5, we ran several of Anglican's inference algorithms. Each algorithm was with sampling budget of 1000 and 10000. The precise settings used are the following:
Importance Sampling (IS): has no settings
Lightweight Metropolis Hastings (LMH): has no settings
Handom-Walle Metropolis Hastings (RMH): has two settings: the probability a of using local
MCMC move and the spread of the local move. We used the default = 0.5, 1 and another setting facilitating local moves with a 0.8, o =
Sequential Monte Carlo (SMC): number of particles E {1, 1000}. The default is and also picked 1000 to see the effect of larger number of particles while not increasing the running time much.
Particle Gibbs (PGibbs): number of particles {2, 1000}. The default is 2 and we also picked 1000 as the other setting for the same reason as for SMC.
Interacting Particle MCMC (IPMCMC): it has the settings "number of particles per sweep", "number of nodes running SMC and CSMC", "number of nodes running CSMC", and "whether to return all particles (or just one)". We left the last three settings their defaults (32 nodes, 16 nodes running CSMC, return all particles). For the number of nodes and the number of nodes we chose the default setting (2) and a higher one (1000) for the same reason as SMC.
full list of Anglican's 's inference algorithms and their settings can be found at https: / /pppbprrog.o..ooooooo 20 times each benchmark
gi thub. io/ /aagnnll i can/inference, Each inference algorithm and the running times were averaged to reduce noise. This took around an hour per benchmark.
D.3 Benchmarks with infinite support
A summary with important information on each benchmark can be found in Table 5, including some statistics about the probabilistic programs. We describe each benchmark more detail the following.
Population models The probabilistic program for the original population model from (26) shown in Fig. 6a. This program also be written without the extra variable New by writing the two
31

---

0, Ni Poisson (Xg) (10
N Poisson(Au);
N2 Poisson(A A22));
N ~ 3inomial(N,o)
N2 +~ Binomial(N1,17)
Disaster ~ Bernoulli(0.1);
N1 (nnomial(N1,61))
N: Poisson(Au);
if Disaster = 1
N2 Binomial(N2,.22)
N +~ Poisson(Xi)
N1 +~ Poisson (X)):
N ~ Binomial 6);
} else
N2 +~ Poisson
New ~ Poisson(A1);
N +~ Poisson(A1)
observe (1) Binomial(AV,,
N == New + N;
observe y1 Binomial(N,
observe y1 Binomial pN;
observe vi (2) Binomial p);
(a) Original model from [26]. (b) Randomly modified arrival rate. (c) Two interacting populations.
Figure The program code for the population model variations.
statements involving it as N + (Ai) Here we used the same parameter values as [26]: 5 0.2636, À = A (0.0257, 0.1163, 0.2104, 0. 0504, 0.0428). and A 2000 is the population size parameter. Note that the largest population size considered by [26]
= 500. For the observation
rate p. Winner and Sheldon [26] consider values from 0.05 to 0.95. For simplicity, we set it to p = 0.2. The y/ are k = 4 simulated data points: (45, 98, 73, 38) for = 2000.
The modified population example is shown in Fig. 6b, where the arrival rate is reduced to X := 0.1 in the case of a natural disaster, which happens with probability 0.1. The rest of the model is the same.
The model of two interacting populations (multitype branching process) is programmed shown in Fig. 6c, where ( 0.9A, (2) 0.1A, == 0.1. In this model, some individuals of the first kind (i.c. in N1) can turn into individuals of the second kind and are added to N2. The other parameters of this model are the same as before. The data points are (1)) (35, 83, 78, 58) and y (3, 6, 10, 4).
Bayesian switchpoint analysis Bayesian switchpoint analysis is about detecting change in the frequency of certain events over time. An example is the frequency of coal mining disasters in the United States from 1851 1962, as discussed in the PyMC3 tutorial [23]. We use the same model as in [23], and its probabilistic program is shown in Fig. 7a. This situation can be modeled with two Poisson distributions with parameters A₁ (before the change) and A2 (after the change). Suppose we have observations ,.,.... In. with 3/+ Poisson(A.) if < T and y/: ~ Poisson A2) if 2 T. The parameters A1, A2 are given exponential priors Exponential(1), The change point T itself is assumed The probabilistic program code is shown in Fig. 7a, to be uniformly distributed: T ~ Uniform((..n)
where n = 111 and y = (4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6. 3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 0. 1,
2, 2, 2, 3, 4, 2, 1, 3, n/a, 2, 1, 1, 1, 0,
2, 0. 1. 1, I. 0, 1, 0. 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0,
2, 1, 0, 0, 0, 1, 1, 0, 2, 3, 1, n/a, 2, 1, 1,
0.
0, 1,
1) where "n/a" indicates missing data (which is omitted in the program).
In fact, this program can be rewritten to the more efficient form shown in Fig. 7b, to minimize the number of variables needed.
Mixture model In the binary mixture model, the frequency events are observed from an equal- weight mixture of two Poisson distributions with different parameters A,, , A2, which are in discrete steps of 0. This is implemented by imposing geometric prior and multiplying the rate by 0.1 in the Poisson distribution: Poisson(0.1 A1). The data the same as for the switchpoint model. The task is to infer the first rate A₁. The probabilistic program is shown in Fig. 8a.
32

---

T' ~ Uniform(1..
A ~ Exponential
if 1 Bernoulli(1/m)) {
observey: Poisson(A);
T Uniform(1..7
A Exponential(1);
Exponential());
observe //2 Poisson(A);
A Exponential(1):
if 1 < <
observe Poisson(A1);
observe yn Poisson(A);
else {
T: 1;
observe yn Poisson(Az):
else (if Bernouli((1/(1 1)) {
observe y1 Poisson(A):
observe 3/2 ~ Poisson(A);
A ~ Exponential(1);
if n T
observe Poisson(A) A);
observe Yn ~ (oi)s
else {
observe yn Poisson(A2):
observe yr. Poisson(A):
2;
else {
aa) Switchpoint model from [23].
}
(b) More efficient version of the switchpoint model.
Figure 7: The program code for the Bayesian switchpoint analysis.
Z = 1;
A1 Geometric(O.1))
A Geometric
Ag Geometric(0.1);
A2 Geometric
Z 0
if 1 ~ Bernoull((0.5)
observe ~ Poisson(O.1 Ar);
observey) Poisson (0.1 A1)
Z Bernoull
} else {
else
observe Poisson(O.1 A2)
observe y1 Poisson A2);
Bernoull(00.8;
if 1 Bernoulli(0.5)
7 = 0
observe you Poisson(O.11 A1)
observe Um ~ Poisson(O.1 -A))
} else {
Z Bernoulli(0.2);
observe Poisson(0.1 A2)
else
observe Vm
Poisson A2):
Bernoulli(0.8);
(a) Program code for the mixture model.
(b) Program code for the HMM.
Figure 8: The program code for the mixture and HMM model.
33

---

0.7
0.5
0.3
0.15
3.03
0.C
3.0
0.
25
1.0
2.5
20
1.3
1.5
Timne
(a) Error the mean
Error the standard de- (c) Error of the skewness viation
(d) Error the kurtosis
Figure 9: Comparison of moments for the original population model.
1.3
2,03
Nokes
20
Time
(a) Error the mean
(b) Error the standard de- (c) Error of the skewness viation
(d) Error the kurtosis
Figure 10: Comparison of moments for the modified population model.
Hidden Markov model The hidden Markov example based on [22, Section 2.2]. The hidden
state % can be and transitions the other state with probability 0.2. The rate of the Poisson- distributed number of observed events depends on the state Z and is 0.1 A1 or 0.1 A2, where we impose a geometric prior on A1 and Ag. As for the mixture model, this discretizes the rates in steps of 0.1. The inference problem is to infer the first rate A₁ The probabilistic program is shown in Fig. 8b, where Y/ are m 30 simulated data points from the true rates \ == 0.5 and X2 == 2.5. Concretely, we have (2, 2. 4, O, O, O, O, O, 1, 1. 0, 2, 4, 3, 3, 5, 1, 2, 3, 1, 3, 3, 0, 0, 0. 2, 6, 1).
D.4 Additional results
In the main text, the data presented mainly in form of the total variation distance (TVD) between the true posterior (computed by our method) and the approximated posterior (computed by Anglican's inference algorithms). In this section, we present additional evidence in two forms:
errors of the posterior moments (mean, standard deviation, skewness, and kurtosis) in Figs. to 14, and
histograms (Fig. 15) of the sampled distribution produced by the best MCMC algorithm, where the sampling budget is picked such that the running time is close to that of our method (e.g. 1000 for the (fast) population model and 10000 for the (slower) mixture model).
The additional plots confirm the conclusions from Section 5.
3.5 BMCNC
0.4
033
7080s SW
PORSE
02
3.3
022
3
(a) Error the mean
(b) Error the standard de- (c) Error of the skewness
(d) Error the kurtosis
viation
Figure 11: Comparison of moments for the two-type population model.
34

---

2.0
14
(a) Error the mean
(b) Error the standard de- (c) Error of the skewness viation
(d) Error the kurtosis
Figure 12: Comparison of moments for the switchpoint model.
1.2
UP 24C
Lis SMC
0.3-
4-n
044
0.3
03
25
23
20
30
(a) Error of the mean
(b) Error of the standard de- (c) Error of the skewness
(d) Error of the kurtosis
viation
Figure 13: Comparison of moments for the mixture model.
vear coeors bwattn
+- IMONN
70
DNC
Ca
C2
32
Tiam
(a) Error the mean
(b) Error the standard de- (c) Error of the skewness
(d) Error the kurtosis
viation
Figure 14: Comparison of moments for the HMM.
-SHL 11000 sampasi -urs (ecert]
0.05
SMAL 11:00 samplss
(eocet)
30%% (PONTT)
samples)
0.08
0.04
prebeaiity
0.06
0.02
0.01
0.01
0.00
0.00
0.00
200
10
20
resul:
resul
(a) Original population model.
(b) Modified population model.
(c) Two-type population model.
0.30
L.M 11.000 sampleess toxact)
0.150
sampless (anect)
-SML 11.000 sampless
lexect)
0.125
0.10
brobabitty
0.20
0.200
0.08
0.075
0.08
6.10
0.050
0.04
0n
6.00
0.000
6.00
(d) Switchpoint model.
(c) Mixture model.
(f) Hidden Markov model.
Figure 15: Histograms of the exact distribution and the MCMC samples with the lowest TVD and similar running time to our exact method.
35

---

Contents
1 Introduction
2 Bayesian Probabilistic Programming
3 Generating Functions
3.1 Translating programs to generating functions
4
Implementation & Optimizations
5
Empirical Evaluation
5.1
Comparison with exact inference methods
8
5.2
Comparison with approximate inference methods
5.3
Benchmarks with infinite-support distributions.
6
Conclusion
10
Details on the Probabilistic Programming L Language SGCL
14
A.1 A minimal grammar
14
A.2 Syntactic sugar
14
A.3 Standard measure semantics
14
B Details on Generating Functions
15
B.1
Generating function semantics
17
B.2 Correctness proof
17
B.3 Possible extensions to the probabilistic programming language
22
Details on implementation and optimizations
23
C.I Automatic differentiation and Taylor polynomials
24
C.2 Optimizing Observations from Compound Distributions.
26
C33 Analysis of the running time
29
C,4 Numerical issues
30
D Details on the Empirical Evaluation
30
D.1 Comparison with exact inference
30
D.2 Comparison with approximate inference
31
D.3 Benchmarks with infinite support
31
D.4 Additional results
34
36


Reviews:
[{"review_id": "wZjBIVNFxB", "reviewer": "Reviewer_c91c", "summary": "A probabilistic programming framework for exact inference in discrete probabilistic programs with infinite support is introduced. The framework levegares probability generating functions as a first class attribute of discrete distributions. The probabilistic programming language is formally defined, inference is theoretically analysed for correctness. The performance and accuracy are empirically evaluated and compared with baselines on real-world examples.", "questions": "You apparently unroll loops in your benchmarks (I assume because the language does not support loops). What prevents you from implementing restricted loops in the language?\n\nThe syntax in the examples is richer than the BNF in the paper. What is '+~' for example? Can you define the full syntax in the paper, at least briefly?\n\nWhy rust? Have you considered other options?\n\n", "limitations": "The paper adequately discusses limitations of the proposed solution.", "rating": 8, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper brings the idea of using GF as a first-class property of a distribution into a practical implementation as a probabilistic programming language. Formal properties of the language and inference are rigorously analysed. The empirical evaluation demonstrates capabilities of the approach, and favorably compares to alternatives. The submission is accompanied by a working and well organized code based.", "weaknesses": "The PPL is implemented in Rust, a language which, while it certainly has its merits, less popular for statistical modelling and Bayesian inference, which may make adoption of the approach more difficult. This is a matter of taste though.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The PPL is implemented in Rust, a language which, while it certainly has its merits, less popular for statistical modelling and Bayesian inference, which may make adoption of the approach more difficult.\n\nIt is correct that Rust is not used in Bayesian inference very much. However, it would not be too difficult to implement, say, a Python API calling Rust under the hood. Many Python libraries work that way (e.g. PyTorch calling C++). Developing the Python interface is mostly an engineering problem, whereas we wanted to focus on the research questions and a proof of concept in our work.\n\n> What prevents you from implementing restricted loops in the language?\n\nThere is no particular reason why we couldn’t have implemented bounded loops in our language. But since this does not affect the expressiveness of the language, we preferred the conceptual simplicity of omitting them.\n\n> The syntax in the examples is richer than the BNF in the paper. What is '+~' for example? Can you define the full syntax in the paper, at least briefly?\n\nYou are correct that the examples use syntax that is slightly richer than the one provided in the paper. We introduced additional constructs to reduce the number of temporary variables needed; this does not affect the expressivity of the language. There are three additional constructs:\n\n* `X +~ D;` stands for `TMP ~ D; X := X + TMP`\n* `X += a * Y + b;` stands for `TMP := a * Y + b; X := X + TMP;`)\n* `if n ~ D { … } else { … }` stands for `TMP ~ D; if TMP = n { … } else { … }`\n\nWe will add an explanation of these constructs to the paper.\n\n> Why Rust? Have you considered other options?\n\nThe main reasons were low-level control and performance. The operations on the Taylor polynomials (to evaluate derivatives of the generating function) need to be fast and we have a few optimizations that exploit the structure of GFs arising from probabilistic programs. C or C++ would have satisfied this criterion as well, but Rust’s language features like memory safety, enums (tagged unions), and pattern matching made the implementation a lot more pleasant, robust, and easier to maintain. The first author’s experience with Rust was another contributing factor.\n"}, {"review_id": "Ldx8YnpsLo", "reviewer": "Reviewer_hi49", "summary": "The paper presents a restricted, first-order probabilistic programming language in which inference for infinite-support latent variables, with continuous or discrete sample spaces, admits a closed form.  For experiments they compare against importance sampling, Lightweight M-H, Random Walk M-H (RMH), sequential Monte Carlo, particle Gibbs sampling, and interacting particle MCMC (IPMCMC) in the Anglican universal PPL.  The authors achieve much lower error than Monte Carlo sampling across experiments, typically with significantly less compute.", "questions": "How could the exact inference algorithm in this paper be used in a larger probabilistic programming framework, such as when probabilistic circuits admitting exact density calculations are applied as proposals in a broader family of inference algorithms?  When the paper mentions outputting posterior probability masses for certain intervals of integer-valued latent variables, why not also include densities for continuous latent variables?\n\nThe authors have addressed these questions in their response.", "limitations": "The authors have admitted the limitations of their work, and in exchange for those limitations in expressivity they gain performance in both error and compute time.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The paper demonstrates that when latent variables are subject only to affine transformations, a first-order (lacking general recursion) probabilistic programming language can admit tractable exact inference for certain important classes of discrete observations. Even when the resulting exact inference engine takes compute time comparable to the Anglican Monte Carlo inference routines, it tends to run in the middle of the pack while providing much lower approximation error.", "weaknesses": "Since the paper focuses on proving and benchmarking something that has an exact answer, it has few weaknesses as such.  It could perhaps have done a better job of making explicit the if-and-only-if conditions that must be met for generating functions to support closed-form inference.\n\nThe authors have addressed this concern in the rebuttal.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> making explicit the if-and-only-if conditions that must be met for generating functions to support closed-form inference\n\nPlease refer to our general rebuttal for the “only-if” aspect (whether additional constructs could be added to the programming language while still preserving a closed form). For the “if”-direction, the syntax of our language guarantees that the generating function has a closed form.\n\n> How could the exact inference algorithm in this paper be used in a larger probabilistic programming framework?\n\nPlease refer to the general rebuttal for this point as well.\n\n> When the paper mentions outputting posterior probability masses for certain intervals of integer-valued latent variables, why not also include densities for continuous latent variables?\n\nWhile there are mathematical ways of recovering the probability density function from a generating function via an inverse Laplace transform, this cannot be automated in practice because it requires solving integrals, which is intractable in general. So outputting posterior densities is not possible, as far as we know.\n"}, {"review_id": "zN5K8tSxSd", "reviewer": "Reviewer_SkQ5", "summary": "This paper proposes a new probabilistic programming language called statistical guarded command language (SGCL) and an exact Bayesian inference method for any discrete statistical model that can be expressed in SGCL. Empirical evaluatios on a variety of different discrete probabilistic models demonstrate the effectiveness of the method, especially for challenging distributions with discrete variables that have infinite support.", "questions": "See weaknesses.", "limitations": "The authors adequately addressed the limitations.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper is well organized and clearly written.\n- The proposed method seems novel, and is very useful especially for challenging discrete distributions involving variables with infinite support.\n", "weaknesses": "- The paper seems to miss discussions on probabilistic circuits which is a large body of work concerning probabilistic models with tractable inference.\n- Given the limitations of the proposed PPL and its specialized use cases, it would also be helpful to add some discussions on how this can potentially be integrated into/work together with other more general PPLs to expand the applicability of the proposed method.\n- The current title and abstract seems slightly misleading, in the sense that the proposed PPL is not generally applicable to any discrete distributions. It would be better to rephrase to make the scope of its applicability clearer, and also emphasize its targeted use case is distributions with variables having infinite support.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The paper seems to miss discussions on probabilistic circuits which is a large body of work concerning probabilistic models with tractable inference.\n\nThank you for pointing this out. Probabilistic circuits (Choi et al., 2020) are indeed a family of probabilistic models that allow for efficient and exact inference, unifying and generalizing many previously introduced tractable probabilistic models (e.g. arithmetic circuits, sum-product networks, and cutset networks). Two exact inference systems that we cite are connected to probabilistic circuits: SPPL (Saad et al., PLDI 2021) extends sum-product networks, and Dice (Holtzen et al., OOPSLA 2020) is based on Weighted Model Counting on binary decision diagrams. After the NeurIPS submission deadline we even learned about work applying generating functions to probabilistic circuits, called “Probabilistic Generating Circuits” (Zhang et al., ICML 2021; Harviainen et al., UAI 2023). It should be noted, however, that the latter circuits only support binary random variables, whereas we support variables with infinite support. We will add a discussion of probabilistic circuits, and probabilistic generating circuits in particular, to the related works section of the paper.\n\n> how this can potentially be integrated into/work together with other more general PPLs\n\nPlease refer to the general rebuttal regarding this point.\n\n> The current title and abstract seems slightly misleading, in the sense that the proposed PPL is not generally applicable to any discrete distributions. It would be better to rephrase to make the scope of its applicability clearer, and also emphasize its targeted use case is distributions with variables having infinite support.\n\nIt is true that the title is not a complete description of our setting, but we believe the title is quite long already, so we would prefer not to add extra clarifications. We’ve made an effort to be as transparent about the limitations and restrictions about our approach as possible (paragraph “Limitations” in the introduction and “Restrictions” in Section 2). In the abstract, we are explicit about which language constructs are supported, so we did not feel this was misleading. We could add a sentence listing unsupported constructs, but this would just be a negation of the supported constructs (no continuous observations, no non-affine functions, etc.), so it seems redundant. There are subtleties (such as which compound distributions we support), but we believe that this would be too detailed for an abstract and we discuss it in the above paragraphs in the paper. However, we will clarify in the abstract that we only support observing events involving *discrete* variables and that our method is not applicable to all discrete distributions."}, {"review_id": "vk09IPKVuX", "reviewer": "Reviewer_JaNH", "summary": "This paper proposes a new approach to exact inference in probabilistic programs. Existing approaches to exact inference attempt to compute marginal and normalized probabilities by efficiently computing the necessary (finite) sums, or using computer algebra to recognize and solve tractable integrals. By contrast, the present work constructs the moment-generating function ($G(x) = \\mathbb{E}_{X \\sim \\mu}[x^X]$) for the distribution represented by the program, and takes its (higher-order) derivatives (using AD) to compute exact posterior marginal probabilities, as well as exact moments. This enables automated exact inference in a broader class of models than handled by existing techniques. The approach is implemented in Rust, and experiments in four example models show that their implementation typically computes exact answers in less time than it takes for Anglican’s generic Monte Carlo algorithms (e.g., sequential Monte Carlo with the prior as the proposal) to converge to accurate estimates.", "questions": "1. Are you 100% sure that systems like Hakaru and $\\lambda$PSI cannot handle your examples? I don’t have a great handle on their capabilities, but I expect Maple’s computer algebra simplifications can handle at least some infinite sums and integrals. Have you tried implementing the example programs and seeing what these tools spit out?\n2. What are the prospects for a generating-function semantics of a functional (rather than imperative) source language, such as Anglican (but suitably restricted)?", "limitations": "The paper is up-front about both the expressiveness limitations of the probabilistic language they support, and the factors in which their algorithm scales exponentially.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The key strengths I see in this paper are:\n\n- It introduces a significant new approach to PPL inference (automatic differentiation on compositionally constructed moment-generating functions), and I believe it’s important enough that most people working in PPLs will want to understand the new technique. The exposition is clear and should be accessible to a broad probabilistic programming audience.\n- The generating function approach resolves a key limitation of other exact inference methods: it can marginalize latent variables with infinite support.\n- Limitations of the language, and performance concerns, are clearly outlined (e.g., in the Restrictions paragraph on L147, and the Limitations paragraph on L83).\n- The appendix contains rigorous correctness proofs and further exposition of the generating function semantics. There is also a Rust implementation that helps clarify the actual implementation of the technique (e.g., the data representation of Taylor polynomials used by the AD algorithm).", "weaknesses": "The key weaknesses I notice are:\n\n- The novel capabilities (and limitations) with respect to the line of work by Klinkenberg et al. could be more clearly described. It would also be good to cite and explain the differences with their recent LAFI abstract, “Exact Probabilistic Inference Using Generating Functions” (https://arxiv.org/pdf/2302.00513.pdf). That workshop abstract is non-archival and much less thorough than the present work, so I don’t mean this as a critique of your work’s novelty. But I think readers will appreciate a clear-eyed description of exactly how their proposed approach differs from yours.\n- It is unclear to what extent the limitations of the language are fundamental limitations of the generating function technique, vs. concessions made to get to a simple first publication/implementation of the idea. Could the technique be extended to programs with loops, or does that seem unlikely? What about compound distributions with more than one random parameter (e.g., Binomial where both $n$ and $p$ is also random)? What about ‘score statements’ that multiply a non-negative (affine?) expression into the likelihood? (It seems this should be possible for scores < 1, because we can `observe 1 ~ Bernoulli(X)`, but can it be made to work for more general scores?)\n- The experiments are all on programs for which existing exact inference methods fail. But it would be nice to understand how the new approach compares to existing methods on the programs that both approaches can handle. (Even if it is much less efficient, I believe the method still warrants publication, but it would be useful to readers to understand the expressiveness / efficiency trade-offs clearly.)\n- L310-312 incorrectly state that WebPPL does not support SMC or particle MCMC. See the SMC and PMCMC implementations here: https://github.com/probmods/webppl/tree/master/src/inference.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The novel capabilities (and limitations) with respect to the line of work by Klinkenberg et al. could be more clearly described. It would also be good to cite and explain the differences with their recent LAFI abstract\n\nThe ideas alluded to in their LAFI abstract are now available in a more fleshed-out form as a preprint (https://arxiv.org/abs/2307.07314). (Note that this was posted on arXiv after the NeurIPS submission deadline.) The main difference is that their work only supports discrete distributions and yields a symbolic representation of the GF for the posterior distribution. They also have some support for loops, which requires the user to provide a loop invariant template. On the other hand, our approach allows for some continuous distributions and is (usually) much more performant and scalable, because we only evaluate the derivatives of the generating function instead of representing them fully symbolically (cf. the comparison with PSI and Prodigy in the general rebuttal). We will add a discussion of this to the paper.\n\n> It is unclear to what extent the limitations of the language are fundamental limitations of the generating function technique, vs. concessions made to get to a simple first publication/implementation of the idea.\n\nWe give a high-level answer to this in the general rebuttal. Details on your specific questions follow here.\n\n> Could the technique be extended to programs with loops, or does that seem unlikely?\n\nSupporting exact inference for programs with loops seems extremely difficult, unless the inference algorithm is given additional information, e.g. a loop invariant (see the work of Klinkenberg et al.). However, finding such a loop invariant for nontrivial programs seems exceedingly difficult. Furthermore, the variables of looping programs may have infinite expected values, so the generating function may not be definable at the point 1. But in our semantics, evaluating at 1 (for marginalization) is very important. So there are many challenges and opportunities for future research here.\n\n> What about compound distributions with more than one random parameter (e.g., Binomial where both n and p is also random)?\n\nIt is hard to say with certainty (somebody might find a closed form that nobody expected), but we tried very hard to find additional compound distributions that admit a closed form and were unsuccessful. The only additional compound distributions we can support are `Binomial(n, X)` and `Gamma(X, a)`, as mentioned in the general rebuttal.\n\n> What about ‘score statements’ that multiply a non-negative (affine?) expression into the likelihood? It seems this should be possible for scores < 1, because we can observe 1 ~ Bernoulli(X), but can it be made to work for more general scores?\n\nAs mentioned in the general rebuttal, `score a^X` for nonnegative real `a` < 1 and variable `X` should be supportable. (For `a` > 1, this runs into problems: the GF may not be defined at 1 anymore.) Nonnegative affine scoring (or even nonnegative polynomial scoring) should be possible as well, using a similar technique as for observing from `Bernoulli(X)`. We didn’t include this because we did not find any real-world example to use it.\n\n> But it would be nice to understand how the new approach compares to existing methods on the programs that both approaches can handle.\n\nIn the general rebuttal, we attached a comparison with Klinkenberg’s approach to generating functions and the PSI solver.\n\n> L310-312 incorrectly state that WebPPL does not support SMC or particle MCMC\n\nThank you very much for pointing this out. We will of course correct this statement in the paper. It seems that the documentation of WebPPL’s inference algorithms (http://docs.webppl.org/en/master/inference/methods.html) is incomplete and there are more undocumented algorithms available (such as Particle MCMC).\n\n> Are you 100% sure that systems like Hakaru and λPSI cannot handle your examples?\n\nWe tested λPSI and it could not solve any of our examples. We did not test Hakaru because it is quite difficult to get running, requires a Maple license, and benchmarks in the PSI paper indicated that PSI was more performant and able to solve more benchmarks than Hakaru. Even if Hakaru is able to solve small examples, as a fully symbolic tool, it is unlikely to scale to our examples given that even Prodigy (which already exploits the generating function structure) times out.\n\n> What are the prospects for a generating-function semantics of a functional (rather than imperative) source language, such as Anglican (but suitably restricted)?\n\nWe believe it should be possible to define a generating-function semantics for a functional language (after all, a functional program can be translated to an imperative one), but it is highly nontrivial: semantics for functional languages usually rely on composition, but the generating functions describe the joint distribution of all variables in the program at any point – we never reason “locally” about only a subset of the variables because this would lose information. So this is an open research question.\n"}]

Author rebuttal:
First of all, we would like to thank the reviewers for their thoughtful reviews and interesting questions about our work. We are delighted that the reviewers think that our approach is “significant” (JaNH), “novel and very useful” (SkQ5), that our analysis is “rigorous” (c91c, JaNH), that our implementation compares favorably”in the experimental evaluation (c91c, hi49), and that our presentation is “well-organized and clearly written” (SkQ5). We will try to answer the reviewers’ questions in the following.

In this general rebuttal, we will address points that came up in more than one review or are of interest to more than one reviewer. Specific questions will be discussed in the individual responses.

**Integration with more general probabilistic programming systems (reviewers hi49 & SkQ5)**: Given the restrictions on the programs where our method can be applied, it is a natural question to ask how it could be integrated with a more general probabilistic programming system. We believe this should be possible. For example, in a sampling-based inference algorithm, one could imagine using generating functions to solve subprograms exactly if this is possible (rather than sampling the variables in those subprograms as well). Investigating the details and how well such an approach would work in practice is an interesting research question.

**Conditions for closed-form generating functions (reviewers hi49 & JaNH)**: we were asked for exact conditions under which the generating function of a probabilistic program admits a closed form. The answer is that we tried to make the programming language in the paper (SGCL) as expressive as possible while still maintaining a closed form for the generating function. There are only the following possible additions to the language (that we know of) that preserve closed forms for the GFs but that we didn’t include:

* more distributions with constant parameters: as long as its generating function exists and that function and its derivatives can be evaluated, the distribution can be supported.
* the compound distribution Binomial(n, X), where the success probability is a variable, but this is already expressible as the sum of n independent Bernoulli(X) variables.
* the compound distribution Gamma(X, a), where X is the shape parameter
* the constructs `X mod 2 = 0`, `X := max(X - 1, 0);`, `X := sum_iid(D, Y);` where `D` is a primitive distribution and `X`, `Y` nonnegative integer variables, as mentioned in Klinkenberg et al.’s work
* `score(a^X);` statements for `a` < 1 or `score(q(X));` where `q` is a polynomial with nonnegative coefficients should also be supportable. (A `score` statement multiplies the likelihood by the given expression).

We did not include these constructs because we could not find a real-world example that can use them. It is obviously impossible to say with certainty that no other constructs preserve closed forms (because there are infinitely many possible constructs we might not have thought of). But we tried very hard to express more concepts (e.g. observations from continuous distributions) in terms of generating functions but could not find anything else that works. We will add these remarks to the appendix of our paper.

**Comparison with existing exact tools (reviewer JaNH)**: Reviewer JaNH said it would be interesting to compare our new approach with existing ones on benchmarks that both methods can handle. Accordingly, we decided to compare the PSI tool for exact symbolic inference with our tool on benchmarks that both can handle (i.e. finite discrete distributions). Furthermore, the Prodigy tool by Klinkenberg et al. (which also uses generating functions, but in combination with computer algebra instead of automated differentiation) recently gained the ability to deal with conditioning as well. (A preprint of this work was posted to the arXiv after the NeurIPS submission deadline: https://arxiv.org/abs/2307.07314.) Thus, we decided to compare PSI, Prodigy, and our tool on the subset of Prodigy’s benchmarks that all three tools can handle. The results can be found in the attached PDF. As you can see, our tool is very much competitive and often fastest on these examples. Furthermore, we tried running Prodigy on our benchmarks, but it either does not support them (because of continuous priors), runs out of memory, or times out after 30 min (whereas our tool finishes in seconds). This confirms our claim about the better scalability of our method compared to fully symbolic approaches.

Meta-review: This paper introduces probabilistic program inference by means of probability generating functions.
Reviewers agree the analysis is rigorous, and the idea is novel in the probabilistic programming context, even if there are some limitations to the technique. Overall a very solid contribution to the field.
Decision: Accept (oral)