<figure><img src="image_1.png" /></figure>

HAL
open science

# Template-Based Static Posterior Inference for Bayesian Probabilistic Programming

Peixin Wang, Hongfei Fu, Tengshun Yang, Guanyan Li, Luke Ong

<figure><img src="image_2.png" /></figure>

## To cite this version:

Peixin Wang, Hongfei Fu, Tengshun Yang, Guanyan Li, Luke Ong. Template-Based Static Posterior Inference for Bayesian Probabilistic Programming. 2023. hal-04022797v2

HAL Id: hal-04022797

https://hal.science/hal-04022797v2

Preprint submitted on 24 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Template-Based Static Posterior Inference for Bayesian
Probabilistic Programming

PEIXIN WANG, University of Oxford, United Kingdom

HONGFEI FU*, Shanghai Jiao Tong University, China

TENGSHUN YANG, SKLCS, Institute of Software, Chinese Academy of Sciences & University of Chinese
Academy of Sciences, China

GUANYAN LI, University of Oxford, United Kingdom

LUKE ONG, University of Oxford, United Kingdom & Nanyang Technological University, Singapore

In Bayesian probabilistic programming, a central problem is to estimate the normalised posterior distribution (NPD) of a probabilistic program with *score* (a.k.a. *observe*) statements. Prominent approximate approaches to address this problem include Markov chain Monte Carlo (MCMC) and variational inference (VI), but neither of them is capable of generating guaranteed outcomes within a finite time limit. Moreover, most existing formal approaches that perform exact inference for NPD are restricted to programs with closed form solutions to NPD or bounded loops/recursion. A recent work (Beutner et al., PLDI 2022) proposed an automated approach that derives guaranteed bounds for NPD over programs with unbounded recursion. However, as this approach requires recursion unrolling, it suffers from the path explosion problem. Furthermore, previous approaches do not consider *score-recursive* probabilistic programs that allow score statements inside loops, which is non-trivial and requires careful treatment to ensure the integrability of the normalising constant in NPD.

In this work, we propose a novel automated approach to derive bounds for NPD via polynomial templates, fixed-point theorems and Optional Stopping Theorem (OST). Our approach can handle probabilistic programs with unbounded while loops and continuous distributions with infinite supports. The novelties in our approach are three-fold: First, we use polynomial templates to circumvent the path explosion problem from recursion unrolling; Second, we derive a novel multiplicative variant of OST that addresses the integrability issue in score-recursive programs; Third, to increase the accuracy of the derived bounds via polynomial templates, we propose a novel technique of truncation that truncates a program into a bounded range of program values. Experiments over a wide range of benchmarks demonstrate that our approach is time-efficient and can derive bounds for NPD that is comparable with (or even tighter than) the recursion-unrolling approach (Beutner et al., PLDI 2022).

CCS Concepts: • **Theory of computation** → **Logic and verification; Automated reasoning; Program verification.**

Additional Key Words and Phrases: Probabilistic programs, Bayesian inference, Quantitative verification, Martingales, Fixed-point Theory, Posterior distributions

# 1 INTRODUCTION

In probabilistic programming, Bayesian statistical probabilistic programming adds specific statements for Bayesian reasoning into probabilistic programs. The principle of Bayesian probabilistic programming aims at first modelling probabilistic models as probabilistic programs and then analyzing the models through their program representation. Compared with traditional approaches

Authors' addresses: Peixin Wang, University of Oxford, United Kingdom, peixin.wang@cs.ox.ac.uk; Hongfei Fu, Shanghai Jiao Tong University, China, fuhf@cs.sjtu.edu.cn; Tengshun Yang, SKLCS, Institute of Software, Chinese Academy of Sciences & University of Chinese Academy of Sciences, China, yangts@ios.ac.cn; Guanyan Li, University of Oxford, United Kingdom, guanyan.li@cs.ox.ac.uk; Luke Ong, University of Oxford, United Kingdom & Nanyang Technological University, Singapore, luke.ong@ntu.edu.sg.

<footnote>*Corresponding author</footnote>
[8, 10, 30, 31] that specify an ad-hoc language, probabilistic programming languages (PPLs) [48] provide a universal framework to perform Bayesian inference. Unlike standard programming languages, PPLs have two specific constructs: sample and score [7].<fnref n="1" /> The former construct allows drawing samples from a (prior) distribution, while the latter records the likelihood of observed data in the form of “score(weight)”.<fnref n="2" /> The universality of PPLs has the potential of enabling scientists and engineers to design and explore sophisticated probabilistic models easily: when using these languages, they no longer need to worry about developing custom inference engines for their models, which is previously a highly-nontrivial task that requires expertise in both statistics and machine learning. Thanks to its universality, Bayesian probabilistic programming has nowadays become an active research subject in both machine learning and programming language communities, and there have been an abundance of PPLs for Bayesian inference, such as Pyro[5], WebPPL[22], Anglican[47], Church[21], etc.

In this work, we consider the central problem of analyzing the normalised posterior distribution (NPD) in Bayesian inference over probabilistic programs. The general setting of this problem is that given a prior distribution $p(z)$ over the latent variables $z \in \mathbb{R}^n$ of interest and the distribution $p(x, z)$ of the probabilistic model represented by a probabilistic program, the task is to estimate the NPD by observing the evidence $x \in \mathbb{R}^m$ with the likelihood $p(x|z)$. Note that the problem can be generally solved by the Bayes' rule $p(z|x) = \frac{p(x|z)p(z)}{p(x)}$, but the main difficulty to apply Bayes' rule is that the normalising constant $p(x) = \int p(x|z)p(z)dz$ is usually intractable to compute.

In the literature, there are two classes of approaches to address the NPD problem. The first is the approximate approaches that estimate the NPD by random simulations, while the second is the formal approaches that aim at deriving guaranteed bounds for NPD. In approximate approaches, two dominant methods are Markov chain Monte Carlo (MCMC) [18] and variational inference (VI) [6]. Although approximate approaches can produce approximate results efficiently, they cannot provide formal guarantee within a finite time limit. Moreover, as shown in [4], approximate approaches may produce inconsistent results between different simulation methods. In formal approaches, there is a large amount of previous works such as (λ)PSI [19, 20], AQUA [25], Hakaru [33] and SPPL [41], aiming to derive exact inference for NPD. However, these methods are restricted to specific kinds of programs, e.g., programs with closed-form solutions to NPD or without continuous distributions, and none of them can handle probabilistic programs with unbounded while-loops/recursion. Recently, Beutner et al. [4] proposed an approach that infers guaranteed bounds for NPD and allows unbounded recursion. The main techniques in this approach are (i) the unrolling of every recursion in a probabilistic program to eliminate the appearance of recursion, (ii) the widening operator of abstract interpretation to eliminate the non-termination case in the unrolling and (iii) the interval semantics to handle continuous distributions.

*Challenges and gaps.* In this work, we focus on developing formal approaches to derive bounds for NPD over probabilistic programs. In the previous formal approaches to address this problem, the most relevant work is Beutner et al. [4], but it has the following drawbacks. The first is that this approach is based on recursion unrolling, and hence may cause path explosion. The second is that this approach cannot handle the situation where score statements with weight greater than 1 appear inside a loop. In the sequel, we call such programs *score-recursive*, and show that score inside a loop may cause unbounded weights and integrability issues and thus requires careful treatment. Staton et al. [45] also noted that for a non-recursive λ-calculus with score, unbounded weights may introduce the possibility of “infinite model evidence errors”. To circumvent the second

<footnote n="1">¹Sometimes observe is used instead of score [23], which has the same implicit effect.</footnote>

<footnote n="2">²The argument “weight” corresponds to the likelihood each time the data is observed.</footnote>
drawback, previous results (e.g., Borgström et al. [7]) allow only 1-bounded weights, and no existing approaches can handle score-recursive programs whose weight of score can be greater than 1.

*Template-based approaches.* To avoid the path explosion problem arising from the recursion unrolling method [4], we consider the template paradigm [9, 10, 13] that first sets up a numerical template for the bound function to be solved, then establishes constraints for the template which is derived from the underlying verification theory (e.g. fixed point theory, Optional Stopping Theorem, etc.), and finally solves the template based on the derived constraints to obtain a concrete bound function. Template-based approaches have been shown to be capable of synthesizing tight bounds for expectation and probabilistic properties over probabilistic programs [50, 53], and are time efficient since they can leverage highly optimized constraint solvers for linear, semidefinite, and convex programming, etc. Therefore, the template paradigm serves as a good candidate to address the NPD problem in Bayesian probabilistic programming. However, simply applying the template paradigm does not suffice by the following reasons:

* The most accurate numerical template by far is polynomial template. However, polynomial template is tight only over a bounded region (see e.g. Weierstrass Approximation Theorem [26]), and in general is not accurate over an unbounded region.

* The common pattern for a template-based approach is to synthesize a single bound function via a single process of template solving. However, having a single template is not enough for synthesizing tight bounds for NPD, as different initial program inputs may need different bound functions to achieve tight bounds.

In this work, we have careful treatment to address the two weakness above. Our detailed contributions are as follows.

*Our contributions.* In this work, we present the following contributions:

* First, to circumvent the path explosion from loop unrolling (that corresponds to recursion unrolling in functional programs) in the previous work [4], we propose a novel approach that uses polynomial templates [11, 12, 53], fixed-point theorems [46] and Optional Stopping Theorem (OST) [54] to synthesize polynomial bounds for NPD.

* Second, to address the integrability issue from score-recursive programs, we present a novel multiplicative variant of the classical OST to derive expectation bounds. This addresses the challenge of handling score-recursive probabilistic programs stated previously.

* Third, to increase the accuracy of the derived bounds, we propose a novel truncation operation that truncates the probabilistic program of concern onto a bounded range of program values, which allows our algorithm to explore high-degree polynomials to improve the accuracy of the derived bounds. Moreover, we devise our algorithm in the way to synthesize multiple bounds for various initial program inputs to further improve the accuracy. The truncation operation and the synthesis of multiple bounds address the weakness of the existing template-based approaches in the literature.

* Finally, experimental results show that our approach can handle a wide range of benchmarks including non-parametric examples such as Pedestrian [4] and score-recursive examples such as phylogenetics [38]. Compared with the previous approach [4], our approach reduces the runtime by up to 6 times, while deriving comparable or even tighter bounds for NPD.

*Technical contributions.* Compared with the previous results on the use of template paradigm and OST, we propose a truncation operation (onto a bounded range of program values) that leverages the approximability of polynomials over a bounded region, and a multiplicative variant of OST that can handle score-recursive programs.
*Limitations*. One limitation is that our approach has the combinatorial explosion that when the degree of the polynomial template increases, the time of the template solving increases exponentially. However, from our experimental results, a moderate choice of the degree (e.g., no more than 10) suffices. Another is that for sampling statements with continuous distributions, our approach currently assumes polynomial density functions. However, with proper approximation via (piecewise) polynomials, our approach can handle any continuous density functions.

## 2 PRELIMINARIES

We first review some basic concepts from probability theory (see standard textbooks such as [35, 54] for a detailed treatment), and then present the Bayesian probabilistic programming language and the normalised posterior distribution (NPD) problem. Throughout the paper, we denote by $\mathbb{N}$, $\mathbb{Z}$ and $\mathbb{R}$ the sets of all natural numbers (including zero), integers, and real numbers, respectively.

### 2.1 Basics of Probability Theory

A measurable space is a pair $(U, \Sigma_U)$, where $U$ is a nonempty set and $\Sigma_U$ is a $\sigma$-algebra on $U$, i.e., a family of subsets of $U$ such that $\Sigma_U \subseteq \mathcal{P}(U)$ contains $\emptyset$ and is closed under complementation and countable union. Elements of $\Sigma_U$ are called measurable sets. A function $f$ from a measurable space $(U_1, \Sigma_{U_1})$ to another measurable space $(U_2, \Sigma_{U_2})$ is measurable if $f^{-1}(A) \in \Sigma_{U_1}$ for all $A \in \Sigma_{U_2}$.

A measure $\mu$ on a measurable space $(U, \Sigma_U)$ is a mapping from $\Sigma_U$ to $[0, \infty]$ such that (i) $\mu(\emptyset) = 0$ and (ii) $\mu$ is countable additive: for every pairwise-disjoint set sequence $\{A_n\}_{n \in \mathbb{N}}$ in $\Sigma_U$, it holds that $\mu(\bigcup_{n \in \mathbb{N}}} A_n) = \sum_{n \in \mathbb{N}} \mu(A_n)$. We call the triple $(U, \Sigma_U, \mu)$ a measure space. If $\mu(U) = 1$, we call $\mu$ a probability measure, and $(U, \Sigma_U, \mu)$ a probability space. The Lebesgue measure $\lambda$ is the unique measure on $(\mathbb{R}, \Sigma_{\mathbb{R}})$ satisfying $\lambda([a, b)) = b - a$ for all valid intervals $[a, b)$ in $\Sigma_{\mathbb{R}}$. For each $n \in \mathbb{N}$, we have a measurable space $(\mathbb{R}^n, \Sigma_{\mathbb{R}^n})$ and a unique product measure $\lambda_n$ on $\mathbb{R}^n$ satisfying $\lambda_n(\prod_{i=1}^n A_i) = \prod_{i=1}^n \lambda(A_i)$ for all $A_i \in \Sigma_{\mathbb{R}}$.

The Lebesgue integral operator $\int$ is a partial operator that maps a measure $\mu$ on $(U, \Sigma_U)$ and a real-valued function $f$ on the same space $(U, \Sigma_U)$ to a real number or infinity, which is denoted by $\int f d\mu$ or $\int f(x)\mu(dx)$. The detailed definition of Lebesgue integral is somewhat technical, see [37, 40] for more details. Given a measurable set $A \in \Sigma_U$, the integral of $f$ over $A$ is defined by $\int_A f(x)\mu(dx) := \int f(x) \cdot [x \in A]\mu(dx)$ where $[-]$ is the Iverson bracket such that $[\phi] = 1$ if $\phi$ is true, and $0$ otherwise. If $\mu$ is a probability measure, then we call the integral as the expectation of $f$, denoted by $\mathbb{E}_{X \sim \mu; A}[f]$, or $\mathbb{E}[f]$ when the scope is clear from the context.

For a measure $v$ on $(U, \Sigma_U)$, a measurable function $f: U \to \mathbb{R}_{\ge 0}$ is the density of $v$ with respect to $\mu$ if $v(A) = \int f(x) \cdot [x \in A]\mu(dx)$ for all measurable $A \in \Sigma_U$, and $\mu$ is called the reference measure (most often $\mu$ is the Lebesgue measure). Common families of probability distributions on the reals, e.g., uniform, normal distributions, are measures on $(\mathbb{R}, \Sigma_{\mathbb{R}})$. Most often these are defined in terms of probability density functions with respect to the Lebesgue measure. That is, for each $\mu_D$ there is a measurable function $pdf_D: \mathbb{R} \to \mathbb{R}_{\ge 0}$ that determines it: $\mu_D(A) := \int_A pdf_D(d\lambda)$. As we will see, density functions such as $pdf_D$ play an important role in Bayesian inference.

Given a probability space $(\Omega, \mathcal{F}, \mathbb{P})$, a random variable is an $\mathcal{F}$-measurable function $X: \Omega \to \mathbb{R} \cup \{+\infty, -\infty\}$. The expectation of a random variable $X$, denoted by $\mathbb{E}(X)$, is the Lebesgue integral of $X$ w.r.t. $\mathbb{P}$, i.e., $\int_X d\mathbb{P}$. A filtration of $(\Omega, \mathcal{F}, \mathbb{P})$ is an infinite sequence $\{\mathcal{F}_n\}_{n=0}^{\infty}$ such that for every $n \ge 0$, the triple $(\Omega, \mathcal{F}_n, \mathbb{P})$ is a probability space and $\mathcal{F}_n \subseteq \mathcal{F}_{n+1} \subseteq \mathcal{F}$. A stopping time w.r.t. $\{\mathcal{F}_n\}_{n=0}^{\infty}$ is a random variable $T: \Omega \to \mathbb{N} \cup \{0, \infty\}$ such that for every $n \ge 0$, the event $\{T \le n\}$ is in $\mathcal{F}_n$.

A discrete-time stochastic process is a sequence $\Gamma = \{X_n\}_{n=0}^{\infty}$ of random variables in $(\Omega, \mathcal{F}, \mathbb{P})$. The process $\Gamma$ is adapted to a filtration $\{\mathcal{F}_n\}_{n=0}^{\infty}$, if for all $n \ge 0$, $X_n$ is a random variable in $(\Omega, \mathcal{F}_n, \mathbb{P})$. A discrete-time stochastic process $\Gamma = \{X_n\}_{n=0}^{\infty}$ adapted to a filtration $\{\mathcal{F}_n\}_{n=0}^{\infty}$ is a martingale (resp. supermartingale, submartingale) if for all $n \ge 0$, $\mathbb{E}(|X_n|) < \infty$ and it holds almost surely (i.e., with
start := sample uniform(0, 3);
pos := start; dist := 0;
l_init: while pos ≥ 0 do
    step := sample uniform(0, 1);
    if prob(0.5) then
        pos := pos - step
    else
        pos := pos + step
fi;
dist := dist + step
od;
score (normal (1.1 , 0.1) , dist);
return start
l_out:

<figure><img src="image_4.png" /><figcaption>Fig. 3. The WPTS of Pedestrian</figcaption></figure>

of our algorithm, we consider the transition-based semantics [9, 12] to our language and treat
each probabilistic program as a weighted probabilistic transition system (WPTS). A WPTS extends a
PTS [9, 12] with weights and an initial probability distribution.

*Definition 2.2 (WPTS).* A weighted probabilistic transition system (WPTS) Π is a tuple

$$
\Pi = (V_p, V_r, L, \ell_{\text{init}}, \ell_{\text{out}}, \mu_{\text{init}}, \mathcal{D}, \mathfrak{T}) \quad (\dagger)
$$

for which:

* $V_p$ and $V_r$ are finite disjoint sets of *program* and resp. *sampling* variables.

* L is a finite set of *locations* with special locations $\ell_{\text{init}}$, $\ell_{\text{out}} \in L$. Informally, $\ell_{\text{init}}$ is the initial location and $\ell_{\text{out}}$ represents program termination.

* $\mu_{\text{init}}$ is the initial probability distribution over $\mathbb{R}^{V_p}$ with a finite support (denoted by supp $(\mu_{\text{init}})$), while $\mathcal{D}$ is a function that assigns a probability distribution $\mathcal{D}(r)$ to each $r \in V_r$. We call each $v \in \text{supp}(\mu_{\text{init}})$ an initial program valuation, and abuse the notation so that $\mathcal{D}$ also denotes the independent joint distribution of all $\mathcal{D}(r)'$s ($r \in V_r$).

* $\mathfrak{T}$ is a finite set of *transitions* where each transition $\tau \in \mathfrak{T}$ is a tuple $\langle \ell, \phi, F_1, \dots, F_k \rangle$ such that (i) $\ell \in L$ is the *source location* of the transition, (ii) $\phi$ is the *guard condition* which is a predicate over variables $V_p$, and (iii) each $F_j := \langle \ell'_j, p_j, \text{upd}_j, wt_j \rangle$ is called a *weighted fork* for which (a) $\ell'_j \in L$ is the *destination location* of the fork, (b) $p_j \in (0, 1]$ is the probability of this fork, (c) $\text{upd}_j : \mathbb{R}^{|V_p|} \times \mathbb{R}^{|V_r|} \rightarrow \mathbb{R}^{|V_p|}$ is an *update function* that takes as inputs the current program and sampling valuations and returns an updated program valuation in the next step, and (d) $wt_j : \mathbb{R}^{|V_p|} \times \mathbb{R}^{|V_r|} \rightarrow [0, \infty)$ is a *score function* that gives the likelihood weight of this fork depending on the current program and sampling valuations.

In a WPTS, we use update and score functions to model the update on the program variables and
resp. the likelihood weight when running a basic block of statements in a program, respectively. If
there is no score statement in the block, then the score function is constantly 1. We always assume
that a WPTS Π is *deterministic* and *total*, i.e., (i) there is no program valuation that simultaneously
satisfies the guard conditions of two distinct transitions from the same source location, and (ii) the
disjunction of the guard conditions of all the transitions from any source location is a tautology.
The transformation from a probabilistic program into its WPTS can be done in a straightforward way (see e.g. [9, 13]).

**Example 2.3.** Fig. 3 shows the WPTS of the program in Fig. 2 which has two locations $\ell_{init}$, $\ell_{out}$. The circle nodes represent locations and square nodes model the forking behavior of transitions. An edge entering a square node is labeled with the condition of its respective transition, while an edge entering a circle node stands for a fork, which is associated with its probability, update functions and score functions that marked by w.<fnref n="3" /> The value of step is initialised to 0. An the initial probability distribution $\mu_{init}$ is determined by the joint distribution of (start, pos, dis, step) where start ~ uniform(0,3) and pos, dis, step observe the Dirac measures Dirac({start}), Dirac({0}) and Dirac({0}), respectively, e.g., the probability of the event "dis $\in$ {0}" equals 1. As step simply receives values from a sampling variable, we neglect it in the WPTS. $\square$

We say that a WPTS is *non-score-recursive* if for all transitions $\tau = (\ell, \phi, F_1, F_2, \dots, F_k)$ in the WPTS with each fork $F_j = (\ell'_j, p_j, upd_j, wt_j) (1 \le j \le k)$, we have that each score function $wt_j$ is constantly 1 (i.e., the multiplicative weight does not change) for every $\ell'_j \ne \ell_{out}$. Otherwise, the WPTS is *score-recursive*. Informally, a non-score-recursive WPTS has non-trivial score functions only on the transitions to the termination of a program, while a score-recursive WPTS has score statements in the execution of the program. For example, the WPTS of the program in Section 3.1 is non-score-recursive as the nontrivial (i.e., score values not equal to 1) score statement only appears to the termination, while the WPTS of the program in Section 3.2 is score recursive since it has score statements inside the loop body. In the case of a non-score-recursive WPTS, we say that the WPTS is *score-bounded* by a positive real $M > 0$ if for every $\tau = (\ell, \phi, F_1, F_2, \dots, F_k)$ in the WPTS with $F_j = (\ell'_j, p_j, upd_j, wt_j) (1 \le j \le k)$, we have that $|wt_j| \le M$ whenever $\ell'_j = \ell_{out}$.

Given a program valuation v and a predicate $\phi$ over variables $V_p$, we say that v satisfies $\phi$ (written as $v \models \phi$) if $\phi$ holds when the variables in $\phi$ are substituted by their values in v. A state is a pair $\Xi = (\ell, v)$ where $\ell \in L$ (resp. $v \in \mathbb{R}^{|V_p|}$) represents the current location (resp. program valuation), respectively, while a weighted state is a triple $\Theta = (\ell, v, w)$ where $(\ell, v)$ is a state and $w \in [0, \infty)$ represents the multiplicative likelihood weight accumulated so far.

Below we specify the semantics of a WPTS. Consider a WPTS $\Pi$ in the form of ($\dagger$). The semantics of $\Pi$ is formalized by the infinite sequence $\Gamma = \{\hat{\Theta}_n = (\hat{\ell}_n, \hat{v}_n, \hat{w}_n)\}_{n \ge 0}$ where each $(\hat{\ell}_n, \hat{v}_n, \hat{w}_n)$ is the random weighted state at the $n$th execution step of the WPTS such that $\hat{\ell}_n$ (resp. $\hat{v}_n, \hat{w}_n$) is the random variable for the location (resp. the random vector for the program valuation, the random variable for the multiplicative likelihood weight) at the $n$th step, respectively. The sequence $\Gamma$ starts with the initial random weighted state $\hat{\Theta}_0 = (\hat{\ell}_0, \hat{v}_0, \hat{w}_0)$ such that $\hat{\ell}_0$ is constantly $\ell_{init}$, $\hat{v}_0 \in \text{supp}(\mu_{init})$ is sampled from the initial distribution $\mu_{init}$ and the initial weight $\hat{w}_0$ is constantly set to $1^4$. Then, given the current random weighted state $\hat{\Theta}_n = (\hat{\ell}_n, \hat{v}_n, \hat{w}_n)$ at the $n$th step, the next random weighted state $\hat{\Theta}_{n+1} = (\hat{\ell}_{n+1}, \hat{v}_{n+1}, \hat{w}_{n+1})$ is determined by: (a) If $\hat{\ell}_n = \ell_{out}$, then $(\hat{\ell}_{n+1}, \hat{v}_{n+1}, \hat{w}_{n+1})$ takes the same weighted state as $(\hat{\ell}_n, \hat{v}_n, \hat{w}_n)$ (i.e., the next weighted state stays at the termination location $\ell_{out}$); (b) Otherwise, $\hat{\Theta}_{n+1}$ is determined by the following procedure:

* First, since the WPTS $\Pi$ is deterministic and total, we take the unique transition $\tau = (\ell_n, \phi, F_1, \dots, F_k)$ such that $\hat{v}_n = \phi$.

* Second, we choose a fork $F_j = (\ell'_j, p_j, upd_j, wt_j)$ with probability $p_j$.

* Third, we obtain a sampling valuation $r \in \text{supp}(D)$ by sampling each $r \in V_r$ independently w.r.t the probability distribution $D(r)$.

<footnote n="3">³Here we omit the update functions if the values of program variables are unchanged.</footnote>

<footnote n="4">⁴This follows the traditional setting in e.g. [4].</footnote>
* Finally, the value of the next random weighted state ($\hat{l}_{n+1}$, $\hat{v}_{n+1}$, $\hat{w}_{n+1}$) is determined as that of ($l'_j$, $upd_j(\hat{v}_n, r)$, $\hat{w}_n \cdot wt_j(\hat{v}_n, r)$).

Given the semantics, a program run of the WPTS Π is a concrete instance of Γ, i.e., an infinite sequence $ω = \{Θ_n\}_{n≥0}$ of weighted states where each $Θ_n = (l_n, v_n, w_n)$ is the concrete weighted state at the nth step in this program run with location $l_n$, program valuation $v_n$ and multiplicative likelihood weight $w_n$. A state $(ℓ, v)$ is called reachable if there exists a program run $ω = \{Θ_n\}_{n≥0}$ such that $Θ_n = (ℓ, v, w_n)$ for some $n$.

**Example 2.4.** Consider the WPTS in Example 2.3. Consider an initial program valuation $(1, 1, 0)$ which means that the initial values of start, pos, dis are $1, 1, 0$, respectively. Then starting from the initial weighted state $(l_{init}, (1, 1, 0), 1)$, a program run w.r.t the WPTS semantics above could be
$$ (l_{init}, (1, 1, 0), 1) \to (l_{init}, (1, 0.5, 0.5), 1) \to (l_{init}, (1, -0.1, 1.1), 1) \to (l_{out}, (1, -0.1, 1.1), 3.9894). \quad \square $$

Given an initial program valuation $v_{init}$ of a WPTS, one could construct a probability space over the program runs by their probabilistic evolution described above and standard constructions such as general state space Markov chains [32]. We denote the probability measure in the probability space by $P_{v_{init}}(-)$ and the expectation operator by $E_{v_{init}}[-]$.

## 2.4 Normalised Posterior Distribution

Before presenting the central problem of Bayesian probabilistic programming, i.e., analyzing normalised posterior distribution with our WPTS models, we introduce some technical concepts.

**Definition 2.5 (Termination).** The termination time of a WPTS Π is the random variable $T$ given by $T(ω) := \min\{n \in N \mid l_n = l_{out}\}$ for every program run $ω = (\{l_n, v_n, w_n\})_{n≥0}$ where $\min\emptyset := \infty$. That is, $T(ω)$ is the number of steps a program run $ω$ takes to reach the termination location $l_{out}$. A WPTS Π is almost-surely terminating (AST) if $P_{v_{init}}(T < \infty) = 1$ for all initial program valuations $v_{init} \in \text{supp}(\mu_{init})$.

**Definition 2.6 (Expected Weights).** Given a WPTS Π in the form of ($\dagger$), a designated initial program valuation $v_{init}$ and a measurable subset $U \in \Sigma_{\mathbb{R}|V_p|}$, the expected weight $[\Pi]_{v_{init}}(U)$ is defined as $[\Pi]_{v_{init}}(U) := E_{v_{init}}[[\hat{v}_T \in U] \cdot \hat{w}_T]$.

By definition, we have that $\hat{v}_T$ (resp. $\hat{w}_T$) is the random vector (resp. variable) of the program valuation (resp. the multiplicative likelihood weight) at termination, respectively. Thus, $[\Pi]_{v_{init}}(U)$ is the expectation of $\hat{w}_T$ that start from the state $(l_{init}, v_{init}, 1)$ and end with $\hat{v}_T \in U$. If $U = \mathbb{R}^{|V_p|}$, the restriction of $\hat{v}_T \in U$ can be removed.

Below we define the normalised posterior distribution (NPD) problem.

**Definition 2.7 (Normalised Posterior Distribution).** Given a WPTS Π in the form of ($\dagger$), the normalised posterior distribution (NPD) posterior$_{Π}$ of Π is defined by:

$$ \text{posterior}_{\Pi}(\mathcal{U}) := [\Pi](\mathcal{U})/Z_{\Pi} \text{ for all measurable subsets } \mathcal{U} \in \Sigma_{\mathbb{R}|V_p|}, $$

where $[\Pi](\mathcal{U}) := \int_V [\Pi]_V(\mathcal{U}) \cdot \mu_{init}(dv)$ is the unnormalised posterior distribution w.r.t. $\mathcal{U}, \mathcal{V} := \text{supp}(\mu_{init})$, and $Z_{\Pi} := [\Pi](\mathbb{R}^{|V_p|})$ is the normalising constant. The WPTS Π is called integrable if we have $0 < Z_{\Pi} < \infty$.

**Interval Bounds for NPD.** In this work, we consider the automated interval-bound analysis for NPD of a WPTS. Formally, we aim to derive an interval $[l, u] \subseteq [0, \infty)$ for an integrable WPTS Π and any measurable sets $U \in \Sigma_{\mathbb{R}|V_p|}$ as tight as possible such that $l \le \text{posterior}_{\Pi}(U) \le u$.

To achieve interval bounds for NPD, below we introduce the construction of a new WPTS Π$_U$ based on the original WPTS Π and a measurable set $U \in \Sigma_{\mathbb{R}|V_p|}$.
<figure><img src="image_5.png" /><figcaption>Fig. 4. An overview of our method</figcaption></figure>

Construction of $\Pi_{q_t}$. Consider a probabilistic program $P$ and its WPTS $\Pi$, given a measurable set $\mathcal{U} \in \Sigma_{\mathbb{R}|V_p|}$, we construct a new program $P_{q_t}$ by adding a conditional branch of the form “if $v_T \notin \mathcal{U}$ then score(0) fi” immediately after the termination of $P$ and obtain the WPTS $\Pi_{q_t}$ of $P_{q_t}$. Therefore, $\Pi$ and $\Pi_{q_t}$ have the same initial probability distribution $\mu_{init}$ and the same finite support $\mathcal{V} = \text{supp}(\mu_{init})$. The following proposition shows that interval-bound analysis for NPD can be reduced to interval-bound analysis for expected weights in the form $[\Pi]_v(\mathbb{R}^{|V_p|})$.

**PROPOSITION 2.8.** Given a WPTS $\Pi$ in the form of $(\dagger)$, a measurable set $\mathcal{U} \in \Sigma_{\mathbb{R}|V_p|}$ and the WPTS $\Pi_{q_t}$ constructed as above, we have that $[\Pi]_v(\mathcal{U}) = [\Pi_{q_t}]_v(\mathbb{R}^{|V_p|})$ for any $v \in \mathcal{V} = \text{supp}(\mu_{init})$. Furthermore, if there exist intervals $[l_1, u_1]$, $[l_2, u_2] \subseteq [0, \infty)$ such that $[\Pi_{q_t}]_v(\mathbb{R}^{|V_p|}) \in [l_1, u_1]$ and $[\Pi]_v(\mathbb{R}^{|V_p|}) \in [l_2, u_2]$ for any $v \in \mathcal{V}$, then we have two intervals $[l_q, u_q]$, $[l_Z, u_Z] \subseteq [0, \infty)$ such that the unnormalised posterior distribution $[\Pi](\mathcal{U}) = [l_q, u_q]$ and the normalising constant $Z_{\Pi} \in [l_Z, u_Z]$. Moreover, if $\Pi$ is integrable, i.e., $[l_Z, u_Z] \subseteq (0, \infty)$, then we can obtain the NPD posterior$_{\Pi}(\mathcal{U}) = [\frac{l_q}{u_Z}, \frac{u_q}{l_Z}]$.<fnref n="5" /> Note that by Definition 2.7, $l_q = \int_V l_1 \cdot \mu_{init}(dv)$, $u_q = \int_V u_1 \cdot \mu_{init}(dv)$, $l_Z = \int_V l_2 \cdot \mu_{init}(dv)$ and $u_Z = \int_V u_1 \cdot \mu_{init}(dv)$.

The proof of Proposition 2.8 is relegated to Appendix A.5. In the following, we will develop approaches to obtain interval bounds for expected weights.

## 3 MOTIVATING EXAMPLES

In this section, we present an overview of our novelties via two motivating examples, namely the Pedestrian example and the Phylogenetic model. The Pedestrian example is non-score-recursive and handled by our fixed-point approach, while the Phylogenetic model is score-recursive and addressed by our OST approach.

Before we go into the details of the examples, we first present the workflow of our approaches. Given a probabilistic program $P$ and a measurable set $\mathcal{U} \in \Sigma_{\mathbb{R}|V_p|}$, the workflow of our approaches is shown in Fig. 4: First, our parser transforms the input program $P$ into its WPTS $\Pi$ in the form of $(\dagger)$; Second, our approaches perform a truncation operation (to be introduced in Section 4.3) that restricts the WPTS $\Pi$ into the bounded range of program values so that whenever the program runs out of the range, our approaches directly have upper and lower bounds to over- and under-approximate the behaviour of the program outside the bounded range, the purpose of which is that a bounded range allows our approach to derive tight bounds via templates; Third, the truncated WPTS is further handled by our fixed-point or OST approach, depending on whether or not $\Pi$ is score-recursive, to generate polynomial bounds for expected weights at the termination of the program; Fourth, the interval bounds for the NPD posterior$_{\Pi}(\mathcal{U})$ are calculated by the interval-bound analysis for expected weights. An illustration of our workflow is given in Fig. 4.

<footnote n="5">⁵The interval bounds derived in this manner may be loose, but they are definitely correct.</footnote>
## 3.1 Pedestrian Random Walk

Consider the pedestrian random walk example [29] in Fig. 2. In this example, a pedestrian is lost on the way home, and she only knows that she is away from her house at most 3 km. Thus, she starts to repeatedly walk a uniformly random distance of at most 1 km in either direction of the road, until reaching her house. Upon she arrives, an odometer tells that she has walked 1.1 km totally. However, this odometer was once broken and the measured distance is normally distributed around the true distance with a standard deviation of 0.1 km. The question of this example is: what is the posterior distribution of the starting point? This example is modeled as a non-parametric probabilistic program whose number of loop iterations is unbounded. In the program, the variable *start* represents the starting point, *pos* records the current position of the pedestrian, *step* records the distance she walks in the next step, and *dis* records the total distance the pedestrian travelled so far. The probabilistic branch in the loop body specifies that the pedestrian walks either forward or backward, both with probability 0.5.

This program is non-score-recursive and was previously handled in [4] by exhaustive recursion unrolling that has the path-explosion problem. To circumvent the path-explosion problem, in this work we propose a fixed-point theorem to establish constraints for this example, and solve the constraints by polynomial templates. To utilize the fact that polynomial approximation is usually accurate over a bounded range, during the solving of the polynomial templates, we restrict the behaviour of the program within a bounded range (e.g., {(*pos*, *dis*) | *pos*, *dis* ∈ [0, 5]} ) and over-approximate the expected weights outside the bounded range by an interval (e.g. [0, 2.1 × 10⁻³³]). Note that here we omit the program variable *start* as its value is determined once and has nothing to do with loop iterations. The solving of the polynomial template follows the previous approaches [9, 10, 13] that consider linear/semidefinite programming. Our approach derives comparable bounds to the approach in Beutner et al. [4] while our runtime is two-thirds of that of Beutner et al. [4].

## 3.2 Phylogenetic Birth Model

Consider a simplified version of the phylogenetic birth model [38] in Fig. 5, where a species arises with a birth-rate *lambda*, and it propagates with a constant likelihood of 1.1 at some time interval.<fnref n="6" /> This example is modelled as a probabilistic loop, where the variables *lambda*, *time*, *amount*, *wait* stand for the birth rate of the species, the remaining propagation time, the current amount of the species and the propagation time to be spent, respectively. The variable *lambda* is associated with a prior distribution, and the NPD problem is to infer its posterior distribution given the species evolution described by the loop. The WPTS of the program is given in Fig. 6.

This program cannot be handled by previous approaches (such as [4, 20]) since it is a score-recursive program with an unbounded while loop and its score weight is greater than 1. To see why such a score-recursive program is nontrivial to tackle, consider a simple loop where in each loop iteration, the loop terminates directly with probability 1/2, and continues to the next loop iteration with the same probability. At the end of each loop iteration, a score command “*score(3)*” is executed. It follows that the normalising constant in NPD is equal to $\sum_{n=1}^{\infty} \mathbb{P}(T = n) \cdot 3^n = \sum_{n=1}^{\infty} (\frac{3}{2})^n = \infty$, so that the infinity makes the posterior distribution invalid. One can observe that in this example the main problem lies at the fact that the scaling speed of the likelihood weight (i.e., 3) is higher than that for program termination (i.e., 1/2).

To tackle the difficulty of having scores greater than 1 in a loop iteration, we propose a novel multiplicative variant of Optional Stopping Theorem (OST) that addresses score statements. Based on our OST variant, we apply truncation and polynomial template solving as in the non-score-recursive case. For example, we can restrict the behaviour of the program within a bounded range

<footnote n="6">⁶For simplicity, we assume constant weights that can be viewed as over-approximation for a continuous density function.</footnote>
such as {($\lambda\text{Lambda}$, $\text{time}$) | $\lambda\text{Lambda} \in [0, 3]$ and $\text{time} \in [0, 10]$}, and over-approximate the expected weights outside the bounded range by an interval bound of polynomial functions derived from the same OST variant and polynomial-template method but without truncation. For the same reason as in Section 3.1, here we can safely omit the program variable *amount*. Our experimental result on this example shows that the derived bounds match the simulation result with $10^6$ samples.

```
lambda := sample uniform(0, 2);
time := 10; amount := 0;

l_init: while time ≥ 0 do
    wait := sample uniform(0, 0.5);
    time := time - wait;
    if prob (0.5 · lambda) then
        birth := sample uniform (0, 0.01);
        amount := amount + birth;
        score (1.1)
    fi;
od;
return lambda
l_out:
```

<figure><img src="image_6.png" /><figcaption>Fig. 6. The WPTS of Phylogenetic Model</figcaption></figure>

# 4 THEORETICAL APPROACHES

In this section, we present our theoretical approaches for interval-bound analysis of expected weights in the form $\llbracket\Pi\rrbracket_v(\mathbb{R}^{|V_P|})$, namely the fixed-point approach, the OST approach. To obtain tight results, we also introduce the truncation operation.

## 4.1 The Fixed-Point Approach

We review some basic concepts of lattice theory. Given a partial order $\sqsubseteq$ on a set $K$ and a subset $K' \subseteq K$, an *upper bound* of $K'$ is an element $u \in K$ that is no smaller than every element of $K'$, i.e., $\forall k' \in K'$. $k' \sqsubseteq u$. Similarly, a *lower bound* for $K'$ is an element $l$ that is no greater than every element of $K'$, i.e. $\forall k' \in K'$. $l \sqsubseteq k'$. The *supremum* of $K'$, denoted by $\biguplus K'$, is an element $u^* \in K$ such that $u^*$ is an upper-bound of $K'$ and for every upper bound $u$ of $K'$, we have $u^* \sqsubseteq u$. Similarly, the *infimum* $\bigwedge K'$ is a lower bound $l^*$ of $K'$ such that for every lower-bound $l$ of $K'$, we have $l \sqsubseteq l^*$. We define $\perp := \bigwedge K$ and $\uparrow := \bigwedge K$. In general, suprema and infima may not exist.

A partially-ordered set $(K, \sqsubseteq)$ is called a *complete lattice* if every subset $K' \subseteq K$ has a supremum and an infimum. Given a partial order $(K, \sqsubseteq)$, a function $f: K \to K$ is called *monotone* if for every $k_1 \sqsubseteq k_2$ in $K$, we have $f(k_1) \sqsubseteq f(k_2)$.

Given a complete lattice $(K, \sqsubseteq)$, a function $f: K \to K$ is called *continuous* if for every increasing chain $k_0 \sqsubseteq k_1 \sqsubseteq \dots$ in $K$, we have $f(\sqcup\{k_n\}_{n=0}^\infty) = \sqcup\{f(k_n)\}_{n=0}^\infty$, and *cocontinuous* if for every decreasing chain $k_0 \supseteq k_1 \supseteq \dots$ of elements of $K$, we have $f(\sqcup\{k_n\}_{n=0}^\infty) = \sqcup\{f(k_n)\}_{n=0}^\infty$. An element $k \in K$ is called a *fixed-point* if $f(k) = k$. Moreover, $k$ is a *pre fixed-point* if $f(k) \sqsubseteq k$ and a *post fixed-point* if $k \sqsubseteq f(k)$. The *least fixed-point* of $f$, denoted by $\text{lfp }f$, is the fixed-point that is no greater than every fixed-point under $\sqsubseteq$. Analogously, the *greatest fixed-point* of $f$, denoted by $\text{gfp }f$, is the fixed-point that is no smaller than all fixed-points.
**THEOREM 4.1 (KLEENE [42]).** Let $(K, \sqsubseteq)$ be a complete lattice and $f: K \to K$ be an continuous function. Then, we have

$$ \text{lfp } f = \sqcup_{i \ge 0} \{f^{(i)}(\perp)\}. $$

Analogously, if $f$ is cocontinuous, then we have

$$ \text{gfp } f = \sqcap_{i \ge 0} \{f^{(i)}(\top)\}. $$

In this work, we apply the well-known Tarski's fixed-point theorem.

**THEOREM 4.2 (TARSKI [46]).** Let $(K, \sqsubseteq)$ be a complete lattice and $f: K \to K$ a monotone function. Then, both lfp $f$ and gfp $f$ exist. Moreover, we have

$$ \text{lfp } f = \sqcap \{x \mid f(x) \sqsubseteq x\} \text{ and gfp } f = \sqcup \{x \mid x \sqsubseteq f(x)\}. $$

Based on Theorem 4.2, we present our fixed-point approach for non-score-recursive WPTS's. Below we fix a non-score-recursive WPTS $\Pi$ in the form of $(+)$. Given a maximum finite value $M \in [0, \infty)$, we define a state function as a function $h: \Lambda \to [-M, M]$ such that for all $v \in \mathbb{R}^{|V_p|}$, we have that $h(\ell_{out}, v) \in [0, M]$. We denote the set of all state functions with maximum value $M$ by $\mathcal{K}_M$. We also use the usual partial order $\le$ on $\mathcal{K}_M$ that is defined in the pointwise fashion, i.e., for any $h_1, h_2 \in \mathcal{K}_M$, $h_1 \le h_2$ iff $h_1(\Xi) \le h_2(\Xi)$ for all $\Xi \in \Lambda$. It is straightforward to verify that $(\mathcal{K}_M, \le)$ is a complete lattice. For example, the top element $\top$ (resp. the bottom element $\bot$) in the set $\mathcal{K}_M$ is the constant function that maps every state $\Xi$ to $M$ (resp. $-M$), and the least upper bounds (resp. greatest lower bound) can be given via the pointwise infimum (resp. supremum), respectively.

To connect the complete lattice $(\mathcal{K}_M, \le)$ of state functions with expected weights (Definition 2.6), we define a special state function, called *expected-weight function* $ew_{\Pi}$, by $ew_{\Pi}(\ell_{init}, v) := [\Pi]_v(\mathbb{R}^{|V_p|})$, and omit the subscript $\Pi$ if it is clear from the context. Informally, $ew_{\Pi}(\ell_{init}, v)$ is the expected weight of all program runs that start from $v$ without the restriction of the subset $\mathcal{U}$. In this work, we consider the following monotone function over the complete lattice $(\mathcal{K}_M, \le)$.

**Definition 4.3 (Expected-Weight Transformer).** Given a finite maximum value $M \in [1, \infty)$, the *expected-weight transformer* $ewt_{\Pi}: \mathcal{K}_M \to \mathcal{K}_M$ is the higher-order function such that for each state function $h \in \mathcal{K}_M$ and state $(\ell, v)$, if $\tau = (\ell, \phi, F_1, \dots, F_k)$ is the unique transition that satisfies $v \models \phi$ and $F_j = (\ell'_j, p_j, upd_j, wt_j)$ for each $1 \le j \le k$, then we have that

$$ ewt_{\Pi}(h)(\ell, v) := \begin{cases} \sum_{j=1}^{k} p_j \cdot \mathbb{E}_r \left[ wt_j(v, r) \cdot h(\ell'_j, upd_j(v, r)) \right] & \text{if } \ell \neq \ell_{out} \\ 1 & \text{otherwise} \end{cases}. \quad (1) $$

Here the expectation $\mathbb{E}_r[-]$ is taken over a sampling valuation $r$ that observes the independent joint probability distributions of all the sampling variables $r \in V_r$.

Informally, given a state function $h$, the *expected-weight transformer* $ewt_{\Pi}$ computes the expected weight $ewt_{\Pi}(h)$ after one step of the WPTS execution. From the monotonicity of the Lebesgue integral (from the definition of expectation), we have that $ewt_{\Pi}$ is monotone. Note that although the Lebesgue integral of a function $f$ usually requires measurability of the function, its preliminary definition is given via a supremum over simple functions not exceeding the function, and hence can be applied to all functions (including non-measurable functions) for which the monotonicity remains to be true. Therefore, we do not need to impose measurability here. In the following, we will omit the subscript $\Pi$ in $ewt_{\Pi}$ if it is clear from the context.

**Definition 4.4 (Potential Weight Functions).** A potential upper weight function (PUWF) is a function $h: L \times \text{Val}_{v_p} \to \mathbb{R}$ that has the properties below:

(C1) for all reachable states $(\ell, v)$ with $\ell \neq \ell_{out}$, we have $ewt(h)(\ell, v) \le h(\ell, v)$;
(C2) for all reachable states $(\ell, \mathbf{v})$ such that $\ell = \ell_{\text{out}}$, we have $h(\ell, \mathbf{v}) = 1$.

Analogously, a potential lower weight function (PLWF) is a function $h: L \times \text{Val}_{V_p} \to \mathbb{R}$ that satisfies the conditions (C1') and (C2), for which the condition (C1') is almost the same as (C1) except for that "$ewt(h)(\ell, \mathbf{v}) \le h(\ell, \mathbf{v})$" is replaced with "$ewt(h)(\ell, \mathbf{v}) \ge h(\ell, \mathbf{v})$".

Informally, a PUWF is a state function that satisfies the pre fixed-point condition of $ewt$ at non-terminating locations, and equals one at the termination location. A PLWF is defined similarly, with the difference that we use the post fixed-point condition instead.

The following result lays the backbone of our fixed-point approach.

**THEOREM 4.5.** Let $\Pi$ be a non-score-recursive WPTS that is score-bounded by a positive real $M > 0$. If the WPTS $\Pi$ has the AST property, then the expected-weight function $ew$ is the unique fixed-point of the higher-order function $ewt$ on the complete lattice ($\mathcal{K}_{\max\{M, 1\}}$, $\le$).

**PROOF SKETCH.** Let $M' := \max\{M, 1\}$. We first prove that the expected-weight function $ew$ is the least fixed-point of the higher-order function $ewt$. Then we prove that the expected-weight transformer $ewt: \mathcal{K}' \to \mathcal{K}'$ is both continuous and cocontinuous. Finally, according to Theorem B.1, we prove the uniqueness that $\text{lfp } ewt(\ell, \mathbf{v}) = \text{gfp } ewt(\ell, \mathbf{v})$, i.e., the fixed-point $ew$ is unique. For the space limitation, we relegate the detailed proof for Theorem 4.5 to Appendix B.2. $\square$

By a combination of Theorem 4.5 and Theorem 4.2, one has that it suffices to derive a pre fixed-point of $ewt$ to obtain an upper bound for $ew$, and a post fixed-point to obtain a lower bound in the case of almost-sure termination of the program.

**THEOREM 4.6 (FIXED-POINT APPROACH).** Let $\Pi$ be a non-score-recursive WPTS that is score-bounded by a finite value $M > 0$. If the WPTS $\Pi$ has the AST property, then $[\Pi]_{v_{\text{init}}}(\mathbb{R}^{|V_p|}) \le h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$ (resp. $[\Pi]_{v_{\text{init}}}(\mathbb{R}^{|V_p|}) \ge h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$) for any PUWF (resp. PLWF) $h$ over $\Pi$ and initial state $(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$, respectively.

**REMARK 1.** Notice that we always require a finite maximal value $M$ in our fixed-point approach. In general, our fixed-point approach cannot be applied to score-recursive Bayesian probabilistic programs, as in such a program a finite maximum weight $M$ may not exist.

## 4.2 The OST Approach

As stated in Remark 1, our fixed-point approach cannot handle score-recursive Bayesian probabilistic programs. To tackle score-recursive programs, we consider the adaption of Optional Stopping Theorem (OST) to our case. OST is a classical theorem in martingale theory that characterizes the relationship between the expected values initially and at a stopping time in a supermartingale. Below we first present the classical form of OST.

**THEOREM 4.7 (OPTIONAL STOPPING THEOREM (OST) [54]).** Let $\{X_n\}_{n=0}^{\infty}$ be a supermartingale adapted to a filtration $\mathcal{F} = \{\mathcal{F}_n\}_{n=0}^{\infty}$, and $k$ be a stopping time w.r.t. the filtration $\mathcal{F}$. Then the following condition is sufficient to ensure that $\mathbb{E}(|X_\kappa|) < \infty$ and $\mathbb{E}(X_\kappa) \le \mathbb{E}(X_0)$:

* $\mathbb{E}(\kappa) < \infty$, and

* (bounded difference) there exists a constant $C > 0$ such that for all $n \ge 0$, $|X_{n+1} - X_n| \le C$ holds almost surely.

In our NPD problem, as the score statements accumulate weights in a multiplicative fashion, the classical OST cannot be applied since the bounded difference condition may be violated. To address this difficulty, we propose a novel variant of OST that tackles the multiplicative feature from score statements.
**THEOREM 4.8 (OST VARIANT).** Let {$X_n$}$_{n=0}^\infty$ be a supermartingale adapted to a filtration $\mathcal{F} =$ {$\mathcal{F}_n$}$_{n=0}^\infty$, and $\kappa$ be a stopping time w.r.t. the filtration $\mathcal{F}$. Then the following condition $(\mathcal{U})$ is sufficient to ensure that $\mathbb{E}(|X_\kappa|) < \infty$ and $\mathbb{E}(X_\kappa) \le \mathbb{E}(X_0)$:

(i) There exist integers $b_1, b_2 > 0$ and real numbers $c_1 > 0, c_2 > c_3 > 0$ such that (i) $\mathbb{P}(\kappa > n) \le c_1 \cdot e^{-c_2 \cdot n}$ for sufficiently large $n \in \mathbb{N}$, and (ii) for all $n \in \mathbb{N}$, $|X_{n+1} - X_n| \le b_1 \cdot n^{b_2} \cdot e^{c_3 \cdot n}$ holds almost surely.

Our OST variant extends the classical OST with the relaxation that we allow the magnitude of the next random variable $X_{n+1}$ to be bounded by that of $X_n$ with a multiplicative factor $e^{c_3}$, which corresponds to the multiplicative feature from score statements. The intuition of the theorem is that to cancel the effect of the multiplicative factor, we require in extra the exponential decrease in $\mathbb{P}(\kappa > n) \le c_1 \cdot e^{-c_2 \cdot n}$. The proof resembles [53, Theorem 5.2] and is relegated to Appendix B.3.

Below we show how our OST variant can be applied to handle score-recursive WPTS's. Fix a WPTS $\Pi$ in the form of (†). In the rest of this subsection, we reuse the expected-weight transformer ewt defined in Definition 4.3 and potential weight functions given in Definition 4.4. The slight difference is that for the expected-weight transformer in the context of a score-recursive WPTS, we consider that the weight function $wt_j$ before termination may not be constantly 1.

To apply our multiplicative OST variant, we impose a bounded update requirement as in [53]. A WPTS $\Pi$ has the *bounded update* property if there exists a constant $\kappa > 0$ such that for every reachable state $(\ell, \mathbf{v})$, if $\tau = (\ell, \phi, F_1, \dots, F_k)$ is the unique transition with each fork $F_j = (\ell'_j, p_j, \text{upd}_j, wt_j)$ such that $\mathbf{v} = |\phi$, then we have that $\forall \mathbf{r} \in \text{supp}(\mathcal{D}) \ \forall x \in V_p, |\text{upd}_j(\mathbf{v}, \mathbf{r})(x) - \mathbf{v}(x)| \le \kappa$.

We apply our OST variant (Theorem 4.8) in a way similar to [53, Theorem 6.10 and Theorem 6.12] to obtain the main theorem for our OST-based approach. To be more precise, we construct a stochastic process from a PUWF (or the negative of a PLWF) and show that the stochastic process is a supermartingale and the condition $(\mathcal{U})$ is fulfilled. The statement of the theorem is as follows.

**THEOREM 4.9 (OST APPROACH).** Let $\Pi$ be a score-recursive WPTS that has the bounded update property. Suppose that there exist real numbers $c_1 > 0$ and $c_2 > c_3 > 0$ such that

(E1) $\mathbb{P}(T > n) \le c_1 \cdot e^{-c_2 \cdot n}$ for sufficiently large $n \in \mathbb{N}$, and

(E2) for each score function wt in $\Pi$ we have $|wt| \le e^{c_3}$.

Then for any PUWF (resp. PLWF) $h$ over $\Pi$, we have that $\llbracket \Pi \rrbracket_{v_{\text{init}}} (\mathbb{R}^{|V_p|}) \le h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$ (resp. $\llbracket \Pi \rrbracket_{v_{\text{init}}} (\mathbb{R}^{|V_p|}) \ge h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$) for any initial state $(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$, respectively.

**PROOF SKETCH.** For the upper bounds, we define the stochastic process {$X_n$}$_{n=0}^\infty$ as $X_n := h(\ell_n, \mathbf{v}_n)$ where $(\ell_n, \mathbf{v}_n)$ is the program state at the $n$th step of a program run. Then we construct a stochastic process {$Y_n$}$_{n=0}^\infty$ such that $Y_n := X_n \cdot \prod_{i=0}^{n-1} W_i$ where $W_i$ is the weight at the $i$th step of the program run. We consider the termination time $T$ of $\Pi$ and prove that {$Y_n$}$_{n=0}^\infty$ satisfies the prerequisites of our OST variant (Theorem 4.8). This proof depends on the assumption that $\Pi$ has concentration and bounded update properties, and the score functions in $\Pi$ are also bounded (Item 2, Theorem 4.9)). Then by applying Theorem 4.8, we obtain that $\mathbb{E}[Y_T] \le \mathbb{E}[Y_0]$. By (C2) in Definition 4.4, we have that $Y_T = h(\ell_T, \mathbf{v}_T) \cdot \prod_{i=0}^{T-1} W_i = \hat{w}_T$. Thus, we have that $\llbracket \Pi \rrbracket_{v_{\text{init}}} (\mathbb{R}^{|V_p|}) = \mathbb{E}_{v_{\text{init}}} [\hat{w}_T] = \mathbb{E} [\prod_{i=0}^{T-1} W_i] \le \mathbb{E}[Y_0] = h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$. For the lower bounds, the proof is similar. The more detailed proof is relegated to Appendix B.4. $\square$

## 4.3 Truncation over WPTS's

In the following, we propose a truncation operation for a WPTS that restricts the value of every program variable in the WPTS to a prescribed bounded range. We consider that a bounded range for a program variable could be either $[-R, R]$ ($R > 0$), or $[0, R]$, $[-R, 0]$ if the value of the program variable is guaranteed to be non-negative or non-positive.
To present our truncation operation, we define the technical notions of truncation and approximation functions. A truncation function $\mathcal{B}$ is a function that maps every program variable $x \in V_p$ to a bounded interval $\mathcal{B}(x)$ in $\mathbb{R}$ that specifies the bounded range of the variable $x$. We denote by $\Phi_{\mathcal{B}}$ the formula $\bigwedge_{x \in V_p} x \in \mathcal{B}(x)$ for a truncation function $\mathcal{B}$. An approximation function is a function $\mathcal{M}: \mathbb{R}^{|V_p|} \to [0, \infty)$ such that each $\mathcal{M}(\mathbf{v})$ ($\mathbf{v} \in \mathbb{R}^{|V_p|}$) is intended to be an over- or under-approximation of the expected weight $\llbracket \Pi \rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_p|})$ outside the bounded range specified by $\Phi_{\mathcal{B}}$. The truncation operation is given by the following definition.

**Definition 4.10 (Truncation Operation).** Let $\Pi$ be a WPTS in the form of $(\dagger)$. Given a truncation function $\mathcal{B}$ and an approximation function $\mathcal{M}$, the truncated WPTS $\Pi_{\mathcal{B},\mathcal{M}}$ w.r.t. $\mathcal{B}$ and $\mathcal{M}$ is defined as $\Pi_{\mathcal{B},\mathcal{M}} := (V_p, V_r, L \cup \{\#\}, \ell_{\text{init}}, \ell_{\text{out}}, \mu_{\text{init}}, \mathcal{D}, \tilde{\mathcal{B}}_{\mathcal{B},\mathcal{M}})$ where $\#$ is a fresh deadlock location and the transition relation $\tilde{\mathcal{B}}_{\mathcal{B},\mathcal{M}}$ is given by

$$ \begin{aligned} \tilde{\mathcal{B}}_{\mathcal{B},\mathcal{M}} &:= \{ \langle \ell, \phi \wedge \Phi_{\mathcal{B}}, F_1, \dots, F_k \rangle \mid \langle \ell, \phi, F_1, \dots, F_k \rangle \in \tilde{\mathcal{B}} \text{ and } \ell \neq \ell_{\text{out}} \} \\ &\cup \{ \langle \ell, \phi \wedge (-\Phi_{\mathcal{B}}), F_1^{\mathcal{M},\#}, \dots, F_k^{\mathcal{M},\#} \rangle \mid \langle \ell, \phi, F_1, \dots, F_k \rangle \in \tilde{\mathcal{B}} \text{ and } \ell \neq \ell_{\text{out}} \} \\ &\cup \{ \langle \ell_{\text{out}}, \text{true}, F_{\ell_{\text{out}}}, \rangle, \langle \# , \text{true}, F_{\#} \rangle \} \end{aligned} \quad (\ddagger) $$

for which (a) we have $F_{\ell} := \langle \ell, 1, id, \bar{1} \rangle$ ($\ell \in \{\ell_{\text{out}}, \#\}$) where id is the identity function and $\bar{1}$ is the constant function that always takes the value 1, and (b) for a fork $F = \langle \ell', p, upd, wt \rangle$ in the original WPTS $\Pi$ we have $F^{\mathcal{M},\#} := F$ if $\ell' = \ell_{\text{out}}$ and $F^{\mathcal{M},\#} := \langle \# , p, upd, \mathcal{M} \rangle$ otherwise.

Thus, the truncated WPTS is obtained from the original one by first restraining each transition to the bounded range $\Phi_{\mathcal{B}}$ and then redirecting to the fresh deadlock location $\#\$ all the situations jumping out of the bounded range and not going to the termination location. To make the truncated WPTS deterministic and total, we add the self-loop $\langle \# , \text{true}, F_{\#} \rangle$. Our main theorem shows that by choosing an appropriate approximation function $\mathcal{M}$ in the truncation, one can obtain upper/lower approximation of the original WPTS.

**THEOREM 4.11.** Let $\Pi$ be a WPTS in the form of $(\dagger)$, $\mathcal{B}$ a truncation function and $\mathcal{M}$ an approximation function. Suppose that the following condition $(\ast)$ holds:

(\ast) for each fork $F^{\mathcal{M},\#} = \langle \# , p, upd, \mathcal{M} \rangle$ in the truncated WPTS $\Pi_{\mathcal{B},\mathcal{M}}$ that is derived from some fork $F = \langle \ell', p, upd, wt \rangle$ with the source location $\ell$ in the original WPTS (see sentence (b) in Definition 4.10), we have that $\llbracket \Pi \rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_p|}) \le \mathcal{M}(\mathbf{v})$ for all $\mathbf{v}$ such that the state $(\ell, \mathbf{v})$ is reachable and $\mathbf{v} \nleq \Phi_{\mathcal{B}}$.

Then $\llbracket \Pi \rrbracket_{\text{v}_{\text{init}}}(\mathbb{R}^{|V_p|}) \le \llbracket \Pi_{\mathcal{B},\mathcal{M}} \rrbracket_{\text{v}_{\text{init}}}(\mathbb{R}^{|V_p|})$ for all initial program valuations $\mathbf{v}_{\text{init}}$. Analogously, if it holds the condition $(\star)$ which is almost the same as $(\ast)$ except for that $\llbracket \Pi \rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_p|}) \le \mathcal{M}(\mathbf{v})$ is replaced with $\llbracket \Pi \rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_p|}) \ge \mathcal{M}(\mathbf{v})$, then we have $\llbracket \Pi \rrbracket_{\text{v}_{\text{init}}}(\mathbb{R}^{|V_p|}) \ge \llbracket \Pi_{\mathcal{B},\mathcal{M}} \rrbracket_{\text{v}_{\text{init}}}(\mathbb{R}^{|V_p|})$ for all initial program valuations $\mathbf{v}_{\text{init}}$.

The theorem above states that if the approximation function gives correct bounds for the expected weights of the original WPTS outside the bounded range, then the bounds for the expected weights of the truncated WPTS are also correct bounds for the expected weights of the original WPTS. The detailed proof is relegated to Appendix C.

**Example 4.12.** Recall the Pedestrian example in Fig. 2, here we make truncation to this example and generate its truncated WPTS. The truncation function $\mathcal{B}$ is defined such that $\mathcal{B}(\text{pos}) = [0, 5]$, $\mathcal{B}(\text{dis}) = [0, 5]$, so $\Phi_{\mathcal{B}} = 0 \le \text{pos} \le 5 \wedge 0 \le \text{dis} \le 5$. Following the truncation operation in $(\ddagger)$ from Definition 4.10, we can obtain $\Pi_{\mathcal{B},\mathcal{M}}$ as shown in Fig. 7. $\square$
<figure><img src="image_7.png" /><figcaption>Fig. 7. The Truncated WPTS of Pedestrian</figcaption></figure>

# 5 ALGORITHMIC APPROACHES

In this section, we present algorithmic implementation of our theoretical approaches in Section 4. We first have some assumptions on the input probabilistic program:

* To enable the exact calculation of the integrals over the probability density functions in a Bayesian probabilistic program possible, we assume that every probability density function in a sampling statement is a polynomial. Likewise, in score-recursive programs, we require that the probability density function in any score statement is polynomial. Our approach can handle non-polynomial density functions in sampling and score statements by having their polynomial approximations, for which we leave as a future work.

* To capture all possible program executions, we assume that the input probabilistic program is accompanied with invariants (see e.g. [8, 44]) to over-approximate reachable states. We follow [15, 44] to consider *affine invariants*. An affine invariant for a WPTS is a map $I$ that assigns to each location $\ell$ a (conjunctive) system $I(\ell)$ of affine inequalities such that for all reachable states $(\ell, v)$, every affine inequality in $I(\ell)$ holds w.r.t the program values given in the program valuation $v$.

We also present the high-level technical setting of our algorithms. Consider an input WPTS $\Pi$ in the form of (†). Recall the normalising constant $Z_\Pi$ in Definition 2.7. To tighten the interval bounds for $Z_\Pi$, we split the set $\mathcal{V} := \text{supp}(\mu_{\text{init}})$ of initial program valuations into $m \ge 1$ disjoint partitions $\mathcal{V}_1, \dots, \mathcal{V}_m$ such that $\mathcal{V} = \bigcup_{i=1}^m \mathcal{V}_i$ and $\mathcal{V}_i \cap \mathcal{V}_j = \emptyset$ for any $i \ne j \in [1, m]$.<fnref n="7" /> Our algorithms tackle each $\mathcal{V}_i$ ($1 \le i \le m$) separately to obtain polynomial interval bounds $[l_i, u_i]$ (where $l_i, u_i$ are polynomials) such that $l_i(v) \le [\Pi]_v(\mathbb{R}^{|V_P|}) \le u_i(v)$ for all $v \in \mathcal{V}_i$. It follows that the interval bound $[l_Z, u_Z]$ for $Z_\Pi$ can be derived by integrals of polynomial bounds over all $\mathcal{V}_i$'s, that is,

$$l_Z := \sum_{i=1}^{m} \int_{\mathcal{V}_i} l_i(v) \mu_{\text{init}}(dv) \le Z_\Pi = \sum_{i=1}^{m} \int_{\mathcal{V}_i} [\Pi]_v(\mathbb{R}^{|V_P|}) \cdot \mu_{\text{init}}(dv) \le \sum_{i=1}^{m} \int_{\mathcal{V}_i} u_i(v) \mu_{\text{init}}(dv) =: u_Z. \quad (\blacksquare)$$

Given a measurable set $\mathcal{U} \in \Sigma_{\mathbb{R}^{|V_P|}}$, the interval bounds $[l_\mathcal{U}, u_\mathcal{U}]$ for the unnormalised posterior distribution $[\Pi](\mathcal{U})$ (in Definition 2.7) can be obtained analogously by deriving polynomial bounds for $[\Pi u_\mathcal{U}]_v(\mathbb{R}^{|V_P|})$ (see Proposition 2.8).

<footnote n="7">⁷If $m=1$, the set $\mathcal{V}$ is not split.</footnote>
Below we present our algorithmic methods to deriving interval bounds for the normalising constant. The unnormalised posterior distribution can be derived in the same manner. Our algorithms are template-based and have the following stages:

**Stage 1: Input.** First, our algorithms receive a WPTS Π parsed from a Bayesian probabilistic program *P* written in our PPL and an affine invariant *I* for the WPTS as the basic input. Besides, the algorithms receive as auxiliary inputs a truncation function *B* and two approximation functions *M*<sub>up</sub>, *M*<sub>low</sub> that respectively fulfill the conditions (*), (★) in the statement of Theorem 4.11. The truncation function is represented by the bounded intervals for each program variable, and the approximation functions are derived either from the score function at the termination of a non-score-recursive program or by applying our OST variant directly to a score-recursive program (without truncation). We also have two parameters *d*, *m*, for which *d* is the degree of our polynomial template and *m* is the number of partitions for the set *V* of initial program valuations (refer to ★).
Moreover, if the program *P* is non-score-recursive and its score function *g* at the termination is non-polynomial, our algorithms take as extra inputs a (piecewise) polynomial approximation *g'* of *g* and an error bound *ε* > 0 such that |g(v) - g'(v)| ≤ *ε* for all v ∈ *B*' (*B*' is the extended bounded range of program variables that derived from a one-step program run and the truncation function *B*, which will be introduced in **Step A1** below).

**Stage 2: Partition and Truncation.** Next, our algorithms fetch the set $\mathcal{V}$ = supp $(\mu_{\text{init}})$ from the WPTS $\Pi$, partition it into $m$ disjoint subsets $\mathcal{V}_1, \dots, \mathcal{V}_m$ and construct a set $\mathcal{W} = \{\mathbf{v}_1, \dots, \mathbf{v}_m\}$ such that each $\mathbf{v}_i \in \mathcal{V}_i$. Our algorithms then perform the truncation operation to $\Pi$ w.r.t the input truncation function $\mathcal{B}$ and approximation functions $M_{\text{up}}, M_{\text{low}}$. To synthesize upper bounds for expected weights, we take the truncation w.r.t $\mathcal{B}$ and $M_{\text{up}}$ to generate a truncated WPTS $\Pi_{\mathcal{B}}, M_{\text{up}}$. For lower bounds, we generate the truncated WPTS $\Pi_{\mathcal{B}}, M_{\text{low}}$.

**Example 5.1.** Recall the Pedestrian example in Section 3.1 and its WPTS Π in Fig. 3. We derive an invariant *I* simply from the loop guard so that *I*(*l*<sub>init</sub>) = *pos* ≥ 0 and *I*(*l*<sub>out</sub>) = *pos* < 0. The truncation function *B* is defined such that *B*(*pos*) = [0, 5] and *B*(*dis*) = [0, 5]. We pick the constant approximation functions $M_{\text{up}} = 2.1 \times 10^{-330}$ and $M_{\text{low}} = 0$ to bound the expected weights beyond the truncated range, so that we obtain two truncated WPTSs $\Pi_{\mathcal{B}}, M_{\text{up}}$ and $\Pi_{\mathcal{B}}, M_{\text{low}}$ which are similar to that in Fig. 7. We also choose the algorithm parameters as $d = 1$ and $m = 30$.<fnref n="8" /> Since the program is non-score-recursive and its score function $g$ at the termination is non-polynomial (i.e., $g(\mathit{dis}) = \mathit{pdf}(\mathit{normal}(1.1, 0.1), \mathit{dis})$), we choose a polynomial approximation $g'$ of $g$ with the error bound $\epsilon = 10^{-5}$. We obtain the set $\mathcal{V} = \{(\mathit{pos}, \mathit{dis}) \mid \mathit{pos} \in [0, 3], \mathit{dis} = 0\}$. Since the value of $\mathit{dis}$ in $\mathcal{V}$ is fixed, we partition $\mathcal{V}$ uniformly into $m = 30$ disjoint subsets on the dimension $\mathit{pos}$, i.e., $\mathcal{V}_1 = \{(\mathit{pos}, \mathit{dis}) \mid \mathit{pos} \in [0, 0.1], \mathit{dis} = 0\}, \dots, \mathcal{V}_{30} = \{(\mathit{pos}, \mathit{dis}) \mid \mathit{pos} \in [2.9, 3], \mathit{dis} = 0\}$. We calculate the midpoints of the dimension $\mathit{pos}$ for all $\mathcal{V}_i$'s, and construct the set $\mathcal{W} = \{(0.05, 0), (0.15, 0), \dots, (2.95, 0)\}$. □

**Stage 3: Template Solving.** Then our algorithms establish *d*-degree polynomial templates for $\Pi_{\mathcal{B}}, M_{\text{up}}$ and $\Pi_{\mathcal{B}}, M_{\text{low}}$, and synthesize polynomial upper and lower bounds for the expected weights $[\Pi_{\mathcal{B}}, M_{\text{up}}]_{\mathbf{v}_i}(\mathbb{R}^{|\Pi_{\mathcal{B}}|})$ and $[\Pi_{\mathcal{B}}, M_{\text{low}}]_{\mathbf{v}_i}(\mathbb{R}^{|\Pi_{\mathcal{B}}|})$ for each initial program valuation $\mathbf{v}_i \in \mathcal{V}_i$ in $\mathcal{W}$ by solving the templates w.r.t the PUWF and PLWF constraints (i.e., (C1), (C2), (C1') from Definition 4.4), respectively. The correctness of this stage follows from Theorem 4.6 and Theorem 4.9.

Note that Theorem 4.6 and Theorem 4.9 have prerequisites and we check these prerequisites in a succinct fashion as follows. For Theorem 4.6, we manually verify the AST property by the approaches of [9, 12], and check the score-bounded property by a direct manual inspection. For

<footnote n="8">⁸We choose $d = 1, m = 30$ here to exemplify our algorithms, the values of $d, m$ for this example are larger in the experiments.</footnote>
**Theorem 4.9**, we check the condition (E1) by approaches such as [13, 50], the condition (E2) by a manual examination, and the bounded update property by simple check on whether the update value is bounded by a constant. The checking of the AST can be automated by template-based approaches [9, 12], and other simple prerequisites can be easily automated (via e.g. SMT solvers). We leave a fuller implementation that incorporate these features as a future work.

Below we present the detailed steps (**Step A1** – **Step A4)** of our template solving.

**Step A1.** Consider the truncated WPTS $\Pi_{\mathcal{B}, M_{\text{up}}}$ (to derive an upper bound) or $\Pi_{\mathcal{B}, M_{\text{low}}}$ (to derive a lower bound). We denote the bounded range from the truncation function $\mathcal{B}$ by $B := \{\mathbf{v} \mid \mathbf{v}(x) \in \mathcal{B}(x) \text{ for all } x \in V_{\text{p}}\}$. In this step, both of our algorithms compute an extended bounded range $B'$ such that any program valuation from a one-step execution of the WPTS under any current program valuation in $B$ falls in $B'$. Formally, the extended bounded range $B'$ satisfies that for every location $\ell$, program valuation $\mathbf{v} \in B$, transition $\tau = \langle \ell, \phi, F_1, F_2, \dots, F_k \rangle$, fork $F_j = \langle \ell'_j, p_j, \text{upd}_j, \text{wt}_j \rangle$ in $\tau$ and sampling valuation $\mathbf{r} \in \text{supp}(\mathcal{D})$, we have that $\text{upd}_j(\mathbf{v}, \mathbf{r}) \in B'$. Our algorithms determine the extended bounded range $B'$ by examining the assignment statements in the truncated WPTS $\Pi_{\mathcal{B}, M_{\text{up}}}$. The purpose to have a superset $B'$ of the original bounded range $B$ is to reduce the runtime in the solving of the template, see **Step A3** below.

*Example 5.2.* Recall Example 5.1, the bounded range from the truncation function $\mathcal{B}$ is denoted by $B = \{(\text{pos}, \text{dis}) \mid \text{pos} \in [0, 5], \text{dis} \in [0, 5]\}$. For $\Pi_{\mathcal{B}, M_{\text{up}}}$ or $\Pi_{\mathcal{B}, M_{\text{low}}}$, the extended bounded range is given by $B' = \{(\text{pos}, \text{dis}) \mid \text{pos} \in [-1, 6], \text{dis} \in [0, 6]\}$. $\square$

**Step A2.** In this step, at each location $\ell \notin \{\ell_{\text{out}}, \#\}$, our algorithms set up a $d$-degree polynomial template $h_\ell$ over the program variables $V_{\text{p}}$. Each template $h_\ell$ is a summation of all monomials of degree no more than $d$ and each monomial is multiplied with a fresh unknown coefficient as a parameter to be resolved. For $\ell \in \{\ell_{\text{out}}, \#\}$, our algorithms assume $h_\ell = 1$.

*Example 5.3.* Recall Example 5.1, the algorithm parameter $d=1$. For either $\Pi_{\mathcal{B}, M_{\text{up}}}$ or $\Pi_{\mathcal{B}, M_{\text{low}}}$, we construct a 1-degree polynomial template at the location $\ell_{\text{init}}$, i.e., $h_{\ell_{\text{init}}}(\text{pos}, \text{dis}) = a_1 \cdot \text{pos} + a_2 \cdot \text{dis} + a_3$ where $a_1, a_2, a_3 \in \mathbb{R}$ are unknown coefficients to be resolved, and let $h_{\ell_{\text{out}}} = h_{\ell_{\#}} = 1$. $\square$

**Step A3.** In this step, our algorithms establish constraints for the templates $h_\ell$'s from (C1), (C1') in Definition 4.4 (as (C2) is satisfied by the construction of $h_\ell$ in the previous step). For upper bounds on expected weights, our algorithms have the following relaxed constraints of (C1) to synthesize a PUWF over $\Pi_{\mathcal{B}, M_{\text{up}}}$:

(D1) For every location $\ell \in L \setminus \{\ell_{\text{out}}, \#\}$ and program valuation $\mathbf{v} \in I(\ell) \cap B$, we have that $\text{ewt}(h)(\ell, \mathbf{v}) \le h(\ell, \mathbf{v})$.

(D2) For every location $\ell \in L \setminus \{\ell_{\text{out}}, \#\}$ and program valuation $\mathbf{v} \in I(\ell) \cap (B' \setminus B)$, we have that $M_{\text{up}}(\mathbf{v}) \le h(\ell, \mathbf{v})$.

For lower bounds over $\Pi_{\mathcal{B}, M_{\text{low}}}$, our algorithms have the relaxed PLWF constraints (D1') and (D2') which are obtained from (D1) and resp. (D2) by replacing "$\text{ewt}(h)(\ell, \mathbf{v}) \le h(\ell, \mathbf{v})$" with "$\text{ewt}(h)(\ell, \mathbf{v}) \ge h(\ell, \mathbf{v})$" in (D1) and resp. "$M_{\text{up}}(\mathbf{v}) \le h(\ell, \mathbf{v})$" with "$M_{\text{low}}(\mathbf{v}) \ge h(\ell, \mathbf{v})$" in (D2), respectively. We have that (D1) and (D2) together ensure (C1) since $M_{\text{up}}(\mathbf{v}) \le h(\ell, \mathbf{v})$ implies that $\text{ewt}(h)(\ell, \mathbf{v}) \le h(\ell, \mathbf{v})$ for every location $\ell \in L \setminus \{\ell_{\text{out}}, \#\}$ and program valuation $\mathbf{v} \in I(\ell) \cap (B' \setminus B)$. The same holds for (D1') and (D2').

Note that in (D1), the calculation of the value $\text{ewt}(h)(\ell, \mathbf{v})$ has the piecewise nature that different sampling valuations $\mathbf{r}$ may cause the next program valuation to be within or outside the bounded range, and to satisfy or violate the guards of the transitions in the WPTS. In our algorithms, we have a fine-grained treatment for (D1) that enumerates all possible situations for a sampling valuation $\mathbf{r}$ that satisfy different guards of the WPTS in the calculation of $\text{ewt}(h)(\ell, \mathbf{v})$, for which we use an
SMT solver (e.g., Z3 [16]) to compute the situations. As for (D2), we use (D2) to avoid handling the piecewise feature in the computation of ewt(h)(ℓ, v) from within/outside the bounded range (i.e., the computation is a direct computation over a single-piece polynomial), so that the amount of computation of ewt(h)(ℓ, v) is reduced by ignoring the piecewise feature. The use of (D2) to reduce the computation is the aim of introducing an extended bounded range in **Step A1**. The same also holds for (D1') and (D2').

If $\Pi$ is non-score-recursive and the score function $g$ at termination is non-polynomial, our algorithms replace $g$ in the expression ewt(h)(ℓ, v) with its polynomial approximation $g'$.

**Example 5.4.** Recall Example 5.1–Example 5.3, the bounded range $D_1 = I(\ell_{init}) \cap B$ in (D1) is defined such that $pos \in [0, 5] \land dis \in [0, 5]$. Consider to derive upper bounds from $\Pi_{\mathcal{G}, M_{up}}$. We make a fine-grained treatment for (D1) by splitting the range $D_1$ and enumerating all possible situations for the sampling valuation **sample uniform**(0, 1) that the next program valuation satisfies or violates the loop guard “$pos \ge 0$”. In detail, $D_1$ is split into $D_{11} = \{(pos, dis) \mid pos \in [0, 1), dis \in [0, 5]\}$ and $D_{12} = \{(pos, dis) \mid pos \in (1, 5], dis \in [0, 5]\}$, where $D_{11}$ stands for the situation that with different sampling valuations the next program valuation may satisfy the loop guard (so that the next location is $\ell_{init}$) or violate the loop guard (so that the next location is directed to $\ell_{out}$), and $D_{12}$ stands for the situation that the next program valuation will definitely satisfy the loop guard and the next location is $\ell_{init}$. For the non-polynomial score function $g$ at the termination of the example, we replace it with a polynomial approximation $g'$ under the error bound $10^{-5}$.

We first show the PUWF constraints for $\Pi_{\mathcal{G}, M_{up}}$. Recall that the program has two probabilistic branches with probability 0.5. When the current program valuation is in $D_{11}$, we observe that

* if the loop takes the branch $pos := pos + step$, then the next value of $pos$ remains to be non-negative and the loop continues, and

* if the loop takes the branch $pos := pos - step$, then it is both possible that the next value of $pos$ satisfy or violate the loop guard, depending on the exact value of $step \in [0, 1]$.

Then we have the constraint (D1.1) over $D_{11}$ that has expectation over a piecewise function on step derived from whether the loop terminates in the next iteration or not, as follows:

$$ (D1.1) \quad \forall pos, dis. \ pos \in [0, 1] \land dis \in [0, 5] \Rightarrow \\ 0.5 \cdot \mathbb{E}_{step} [ [pos - step \ge 0] \cdot h_{\ell_{init}}(pos - step, dis + step) + [pos - step < 0] \cdot g'(dis)] \\ + 0.5 \cdot \mathbb{E}_{step} [h_{\ell_{init}}(pos + step, dis + step)] \le h_{\ell_{init}}(pos, dis). $$

When the current program valuation is in $D_{12}$, we observe that the next program valuation is guaranteed to satisfy the loop guard, and hence we have the constraint $D_{12}$ as follows:

$$ (D1.2) \quad \forall pos, dis. \ pos \in (1, 5] \land dis \in [0, 5] \Rightarrow \\ 0.5 \cdot \mathbb{E}_{step} [h_{\ell_{init}}(pos - step, dis + step)] + 0.5 \cdot \mathbb{E}_{step} [h_{\ell_{init}}(pos + step, dis + step)] \le h_{\ell_{init}}(pos, dis). $$

The range $D_2 = I(\ell_{init}) \cap (B' \setminus B)$ in (D2) is represented by the disjunctive formula $\Phi := (pos \in [0, 6] \land dis \in [5, 6]) \lor (pos \in [5, 6] \land dis \in [0, 6])$. For the program valuation in $D_2$, its next location is $\ell_{\#}$. Recall the constant approximation functions $M_{up} = 2.1 \times 10^{-330}$ and $M_{low} = 0$. From (D2), we have two constraints from the disjunctive clauses in $\Phi$ as follows:

$$ (D2.1) \quad \forall pos, dis. \ pos \in [0, 6] \land dis \in [5, 6] \Rightarrow 2.1 \times 10^{-330} \le h_{\ell_{init}}(pos, dis). $$

$$ (D2.2) \quad \forall pos, dis. \ pos \in [5, 6] \land dis \in [0, 6] \Rightarrow 2.1 \times 10^{-330} \le h_{\ell_{init}}(pos, dis). $$

For $\Pi_{\mathcal{G}, M_{low}}$, the PLWF constraints (D1'.1) to (D2'.2) are obtained from the PUWF constraints above by replacing “$\le$” with “$\ge$” and $2.1 \times 10^{-330}$ with 0. $\square$

**Step A4.** In this step, for every initial program valuation $v_i \in V_i$ in $W$ where $V_i$'s and $W$ are obtained from **Stage 2**, our algorithms solve the unknown coefficients in the templates $h_\ell$
($\ell \in L \setminus \{\ell_{out}, \#\}$) via the well-established methods of Putinar's Positivstellensatz [36] or Handelman's Theorem [24]. To be more detailed, our algorithms minimize (resp. maximize) the objective function $h_{\ell_{init}}(\mathbf{v}_i)$ for each $\mathbf{v}_i$ in $\mathcal{W}$ which subjects to the PUWF (resp. PLWF) constraints from the previous step to derive polynomial upper (resp. lower) bounds over $\Pi_{\mathcal{B}}, M_{up}$ (resp. $\Pi_{\mathcal{B}}, M_{low}$), respectively.

Note that the PUWF or PLWF constraints from **Step A4** can be represented as a conjunction of formulas in the form $\forall \mathbf{v} \in \mathcal{P}. (g(\mathbf{v}) \ge 0)$ where the set $\mathcal{P}$ is defined by a conjunction of polynomial inequalities in the program variables and $g$ is a polynomial over $V_p$ whose coefficients are affine expressions in the unknown coefficients from the templates, and such formulas can be guaranteed by the sound forms of Putinar's Positivstellensatz and Handelman's Theorem. The application of Putinar's Positivstellensatz results in semidefinite constraints and can be solved by semidefinite programming (SDP), while the application of Handelman's Theorem is restricted to the affine case (i.e., every condition and assignment in the WPTS is affine) and leads to linear constraints and can be solved by linear programming (LP). We refer to Appendix C for the details on the application of Putinar's Positivstellensatz and Handelman's Theorem.

If the polynomial approximation is applied in the calculation of $\mathrm{ewt}(h)(\ell, \mathbf{v})$ in **Step A3**, we show that the error of the final results will be bounded by the polynomial approximation error bound $\epsilon$. See Theorem C.1 in the Appendix C.

**Example 5.5.** Recall Example 5.1–Example 5.4, we have that $\mathcal{V}_1 = \{(\mathrm{pos}, \mathrm{dis}) \mid \mathrm{pos} \in [0, 0.1], \mathrm{dis} = 0\}, \dots, \mathcal{V}_{30} = \{(\mathrm{pos}, \mathrm{dis}) \mid \mathrm{pos} \in [2.9, 3], \mathrm{dis} = 0\}$ and the set $\mathcal{W} = \{(0.05, 0), (0.15, 0), \dots, (2.95, 0)\}$. Pick a point $\mathbf{v}_1 = (0.05, 0) \in \mathcal{V}_1$ from $\mathcal{W}$, to synthesize the upper bound for $[\Pi_{\mathcal{B}}, M_{up}]_{\mathbf{v}_1}(\mathbb{R}^{|V_p|})$ as tight as possible, we solve the following optimization problem whose objective function is $h_{\ell_{init}}(0.05, 0)$, i.e.,

Minimize $0.05 \cdot a_1 + a_3$

Subject to constraints (D1.1)–(D2.2)

The lower bound for $[\Pi_{\mathcal{B}}, M_{low}]_{\mathbf{v}_1}(\mathbb{R}^{|V_p|})$ is solved by

Maximize $0.05 \cdot a_1 + a_3$

Subject to constraints (D1'.1)–(D2'.2)

Although the constraints are universally quantified, the universal quantifiers can be soundly (but not completely) removed and relaxed into semidefinite constraints over the unknown coefficients $a_i$'s by applying Putinar's Positivstellensatz, where we over-approximate all strict inequalities (e.g., "<") by non-strict ones (e.g., "$\le$"). Then we call a SDP solver to solve the two optimization problems and find the solutions of $a_i$'s, which will generate two polynomial bound functions $Up_1, Lw_1$ such that $[\Pi_{\mathcal{B}}, M_{up}]_{\mathbf{v}}(\mathbb{R}^{|V_p|}) \le Up_1(\mathbf{v}) + \epsilon$ and $Lw_1(\mathbf{v}) - \epsilon \le [\Pi_{\mathcal{B}}, M_{low}]_{\mathbf{v}}(\mathbb{R}^{|V_p|})$ for all program valuation $\mathbf{v} \in \mathcal{V}_1$ where the error bound $\epsilon = 10^{-5}$. $\square$

**Stage 4: Integration.** As a consequence of **Stage 3**, our algorithms obtain a group of concrete polynomial upper bound functions (resp. lower bound functions) $Upper = \{Up_1, \dots, Up_m\}$ (resp. $Lower = \{Lw_1, \dots, Lw_m\}$) for the expected weights $\{[\Pi_{\mathcal{B}}, M_{up}]_{\mathbf{v}_1}(\mathbb{R}^{|V_p|}), \dots, [\Pi_{\mathcal{B}}, M_{up}]_{\mathbf{v}_m}(\mathbb{R}^{|V_p|})\}$ (resp. $\{[\Pi_{\mathcal{B}}, M_{low}]_{\mathbf{v}_1}(\mathbb{R}^{|V_p|}), \dots, [\Pi_{\mathcal{B}}, M_{low}]_{\mathbf{v}_m}(\mathbb{R}^{|V_p|})\}$), respectively. Then our algorithms integrate these polynomial upper and lower bound functions to derive the upper and lower bounds for $[\Pi_{\mathcal{B}}, M_{up}](\mathbb{R}^{|V_p|})$ and $[\Pi_{\mathcal{B}}, M_{low}](\mathbb{R}^{|V_p|})$, respectively. For $[\Pi_{\mathcal{B}}, M_{up}](\mathbb{R}^{|V_p|})$, we have that

$$[\Pi_{\mathcal{B}}, M_{up}](\mathbb{R}^{|V_p|}) \le \sum_{i=1}^{m} \int_{V_i} Up_i(\mathbf{v})dv =: u_Z \quad (2)$$
and the lower bound for $\llbracket \Pi_B, M_{\text{low}} \rrbracket (\mathbb{R}^{|V_P|})$ is given by

$$l_Z := \sum_{i=1}^{m} \int_{V_i} L w_i(v) dv \le \llbracket \Pi_B, M_{\text{low}} \rrbracket (\mathbb{R}^{|V_P|}) \quad (3)$$

If $\Pi$ is non-score-recursive and the score function $g$ at termination is non-polynomial, our algorithms integrate the approximation error $\zeta = \text{volume}(\mathcal{V}) \cdot \epsilon$ caused by polynomial approximation of $g$ to the two bounds, i.e, $l'_Z = l_Z - \zeta$ and $u'_Z = u_Z + \zeta$, where volume($\mathcal{V}$) is the volume of $\mathcal{V}$. In practice, to ensure the tightness of the interval bounds, we can control the amount of $\epsilon$ so that the approximation error $\zeta$ is at least one magnitude smaller than the values of $l_Z, u_Z$.

**THEOREM 5.6 (SOUNDNESS).** If our algorithms find valid solutions for the unknown coefficients of the templates, they return correct interval bounds for the normalising constant $Z_{\Pi}$.

**PROOF SKETCH.** By Theorem 4.6 (resp. Theorem 4.9), if the algorithms successfully find valid solutions for the unknown coefficients of the templates, we can obtain the polynomial upper bound $u_Z$ (resp. lower bound $l_Z$) for $\llbracket \Pi_B, M_{\text{up}} \rrbracket (\mathbb{R}^{|V_P|})$ (resp. $\llbracket \Pi_B, M_{\text{low}} \rrbracket (\mathbb{R}^{|V_P|})$), respectively. Then by Theorem 4.11, the polynomial upper bound $u_Z$ is also the upper bound for $\llbracket \Pi \rrbracket (\mathbb{R}^{|V_P|})$ and $l_Z$ is the lower bound for $\llbracket \Pi \rrbracket (\mathbb{R}^{|V_P|})$. The same holds for the bounds $l'_Z, u'_Z$ if polynomial approximation of non-polynomial score functions happens. $\square$

# 6 EXPERIMENTAL RESULTS

In this section, we present the experimental valuation of our approach over a variety of benchmarks. First, we show that our approach can handle novel examples that cannot be addressed by other existing tools such as [19, 20, 25, 33]. Then we compare our approach with the state-of-the-art tool GuBPI [4]. Finally, even though the problem of path probability estimations is not the focus of our work, we demonstrate that our approach can work well for this problem, and we also compare the performance of our approach with GuBPI.

We have implemented a parser from probabilistic programs to WPTS's in F#, our algorithms in Matlab, and used the LP solver in Matlab (resp. Mosek [1]) for solving linear (resp. semidefinite) programming, respectively. All results were obtained on an Intel Core i7 (2.3 GHz) machine with 16 GB of memory, running MS Windows 10.

## 6.1 Experimental Setup

Below we clarify our experimental setting for the input programs with invariant, the input truncation and approximation functions and the polynomial approximation for the score statement at program termination in the non-score-recursive case.

*Program Input.* All the benchmarks in our experiments are of the form “while $\phi$ do $S_1$ od $S_2$” where $S_1, S_2$ are program statements without loops and $\phi$ is the loop guard which is a predicate over $V_P$. We set two locations for all benchmarks, i.e., $\ell_{\text{init}}$ before **while** and $\ell_{\text{out}}$ after $S_2$. Our algorithms also work for general programs with multiple locations. To reduce the influence from the choice of the invariant, we derive an invariant $I$ simply from the loop guard $\phi$ so that $I(\ell_{\text{init}}) = \phi$ and $I(\ell_{\text{out}}) = \neg \phi$.

*Partition of Initial Program Valuations.* We partition the set $\mathcal{V}$ of initial program valuations into $m$ disjoint subsets by splitting the quantity of one dimension uniformly.

*Truncation Function.* For the truncation function $\mathcal{B}$, we first restrict the range of all program variables to the domain satisfying the loop guard $\phi$. That is, $\Phi_1 = \{v | v \models \phi\}$. Then we empirically specify a large enough bounded range $\Phi_2$ over program variables that captures the major behaviour
of the program. That is, $\Phi_2 = \{\mathbf{v} \mid \mathbf{v} \in [\zeta_1 - \delta, \zeta_2 + \delta], [\zeta_1, \zeta_2] = \text{supp}(\mu_{init})\}$ for some constant $\delta > 0$. Finally, $\mathcal{B}$ is defined by $\mathcal{B}(x) := \Phi_1(x) \wedge \Phi_2(x)$ for all $x \in V_p$.

*Approximation Functions.* To derive the approximation functions $M_{up}$, $M_{low}$, we consider two cases:

* in the case that $\Pi$ is non-score-recursive, $M_{up}$, $M_{low}$ are obtained as follows.

- If the score function at termination is non-polynomial, we determine $M_{up}$, $M_{low}$ by bounds and also the monotonicity of the score function over the set $B'' = B' \setminus B$ (see definitions of $B$, $B'$ in **Step A1** in Section 5).

- Otherwise, we directly use the polynomial score function as $M_{up}$, $M_{low}$.

* in the case that $\Pi$ is score-recursive, $M_{up}$, $M_{low}$ are derived by using our OST approach in Section 4.2 but without truncation operation, i.e., by directly applying Theorem 4.9.

*Polynomial Approximation.* In the case that $\Pi$ is non-score-recursive and the score function $g$ at termination is non-polynomial, there are two situations that polynomial approximation should be applied. The first is that $g$ is handled by polynomial approximation of it over the set $B'$, which will be later used in (D1) in **Step A3** in Section 5. The second is that $M_{up}$, $M_{low}$ are determined by polynomial approximation of $g$ over the set $B''$, which will be used in (D2) in **Step A3** in Section 5. Concretely, we apply polynomial interpolation to approximate the non-polynomial score function $g$. That is, given an error bound $\epsilon > 0$ we aim to find a polynomial approximation $g'$ for $g$ such that $|g(\mathbf{v}) - g'(\mathbf{v})| \le \epsilon'$ for all program valuations $\mathbf{v} \in B'$ (or $\mathbf{v} \in B''$). The above operation can be easily extended to produce piecewise polynomial approximations if $B$, $B''$ are split into multiple partitions. For the first situation, the approximation error caused by polynomial interpolation is taken into account when calculating interval bounds for NPD (see Theorem C.1 in Appendix C).

## 6.2 Results

*NPD - Novel Examples.* We consider 10 novel examples adapted from the literature, where all 7 examples with prefix “Pd” or “RDWALK” are from [4], the two “RACE” examples are from [50], and the last example is from statistical phylogenetics [38] (see also Section 3). Concretely, the “RACE(V2)” and “BIRTH” examples are both score-recursive probabilistic programs with weights greater than 1, and thus their integrability condition should be verified by the existence of suitable concentration bounds (see Theorem 4.8); other examples are non-score-recursive probabilistic while loops with unsupported types of scoring by previous tools (e.g., polynomial scoring). Therefore, no existing tools w.r.t. NPD can tackle these novel examples. The results are reported in Table 1, where the first column is the name of each example, the second column contains the parameter of each example used in our approach (i.e., the degree of the polynomial template and the bounded range of program variables), the third column is the used solver, and the fourth and fifth columns correspond to the runtime of upper and lower bounds computed by our approach, respectively. We set the partition number $m = 60$. Our runtime is reasonable, that is, most examples can obtain tight bounds within 100 seconds, and the simulation results by Pyro [5] ($10^6$ samples per case) match our derived bounds. Due to space limitation, we display part of the comparison in Fig. 8, see Appendix D for other figures.

*NPD - Comparison with GuBPI [4].* Since the parameters used in GuBPI and our approach are completely different, it is infeasible to compare the two approaches directly. Instead, we choose the parameters to our algorithms that can achieve at least comparable results with GuBPI. The main parameters are shown in Table 2 and we set the partition number $m = 60$. We consider the Pedestrian example “Pd” from [4] (see also Section 3.1), and its variants. More concretely, we enlarged the standard deviation of the observed normal distribution to be 5 for all other 6 examples
<table><caption>Table 1. Results for Novel Examples</caption><thead><tr><th rowspan="2">Benchmark</th><th rowspan="2">Parameters</th><th rowspan="2">Solver</th><th>Upper</th><th>Lower</th></tr><tr><th>Time (s)</th><th>Time (s)</th></tr></thead><tbody><tr><td>P<sub>D</sub>(v1)</td><td>d = 6, pos, dis ∈ [0, 5]</td><td>SDP</td><td>54.65</td><td>52.39</td></tr><tr><td>R<sub>ACE</sub>(v1)</td><td>d = 6, h, t ∈ [0, 5]</td><td>LP</td><td>87.43</td><td>86.27</td></tr><tr><td>R<sub>ACE</sub>(v2)*</td><td>d = 6, h, t ∈ [0, 5]</td><td>LP</td><td>81.19</td><td>81.18</td></tr><tr><td>R<sub>DWALK</sub>(v1)</td><td>d = 6, x, y ∈ [0, 5]</td><td>LP</td><td>46.65</td><td>47.72</td></tr><tr><td>R<sub>DWALK</sub>(v2)</td><td>d = 6, x, y ∈ [0, 5]</td><td>LP</td><td>97.45</td><td>103.65</td></tr><tr><td>R<sub>DWALK</sub>(v3)</td><td>d = 6, x, y ∈ [0, 5]</td><td>LP</td><td>250.70</td><td>252.49</td></tr><tr><td>R<sub>DWALK</sub>(v4)</td><td>d = 6, x, y ∈ [0, 5]</td><td>LP</td><td>98.75</td><td>98.21</td></tr><tr><td>P<sub>DMB</sub>(v3)</td><td>d = 4, pos, dis ∈ [0, 5]</td><td>LP</td><td>15.26</td><td>14.57</td></tr><tr><td>P<sub>DMB</sub>(v4)</td><td>d = 4, pos, dis ∈ [0, 5]</td><td>LP</td><td>16.07</td><td>16.12</td></tr><tr><td>B<sub>IRTH</sub>*</td><td>d = 6, lambda ∈ [0, 3], time ∈ [0, 10]</td><td>LP</td><td>14.72</td><td>16.66</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>(a) P<sub>D</sub>(v1)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) R<sub>ACE</sub>(v2)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(a) P<sub>D</sub></figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) P<sub>DLD</sub></figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) R<sub>DWALK</sub>(v1)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) P<sub>DMB</sub>(v3)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) P<sub>DBETA</sub>(v1)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(d) P<sub>DBETA</sub>(v2)</figcaption></figure>

<caption>Fig. 8. NPD Bounds of Novel Examples</caption>

<caption>Fig. 9. NPD Bounds of Comparison</caption>

The red and the blue lines mark the upper and lower bounds of our results; the black bold stars mark the simulation results; the brown and green dotted lines mark the upper and lower bounds generated by GuBPI (we denote by -0.1 the infinity bounds).

whose prefix name are “P<sub>D</sub>”; for the four “P<sub>DBETA</sub>” examples, we also add different beta distributions in the loop bodies. The last example is from [19]. We report the results in Table 2 whose layout is similar to Table 1 except that the column “#” displays whether or not the bounds are trivial, i.e., [0, ∞]. We also compare our results with GuBPI’s and simulation results (10⁶ samples per case), and show part of the comparison in Fig. 9, see Appendix D for other figures. Our runtime is up to 6 times faster than GuBPI while we can still obtain tighter or comparable bounds for all examples. Specifically, for the first example “P<sub>D</sub>”, our upper bounds are a bit higher than GuBPI’s when the value of start falls into [0, 0.7] (which is not suprising as the deviation of the normal distribution

<footnote>* It is a score-recursive probabilistic program with weights greater than 1.</footnote>
<table><caption>Table 2. Comparison with GuBPI</caption><thead><tr><th rowspan="2">Benchmark</th><th colspan="4">Our Tool</th><th colspan="2">GuBPI</th></tr><tr><th>Parameters</th><th>Solver</th><th>Time (s)</th><th>#</th><th>Time (s)</th><th>#</th></tr></thead><tbody><tr><td>P<sub>D</sub></td><td>d = 10, pos, dis ∈ [0, 5]</td><td>SDP</td><td>3176.685</td><td>•</td><td>5266.063</td><td>•</td></tr><tr><td>P<sub>DLD</sub></td><td>d = 6, pos, dis ∈ [0, 5]</td><td>LP</td><td>41.99</td><td>•</td><td>648.151</td><td>•</td></tr><tr><td>P<sub>D</sub>B<sub>E</sub>T<sub>A</sub>(v1)</td><td>d = 6, pos, dis ∈ [0, 5]</td><td>LP</td><td>99.86</td><td>•</td><td>645.055</td><td>○</td></tr><tr><td>P<sub>D</sub>B<sub>E</sub>T<sub>A</sub>(v2)</td><td>d = 6, pos, dis ∈ [0, 5]</td><td>LP</td><td>228.43</td><td>•</td><td>653.237</td><td>○</td></tr><tr><td>P<sub>D</sub>B<sub>E</sub>T<sub>A</sub>(v3)</td><td>d = 6, pos, dis ∈ [0, 5]</td><td>LP</td><td>101.36</td><td>•</td><td>657.645</td><td>○</td></tr><tr><td>P<sub>D</sub>B<sub>E</sub>T<sub>A</sub>(v4)</td><td>d = 6, pos, dis ∈ [0, 5]</td><td>LP</td><td>208.86</td><td>•</td><td>686.207</td><td>○</td></tr><tr><td>P<sub>D</sub>M<sub>B</sub>(v5)</td><td>d = 6, pos, dis ∈ [0, 5]</td><td>LP</td><td>88.41</td><td>•</td><td>391.772</td><td>•</td></tr><tr><td>P<sub>A</sub>R<sub>A</sub>-R<sub>E</sub>C<sub>R</sub>U<sub>R</sub></td><td>d = 8, p ∈ [0, 1]</td><td>LP</td><td>36.61</td><td>•</td><td>253.728</td><td>•</td></tr></tbody></table>

* ○ marks the trivial bound [0, ∞], while ● marks the non-trivial ones.

in this example is quite small, i.e., 0.1, and our approach constructs over-approximation constraints while GuPBI uses recursion unrolling to search for the feasible space exhaustively), but our lower bounds are greater than GuBPI's, and our NPD bounds are tighter in the following.<fnref n="9" /> For all 6 variants of "P<sub>D</sub>" where the deviation of the normal distribution is enlarged, our NPD bounds are tighter than GuBPI's, in particular, our upper bounds are much lower than GuBPI's. For the four "P<sub>D</sub>B<sub>E</sub>T<sub>A</sub>" examples, we also found that GuBPI produced zero-valued unnormalized lower bounds, thus its results w.r.t. NPD are trivial, i.e., [0, ∞]. However, we can still produce non-trivial results and our runtime is at least 2 times faster than GuBPI.

*Path Probability Estimation.* We consider five recursive examples in [4], which were also cited from the PSI repository [19]. Since all five examples are non-parametric and with unbounded numbers of loop iterations, PSI cannot handle them correctly as mentioned in [4]. We estimated the path probability of certain events, i.e., queries over program variables, and report the results in Table 3. For the first three examples, we obtained tighter lower bounds than GuBPI and same upper bounds, while our runtime is at least 2 times faster than GuBPI. Moreover, we found a potential error of GuBPI. That is, the fourth example "CAV-EX-5" in Table 3 is an AST program with no scores, which means its normalizing constant should be exactly one. However, the upper bound of the normalising constant obtained by GuBPI is smaller than 1 (i.e., 0.6981). A stochastic simulation using 10⁶ samples yielded the results that fall within our bounds but violate those computed by GuBPI. Thus, GuBPI possibly omitted some valid program runs of this example and produces wrong results. All our results match the simulation results (10⁶ samples per case).

# 7 RELATED WORKS

Below we compare our results with the most related work in the literature.

*Static analysis in Bayesian probabilistic programming.* There are a lot of works on NPD inference for probabilistic programs, such as (λ)PSI [19, 20], AQUA [25], Hakaru [33] and SPPL [41]. However, these methods are restricted to specific kinds of programs, e.g., programs with closed-form solutions

<footnote n="9">⁹When the value of *start* approaches 3, our NPD bounds are close to zero, but the upper bounds may be lower than zero, which is caused by numerical issues of semi-definite programming. The problem of numerical issues is orthogonal to our work and remains to be addressed in both academic and industrial fields.</footnote>
<table><caption>Table 3. Results for Path Probability Estimation</caption><thead><tr><th rowspan="2">Benchmark</th><th rowspan="2">Query</th><th colspan="3">Our Tool</th><th colspan="2">GuBPI</th><th rowspan="2">Simul</th></tr><tr><th>Parameters</th><th>Time (s)</th><th>Bounds</th><th>Time (s)</th><th>Bounds</th></tr></thead><tbody><tr><td rowspan="2">CAV-EX-7</td><td>Q1</td><td>6, [0, 30], [0, 4]</td><td>15.062</td><td>[0.9698, 1.0000]</td><td>38.834</td><td>[0.7381, 1.0000]</td><td>0.9938</td></tr><tr><td>Q2</td><td>6, [0, 40], [0, 4]</td><td>16.321</td><td>[0.9985, 1.0000]</td><td>37.651</td><td>[0.7381, 1.0000]</td><td>0.9993</td></tr><tr><td rowspan="2">ADDUNI(L)</td><td>Q1</td><td>6, [0, 10], [0, 1]</td><td>8.85</td><td>[0.9940, 1.0000]</td><td>21.064</td><td>[0.9375, 1.0000]</td><td>0.9991</td></tr><tr><td>Q2</td><td>6, [0, 15], [0, 1]</td><td>8.80</td><td>[0.9995, 1.0000]</td><td>14.941</td><td>[0.9375, 1.0000]</td><td>0.9999</td></tr><tr><td rowspan="2">RDBOX</td><td>Q1</td><td>4, [-0.8, 0.8], [0, 10]</td><td>25.87</td><td>[0.9801, 1.0000]</td><td>173.535</td><td>[0.9462, 1.0000]</td><td>0.9999</td></tr><tr><td>Q1</td><td>6, [20, ∞], [0, 10]</td><td>33.17</td><td>[0.8123, 0.9707]</td><td>229.623</td><td>[0.5768, 0.6374]</td><td>0.9098</td></tr><tr><td rowspan="2">CAV-EX-5*</td><td>Q2</td><td>6, [20, ∞], [0, 20]</td><td>54.373</td><td>[0.8970, 1.0000]</td><td>224.504</td><td>[0.5768, 0.6375]</td><td>0.9645</td></tr><tr><td>Q1</td><td>8, [1, ∞], [0, 0.1]</td><td>7.255</td><td>[0.0023, 0.0023]</td><td>33.246</td><td>[0.0023, 0.0024]</td><td>0.0023</td></tr><tr><td rowspan="2">GWALK**</td><td>Q2</td><td>8, [1, ∞], [0, 0.2]</td><td>8.197</td><td>[0.0025, 0.0025]</td><td>31.728</td><td>[0.0025, 0.0025]</td><td>0.0025</td></tr></tbody></table>

to NPD or without continuous distributions, and none of them can handle probabilistic programs with unbounded while-loops/recursion. As far as we know, the most revelant work on static analysis of posterior distribution over unbounded loops/recursion is the approach [4] that infers the bounds for posterior distributions by recursion unrolling and bounding the non-termination case via the widening operator of abstract interpretation. By unrolling recursion to arbitrary depth, this approach can achieve high precision on the derived bounds. However, a major drawback of this approach is that the recursion unrolling may cause path explosion. Our approach circumvents the path explosion problem by constraint solving. Another major drawback is that this approach cannot handle score-recursive programs as simply applying the approach to score-recursive programs leads to the trivial bound [0, ∞], and we address this issue by a novel OST variant.

**MCMC and variational inference.** As mentioned previously, statistical approaches such as MCMC [18, 39] and variational inference [6] cannot provide formal guarantee on the bounds for posterior distributions in a finite time limit. In contrast, our approach has formal guarantee on the derived bounds.

**Static analysis of probabilistic programs.** In recent years, there have been an abundance of works on static analysis of probabilistic programs. Most of them address fundamental aspects such as termination [9, 12, 17], sensitivity [2, 52], expectation [34, 53], tail bounds [27, 49, 51], assertion probability [43, 50], etc. Compared with these results, we have:

* Our work focuses on normalized posterior distribution in Bayesian probabilistic programming, and hence is an orthogonal objective.

* Our algorithm follows the previous works on the synthesis of polynomial templates [9, 10, 13, 53], but we have a truncation operation to increase the accuracy which to our best knowledge is novel.

* Our approach extends the classical OST as the previous works [49, 53] do, but we consider a multiplicative variant, while the work [53] considers only an additive variant, and the work [49] considers a general extension through the uniform integrability condition and an implementation via polynomial functions, but does not have a detailed treatment for a multiplicative variant.

A very recent work [3] considers the synthesis of piecewise bounds for probabilistic programs. We focus on non-piecewise polynomial bounds and hence is orthogonal. A promising future direction would be also to consider piecewise bounds in the NPD problem.

<footnote>* GuBPI's result contradicts ours, and we found GuBPI produces wrong results for this example.</footnote>

<footnote>** As we care about path probabilities, we compared bounds of unnormalized distributions for this example (the NPD can be derived in the same manner above).</footnote>
# 8 CONCLUSION

In this work, we considered the formal analysis of normalized posterior distribution in Bayesian probabilistic programming. Our contribution is a novel template-based approach that circumvents the path explosion problem from loop/recursion unrolling and addresses score-recursive programs via a novel variant of Optional Stopping Theorem. A future direction would be to investigate whether our approach can be further improved by piecewise polynomials.
REFERENCES

[1] MOSEK ApS. 2022. *The MOSEK optimization toolbox for MATLAB manual. Version 10.0.* http://docs.mosek.com/10.0/toolbox/index.html

[2] Gilles Barthe, Thomas Espitau, Benjamin Grégoire, Justin Hsu, and Pierre-Yves Strub. 2018. Proving expected sensitivity of probabilistic programs. *Proc. ACM Program. Lang.* 2, POPL (2018), 57:1-57:29. https://doi.org/10.1145/3158145

[3] Kevin Batz, Mingshuai Chen, Sebastian Junges, Benjamin Lucien Kaminski, Joost-Pieter Katoen, and Christoph Matheja. 2023. Probabilistic Program Verification via Inductive Synthesis of Inductive Invariants. In *Tools and Algorithms for the Construction and Analysis of Systems - 29th International Conference, TACAS 2023, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2022, Paris, France, April 22-27, 2023, Proceedings, Part II (Lecture Notes in Computer Science, Vol. 13994)*, Sriram Sankaranarayanan and Natasha Sharygina (Eds.). Springer, 410-429. https://doi.org/10.1007/978-3-031-30820-8_25

[4] Raven Beutner, C.-H. Luke Ong, and Fabian Zaiser. 2022. Guaranteed bounds for posterior inference in universal probabilistic programming. In *PLDI '22: 43rd ACM SIGPLAN International Conference on Programming Language Design and Implementation, San Diego, CA, USA, June 13-17, 2022*, Ranjit Jhala and Isil Dillig (Eds.). ACM, 536-551. https://doi.org/10.1145/3519939.3523721

[5] Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul A. Szerlip, Paul Horsfall, and Noah D. Goodman. 2019. Pyro: Deep Universal Probabilistic Programming. *J. Mach. Learn. Res.* 20 (2019), 28:1-28:6. http://jmlr.org/papers/v20/18-403.html

[6] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. 2017. Variational inference: A review for statisticians. *Journal of the American statistical Association* 112, 518 (2017), 859-877.

[7] Johannes Borgström, Ugo Dal Lago, Andrew D. Gordon, and Marcin Szymczak. 2016. A lambda-calculus foundation for universal probabilistic programming. In *Proceedings of the 21st ACM SIGPLAN International Conference on Functional Programming, ICFP 2016, Nara, Japan, September 18-22, 2016*, Jacques Garrigue, Gabriele Keller, and Eijiro Sumii (Eds.). ACM, 33-46. https://doi.org/10.1145/2951913.2951942

[8] Aleksandar Chakarov and Sriram Sankaranarayanan. 2013. Probabilistic Program Analysis with Martingales. In *CAV 2013*, 511-526.

[9] Aleksandar Chakarov and Sriram Sankaranarayanan. 2013. Probabilistic Program Analysis with Martingales. In *Computer Aided Verification - 25th International Conference, CAV 2013, Saint Petersburg, Russia, July 13-19, 2013. Proceedings (Lecture Notes in Computer Science, Vol. 8044)*, Natasha Sharygina and Helmut Veith (Eds.). Springer, 511-526. https://doi.org/10.1007/978-3-642-39799-8_34

[10] Krishnendu Chatterjee, Hongfei Fu, and Amir Kafshdar Goharshady. 2016. Termination Analysis of Probabilistic Programs Through Positivstellensatz's. In *CAV 2016*, 3-22.

[11] Krishnendu Chatterjee, Hongfei Fu, Amir Kafshdar Goharshady, and Ehsan Kafshdar Goharshady. 2020. Polynomial invariant generation for non-deterministic recursive programs. In *Proceedings of the 41st ACM SIGPLAN International Conference on Programming Language Design and Implementation, PLDI 2020*, London, UK, June 15-20, 2020, Alastair F. Donaldson and Emina Torlak (Eds.). ACM, 672-687. https://doi.org/10.1145/3385412.3385969

[12] Krishnendu Chatterjee, Hongfei Fu, Petr Novotný, and Rouzbeh Hasheminezhad. 2016. Algorithmic analysis of qualitative and quantitative termination problems for affine probabilistic programs. In *Proceedings of the 43rd Annual ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages, POPL 2016*, St. Petersburg, FL, USA, January 20-22, 2016, Rastislav Bodík and Rupak Majumdar (Eds.). ACM, 327-342. https://doi.org/10.1145/2837614.2837639

[13] Krishnendu Chatterjee, Hongfei Fu, Petr Novotný, and Rouzbeh Hasheminezhad. 2018. Algorithmic Analysis of Qualitative and Quantitative Termination Problems for Affine Probabilistic Programs. *ACM Trans. Program. Lang. Syst.* 40, 2 (2018), 7:1-7:45. https://doi.org/10.1145/3174800

[14] Krishnendu Chatterjee, Petr Novotný, and Dörde Žikelić. 2017. Stochastic invariants for probabilistic termination. In *POPL 2017*, 145-160.

[15] Michael Colón, Sriram Sankaranarayanan, and Henny Sipma. 2003. Linear Invariant Generation Using Non-linear Constraint Solving. In *Computer Aided Verification, 15th International Conference, CAV 2003, Boulder, CO, USA, July 8-12, 2003, Proceedings (Lecture Notes in Computer Science, Vol. 2725)*, Warren A. Hunt Jr. and Fabio Somenzi (Eds.). Springer, 420-432. https://doi.org/10.1007/978-3-540-45069-6_39

[16] Leonardo De Moura and Nikolaj Bjørner. 2008. Z3: An Efficient SMT Solver. In *Proceedings of the Theory and Practice of Software, 14th International Conference on Tools and Algorithms for the Construction and Analysis of Systems* (Budapest, Hungary) (TACAS'08/ETAPS'08). Springer-Verlag, Berlin, Heidelberg, 337-340.

[17] Hongfei Fu and Krishnendu Chatterjee. 2019. Termination of Nondeterministic Probabilistic Programs. In *Verification, Model Checking, and Abstract Interpretation - 20th International Conference, VMCAI 2019, Cascais, Portugal, January 13-15, 2019, Proceedings (Lecture Notes in Computer Science, Vol. 11388)*, Constantin Enea and Ruzica Piskac (Eds.). Springer, 468-490. https://doi.org/10.1007/978-3-030-11245-5_22
[18] Dani Gamerman and Hedibert F Lopes. 2006. *Markov chain Monte Carlo: stochastic simulation for Bayesian inference*. CRC press.

[19] Timon Gehr, Sasa Misailovic, and Martin T. Vechev. 2016. PSI: Exact Symbolic Inference for Probabilistic Programs. In *Computer Aided Verification - 28th International Conference, CAV 2016, Toronto, ON, Canada, July 17-23, 2016, Proceedings, Part I (Lecture Notes in Computer Science, Vol. 9779)*, Swarat Chaudhuri and Azadeh Farzan (Eds.). Springer, 62–83.
https://doi.org/10.1007/978-3-319-41528-4_4

[20] Timon Gehr, Samuel Steffen, and Martin T. Vechev. 2020. λPSI: exact inference for higher-order probabilistic programs. In *Proceedings of the 41st ACM SIGPLAN International Conference on Programming Language Design and Implementation, PLDI 2020, London, UK, June 15-20, 2020*, Alastair F. Donaldson and Emina Torlak (Eds.). ACM, 883–897. https://doi.org/10.1145/3385412.3386006

[21] Noah D. Goodman, Vikash K. Mansinghka, Daniel M. Roy, Kallista A. Bonawitz, and Joshua B. Tenenbaum. 2008. Church: a language for generative models. In *UIAI 2008, Proceedings of the 24th Conference in Uncertainty in Artificial Intelligence, Helsinki, Finland, July 9-12, 2008*, David A. McAllester and Petri Myllymäki (Eds.). AUAI Press, 220–229.

[22] Noah D Goodman and Andreas Stuhlmüller. 2014. The Design and Implementation of Probabilistic Programming Languages. http://dippl.org.

[23] Andrew D Gordon, Thomas A Henzinger, Aditya V Nori, and Sriram K Rajamani. 2014. Probabilistic programming. In *Future of Software Engineering Proceedings*. 167–181.

[24] David Handelman. 1988. Representing polynomials by positive linear functions on compact convex polyhedra. *Pacific J. Math.* 132, 1 (1988), 35–62.

[25] Zixin Huang, Saikat Dutta, and Sasa Misailovic. 2021. AQUA: Automated Quantized Inference for Probabilistic Programs. In *Automated Technology for Verification and Analysis - 19th International Symposium, ATVA 2021, Gold Coast, QLD, Australia, October 18-22, 2021, Proceedings (Lecture Notes in Computer Science, Vol. 12971)*, Zhe Hou and Vijay Ganesh (Eds.). Springer, 229–246. https://doi.org/10.1007/978-3-030-88885-5_16

[26] H Jeffreys. 1988. "Weierstrass’s theorem on approximation by polynomials" and "Extension of Weierstrass’s approximation theory". *Methods of Mathematical Physics* (1988), 446–448.

[27] Satoshi Kura, Natsuki Urabe, and Ichiro Hasuo. 2019. Tail probabilities for randomized program runtimes via martingales for higher moments. In *International Conference on Tools and Algorithms for the Construction and Analysis of Systems*. Springer, 135–153.

[28] Wonyeol Lee, Hangyeol Yu, Xavier Rival, and Hongseok Yang. 2020. Towards verified stochastic variational inference for probabilistic programs. *Proc. ACM Program. Lang.* 4, POPL (2020), 16:1–16:33. https://doi.org/10.1145/3371084

[29] Carol Mak, C.-H. Luke Ong, Hugo Paquet, and Dominik Wagner. 2021. Densities of Almost Surely Terminating Probabilistic Programs are Differentiable Almost Everywhere. In *Programming Languages and Systems - 30th European Symposium on Programming, ESOP 2021, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2021, Luxembourg City, Luxembourg, March 27 - April 1, 2021, Proceedings (Lecture Notes in Computer Science, Vol. 12648)*, Nobuko Yoshida (Ed.). Springer, 432–461. https://doi.org/10.1007/978-3-030-72019-3_16

[30] Annabelle McIver and Carroll Morgan. 2004. Developing and Reasoning About Probabilistic Programs in *pGCL*. In *Refinement Techniques in Software Engineering, First Pernambuco Summer School on Software Engineering, PSSE 2004, Recife, Brazil, November 23-December 5, 2004, Revised Lectures (Lecture Notes in Computer Science, Vol. 3167)*, Ana Cavalcanti, Augusto Sampaio, and Jim Woodcock (Eds.). Springer, 123–155. https://doi.org/10.1007/11889229_4

[31] Annabelle McIver and Carroll Morgan. 2005. *Abstraction, Refinement and Proof for Probabilistic Systems*. Springer. https://doi.org/10.1007/b138392

[32] Sean P Meyn and Richard L Tweedie. 2012. *Markov chains and stochastic stability*. Springer Science & Business Media.

[33] Praveen Narayanan, Jacques Carette, Wren Romano, Chung-chieh Shan, and Robert Zinkov. 2016. Probabilistic Inference by Program Transformation in Hakaru (System Description). In *Functional and Logic Programming - 13th International Symposium, FLOPS 2016, Kochi, Japan, March 4-6, 2016, Proceedings (Lecture Notes in Computer Science, Vol. 9613)*, Oleg Kiselyov and Andy King (Eds.). Springer, 62–79. https://doi.org/10.1007/978-3-319-29604-3_5

[34] Van Chan Ngo, Quentin Carbonneaux, and Jan Hoffmann. 2018. Bounded expectations: resource analysis for probabilistic programs. In *Proceedings of the 39th ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI 2018, Philadelphia, PA, USA, June 18-22, 2018*, Jeffrey S. Foster and Dan Grossman (Eds.). ACM, 496–512. https://doi.org/10.1145/3192366.3192394

[35] David Pollard. 2002. *A user's guide to measure theoretic probability*. Number 8. Cambridge University Press.

[36] Mihai Putinar. 1993. Positive Polynomials on Compact Semi-algebraic Sets. *Indiana University Mathematics Journal* 42, 3 (1993), 969–984. http://www.jstor.org/stable/24897130

[37] RA Rankin. 1968. Real and Complex Analysis. By W. Rudin. Pp. 412. 84s. 1966. (McGraw-Hill, New York.). *The Mathematical Gazette* 52, 382 (1968), 412–412.

[38] Fredrik Ronquist, Jan Kudlicka, Viktor Senderov, Johannes Borgström, Nicolas Lartillot, Daniel Lundén, Lawrence Murray, Thomas B Schön, and David Broman. 2021. Universal probabilistic programming offers a powerful approach
to statistical phylogenetics. *Communications biology* 4, 1 (2021), 1–10.

[39] Reuven Y Rubinstein and Dirk P Kroese. 2016. *Simulation and the Monte Carlo method*. Vol. 10. John Wiley & Sons.

[40] Walter Rudin et al. 1976. *Principles of mathematical analysis*. Vol. 3. McGraw-hill New York.

[41] Feras A. Saad, Martin C. Rinard, and Vikash K. Mansinghka. 2021. SPPL: probabilistic programming with fast exact symbolic inference. In *PLDI '21: 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation, Virtual Event, Canada, June 20-25, 2021*, Stephen N. Freund and Eran Yahav (Eds.). ACM, 804–819. https://doi.org/10.1145/3453483.3454078

[42] Davide Sangiorgi. 2011. *Introduction to bisimulation and coinduction*. Cambridge University Press.

[43] Sriram Sankaranarayanan, Aleksandar Chakarov, and Sumit Gulwani. 2013. Static analysis for probabilistic programs: inferring whole program properties from finitely many paths. In *ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI '13, Seattle, WA, USA, June 16-19, 2013*, Hans-Juergen Boehm and Cormac Flanagan (Eds.). ACM, 447–458. https://doi.org/10.1145/2491956.2462179

[44] Sriram Sankaranarayanan, Henny B. Sipma, and Zohar Manna. 2004. Constraint-Based Linear-Relations Analysis. In *Static Analysis, 11th International Symposium, SAS 2004, Verona, Italy, August 26-28, 2004, Proceedings (Lecture Notes in Computer Science, Vol. 3148)*, Roberto Giacobazzi (Ed.). Springer, 53–68. https://doi.org/10.1007/978-3-540-27864-1_7

[45] Sam Staton, Hongseok Yang, Frank D. Wood, Chris Heunen, and Ohad Kammar. 2016. Semantics for probabilistic programming: higher-order functions, continuous distributions, and soft constraints. In *Proceedings of the 31st Annual ACM/IEEE Symposium on Logic in Computer Science, LICS '16, New York, NY, USA, July 5-8, 2016*, Martin Grohe, Eric Koskinen, and Natarajan Shankar (Eds.). ACM, 525–534. https://doi.org/10.1145/2933575.2935313

[46] Alfred Tarski. 1955. A lattice-theoretical fixpoint theorem and its applications. *Pacific journal of Mathematics* 5, 2 (1955), 285–309.

[47] David Tolpin, Jan-Willem van de Meent, and Frank D. Wood. 2015. Probabilistic Programming in Anglican. In *Machine Learning and Knowledge Discovery in Databases - European Conference, ECML PKDD 2015, Porto, Portugal, September 7-11, 2015, Proceedings, Part III (Lecture Notes in Computer Science, Vol. 9286)*, Albert Bifet, Michael May, Bianca Zadrozny, Ricard Gavaldà, Dino Pedreschi, Francesco Bonchi, Jaime S. Cardoso, and Myra Spiliopoulou (Eds.). Springer, 308–311. https://doi.org/10.1007/978-3-319-23461-8_36

[48] Jan-Willem van de Meent, Brooks Paige, Hongseok Yang, and Frank Wood. 2018. An Introduction to Probabilistic Programming. CoRR abs/1809.10756 (2018). arXiv:1809.10756

[49] Di Wang, Jan Hoffmann, and Thomas W. Reps. 2021. Central moment analysis for cost accumulators in probabilistic programs. In *PLDI '21: 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation, Virtual Event, Canada, June 20-25, 2021*, Stephen N. Freund and Eran Yahav (Eds.). ACM, 559–573. https://doi.org/10.1145/3453483.3454062

[50] Jinyi Wang, Yican Sun, Hongfei Fu, Krishnendu Chatterjee, and Amir Kafshdar Goharshady. 2021. Quantitative analysis of assertion violations in probabilistic programs. In *PLDI '21: 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation, Virtual Event, Canada, June 20-25, 2021*, Stephen N. Freund and Eran Yahav (Eds.). ACM, 1171–1186. https://doi.org/10.1145/3453483.3454102

[51] Peixin Wang. 2022. Tail-Bound Cost Analysis over Nondeterministic Probabilistic Programs. *Journal of Shanghai Jiaotong University (Science)* (2022), 1–11.

[52] Peixin Wang, Hongfei Fu, Krishnendu Chatterjee, Yuxin Deng, and Ming Xu. 2020. Proving expected sensitivity of probabilistic programs with randomized variable-dependent termination time. *Proc. ACM Program. Lang.* 4, POPL (2020), 25:1–25:30. https://doi.org/10.1145/3371093

[53] Peixin Wang, Hongfei Fu, Amir Kafshdar Goharshady, Krishnendu Chatterjee, Xudong Qin, and Wenjun Shi. 2019. Cost Analysis of Nondeterministic Probabilistic Programs. In *Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation (Phoenix, AZ, USA) (PLDI 2019)*. Association for Computing Machinery, New York, NY, USA, 204–220. https://doi.org/10.1145/3314221.3314581

[54] David Williams. 1991. *Probability with martingales*. Cambridge university press.
# A SUPPLEMENTARY MATERIAL FOR SECTION 2

## A.1 Basics of Probability Theory

A measurable space is a pair $(U, \Sigma_U)$, where $U$ is a nonempty set and $\Sigma_U$ is a $\sigma$-algebra on $U$, i.e., a family of subsets of $U$ such that $\Sigma_U \subseteq \mathcal{P}(U)$ contains $\emptyset$ and is closed under complementation and countable union. Elements of $\Sigma_U$ are called measurable sets. A function $f$ from a measurable space $(U, \Sigma_U)$ to another measurable space $(U_2, \Sigma_{U_2})$ is measurable if $f^{-1}(A) \in \Sigma_{U_2}$ for all $A \in \Sigma_{U_2}$.

A measure $\mu$ on a measurable space $(U, \Sigma_U)$ is a mapping from $\Sigma_U$ to $[0, \infty]$ such that (i) $\mu(\emptyset) = 0$ and (ii) $\mu$ is countably additive: for every pairwise-disjoint set sequence $\{A_n\}_{n \in \mathbb{N}}$ in $\Sigma_U$, it holds that $\mu(\bigcup_{n \in \mathbb{N}} A_n) = \sum_{n \in \mathbb{N}} \mu(A_n)$. We call the triple $(U, \Sigma_U, \mu)$ a measure space. If $\mu(U) = 1$, we call $\mu$ a probability measure, and $(U, \Sigma_U, \mu)$ a probability space. The Lebesgue measure $\lambda$ is the unique measure on $(\mathbb{R}, \Sigma_{\mathbb{R}})$ satisfying $\lambda([a, b)) = b - a$ for all valid intervals $[a, b)$ in $\Sigma_{\mathbb{R}}$. For each $n \in \mathbb{N}$, we have a measurable space $(\mathbb{R}^n, \Sigma_{\mathbb{R}^n})$ and a unique product measure $\lambda_n$ on $\mathbb{R}^n$ satisfying $\lambda_n(\prod_{i=1}^n A_i) = \prod_{i=1}^n \lambda(A_i)$ for all $A_i \in \Sigma_{\mathbb{R}^n}$.

The Lebesgue integral operator $\int$ is a partial operator that maps a measure $\mu$ on $(U, \Sigma_U)$ and a real-valued function $f$ on the same space $(U, \Sigma_U)$ to a real number or infinity, which is denoted by $\int f \, d\mu$ or $\int f(x)\mu(dx)$. The detailed definition of Lebesgue integral is somewhat technical, see [37, 40] for more details. Given a measurable set $A \in \Sigma_U$, the integral of $f$ over $A$ is defined by $\int_A f(x)\mu(dx) := \int f(x) \cdot [x \in A]\mu(dx)$ where $[-]$ is the Iverson bracket such that $[\phi] = 1$ if $\phi$ is true, and 0 otherwise. If $\mu$ is a probability measure, then we call the integral as the expectation of $f$, denoted by $\mathbb{E}_{X \sim \mu; A}[f]$, or $\mathbb{E}[f]$ when the scope is clear from the context.

For a measure $\nu$ on $(U, \Sigma_U)$, a measurable function $f: U \to \mathbb{R}_{\ge 0}$ is the density of $\nu$ with respect to $\mu$ if $\nu(A) = \int f(x) \cdot [x \in A]\mu(dx)$ for all measurable $A \in \Sigma_U$, and $\mu$ is called the reference measure (most often $\mu$ is the Lebesgue measure). Common families of probability distributions on the reals, e.g., uniform, normal distributions, are measures on $(\mathbb{R}, \Sigma_{\mathbb{R}})$. Most often these are defined in terms of probability density functions with respect to the Lebesgue measure. That is, for each $\mu_D$ there is a measurable function $\text{pdf}_D: \mathbb{R} \to \mathbb{R}_{\ge 0}$ that determines it: $\mu_D(A) := \int_A \text{pdf}_D(d\lambda)$. As we will see, density functions such as $\text{pdf}_D$ play an important role in Bayesian inference.

Given a probability space $(\Omega, \mathcal{F}, \mathcal{P})$, a random variable is an $\mathcal{F}$-measurable function $X: \Omega \to \mathbb{R} \cup \{+\infty, -\infty\}$. The expectation of a random variable $X$, denoted by $\mathbb{E}(X)$, is the Lebesgue integral of $X$ w.r.t. $\mathcal{P}$, i.e., $\int_X \mathcal{P} \, d\mathcal{P}$. A filtration of $(\Omega, \mathcal{F}, \mathcal{P})$ is an infinite sequence $\{\mathcal{F}_n\}_{n=0}^{\infty}$ such that for every $n \ge 0$, the triple $(\Omega, \mathcal{F}_n, \mathcal{P})$ is a probability space and $\mathcal{F}_n \subseteq \mathcal{F}_{n+1} \subseteq \mathcal{F}$. A stopping time w.r.t. $\{\mathcal{F}_n\}_{n=0}^{\infty}$ is a random variable $T: \Omega \to \mathbb{N} \cup \{0, \infty\}$ such that for every $n \ge 0$, the event $\{T \le n\}$ is in $\mathcal{F}_n$.

A discrete-time stochastic process is a sequence $\Gamma = \{X_n\}_{n=0}^{\infty}$ of random variables in $(\Omega, \mathcal{F}, \mathcal{P})$. The process $\Gamma$ is adapted to a filtration $\{\mathcal{F}_n\}_{n=0}^{\infty}$, if for all $n \ge 0$, $X_n$ is a random variable in $(\Omega, \mathcal{F}_n, \mathcal{P})$. A discrete-time stochastic process $\Gamma = \{X_n\}_{n=0}^{\infty}$ adapted to a filtration $\{\mathcal{F}_n\}_{n=0}^{\infty}$ is a martingale (resp. supermartingale, submartingale) if for all $n \ge 0$, $\mathbb{E}(|X_n|) < \infty$ and it holds almost surely (i.e., with probability 1) that $\mathbb{E}[X_{n+1} | \mathcal{F}_n] = X_n$ (resp. $\mathbb{E}[X_{n+1} | \mathcal{F}_n] \le X_n$, $\mathbb{E}[X_{n+1} | \mathcal{F}_n] \ge X_n$). See [54] for details. Applying martingales to qualitative and quantitative analysis of probabilistic programs is a well-studied technique [8, 10, 14].

## A.2 Details for WPTS Semantics

We denote by $\Lambda$ the set of all states, by $\Delta$ the set of all weighted states, and by $\Sigma_\Delta$ the product $\sigma$-algebra (on $\Delta$) among the discrete $\sigma$-algebra $(L, 2^L)$ for locations, the $\sigma$-algebra $\Sigma_{\mathbb{R}|\nu_P}$ for program valuations, and the $\sigma$-algebra $\Sigma_{\mathbb{R}}$ for the multiplicative likelihood weight. We define $\Sigma_\Delta^n$ (for $n \ge 1$) as the set $\{A_1 \times \cdots \times A_n \mid \forall 1 \le i \le n. (A_i \in \Sigma_\Delta)\}$, and $\Delta^\infty$ as the set of all infinite sequences of weighted states.
The probability space for the WPTS Π is defined such that its sample space is the set of all program runs, its σ-algebra is generated by the countable union $\bigcup_{n \ge 1} \{B \times \Lambda^{\infty} | B \in \Sigma_{\Lambda}^{n}\}$, and its probability measure $P$ is the unique one such that (i) $P(A \times \Lambda^{\infty}) = \mu_{\text{init}}(\{v | (\ell_{\text{init}}, v, 1) \in A\})$ for all $A \in \Sigma_{\Lambda}$, and (ii) $P(A \times B \times \Lambda^{\infty})$ (for every $A \in \Sigma_{\Lambda}$, $B \in \Sigma_{\Lambda}^{n}$ ($n \ge 1$)) equals the probability w.r.t the sampling of $\mu_{\text{init}}$ (for the initial program valuation) and $\mathcal{D}$ (for a sampling valuation in each step until the ($n+1$)-th step) that a program run $\{\Theta_n\}_{n \ge 0}$ is subject to $\Theta_0 \in A$ and $(\Theta_1, \dots, \Theta_{n+1}) \in B$. For each program valuation $v$, we denote by $P_v$ the probability measure of Π when the initial distribution is changed to the Dirac distribution at $v$.

## A.3 Sampling-based Semantics

We recall one prominent semantics in the literature, i.e., the sampling-based semantics [7, 45]. We show that the transition-based semantics in our work is equivalent to the widely-used sampling-based semantics in Bayesian statistical probabilistic programming.

The sampling-based semantics by Borgström et al. [7] interprets a probabilistic program as a deterministic program parameterized by a sequence of random draws sampled during the execution of the program.

A sampling trace is a finite sequence $t = \langle r_1, \dots, r_n \rangle$ of real numbers, and we define $\mathcal{T} := \bigcup_{n \in \mathbb{N}} \mathbb{R}^n$ as the set of all sampling traces. Given a probabilistic program $P$, a configuration $\sigma$ under the semantics is a tuple $\langle v, S, w, t \rangle$ where $v \in \text{Val}_{v_p}$, $S$ is the statement to be executed, $w \in [0, \infty)$ is the global weight variable whose value expresses how well the current computation matches the observations, and $t$ is a sampling trace. We denote by $\Sigma$ the set of all configurations.

The semantics operates on the configurations, where an execution of the program is initialized with $\sigma_0 = \langle v_0, P, 1, t \rangle$, and the termination configurations have the form of $\langle \_, \text{skip}, \_, [] \rangle$, for which _ is a “wildcard” character that matches everything and [] is an empty set. Fig. 10 shows the corresponding one-step reduction relation $\rightarrow$ (note that $\downarrow$ is the usual big-step semantics for deterministic Boolean and arithmetic expressions, so we omit it here).

Let $\rightarrow^*$ be the reflexive transitive closure of the one-step reduction $\rightarrow$ in Fig. 10. Given a probabilistic program $P$, we call a sampling trace $t$ terminating if $\langle v, P, 1, t \rangle \rightarrow^* \langle v', \text{skip}, w, [] \rangle$ for some valuations $v, v' \in \text{Val}_{v_p}$ and weight $w \in \mathbb{R}_{\ge 0}$, i.e., the program $P$ terminates under the samples drawn as in $t$.

Below we define the notion of posterior distributions by the sampling-based semantics. From the one-step reduction rules (in Fig. 10), we can reason about the global behavior of probabilistic programs in terms of the sampling traces they produce. That is, given a probabilistic program $P$, and a terminating trace $t$ such that $\langle v, P, 1, t \rangle \rightarrow^* \langle v', \text{skip}, w, [] \rangle$ for valuations $v, v' \in \text{Val}_{v_p}$ and weight $w \in \mathbb{R}_{\ge 0}$, we define the *value function* $\text{val}_P$ and the *weight function* $\text{wt}_P$ as follows:

$$ \text{val}_P(v, t) = \begin{cases} v' & \text{if } \langle v, P, 1, t \rangle \rightarrow^* \langle v', \text{skip}, w, [] \rangle \\ \text{unspecified} & \text{otherwise,} \end{cases} \quad (4) $$

$$ \text{wt}_P(v, t) = \begin{cases} w & \text{if } \langle v, P, 1, t \rangle \rightarrow^* \langle v', \text{skip}, w, [] \rangle \\ 0 & \text{otherwise.} \end{cases} \quad (5) $$

Moreover, we denote the return variable by $\text{val}_{P, \text{ret}}(v, t)$, i.e., $\text{val}_{P, \text{ret}}(v, t) := v'[\text{ret}]$. We also consider the measure space $(\mathcal{T}, \Sigma_{\mathcal{T}}, \mu_{\mathcal{T}})$ where $\mathcal{T} = \bigcup_{n \in \mathbb{N}} \mathbb{R}^n$ (as mentioned previously), $\Sigma_{\mathcal{T}} := \{\bigcup_{n \in \mathbb{N}} U_n \mid U_n \in \Sigma_{\mathbb{R}^n}\}$ and $\mu_{\mathcal{T}}(U) := \sum_{n \in \mathbb{N}} \lambda_n(U \cap \mathbb{R}^n)$. By definition, the measure space $(\mathcal{T}, \Sigma_{\mathcal{T}}, \mu_{\mathcal{T}})$ specifies the probability values for sets of sampling traces.
$$
\begin{array}{r@{\quad}l@{\quad}l}
\multicolumn{3}{c}{\langle \mathbf{v}, \text{skip}, w, t \rangle \to \langle \mathbf{v}, \text{skip}, w, t \rangle} \\
& \mathbf{v} \vdash E \Downarrow r & \\
& \langle \mathbf{v}, x := E, w, t \rangle \to \langle \mathbf{v}[x \mapsto r], \text{skip}, w, t \rangle \\
& \mathbf{v} \vdash B \Downarrow b, b = \text{true} & \\
& \langle \mathbf{v}, \text{if } B \text{ then } S_1 \text{ else } S_2 \text{ fi}, w, t \rangle \to \langle \mathbf{v}, S_1, w, t \rangle \\
& \mathbf{v} \vdash B \Downarrow b, b = \text{false} & \\
& \langle \mathbf{v}, \text{if } B \text{ then } S_1 \text{ else } S_2 \text{ fi}, w, t \rangle \to \langle \mathbf{v}, S_2, w, t \rangle \\
& w' = \text{pdf}_D(r) \ge 0 & \\
& \langle \mathbf{v}, x := \text{sample } D, w, r \:: \mathbf{t} \rangle \to \langle \mathbf{v}[x \mapsto r], \text{skip}, w \cdot w', \mathbf{t} \rangle \\
& \mathbf{v} \vdash B \Downarrow b, b = \text{true} & \\
& \langle \mathbf{v}, \text{while } B \text{ do } S \text{ od}, w, t \rangle \to \langle \mathbf{v}, S; \text{while } B \text{ do } S \text{ od}, w, t \rangle \\
& \mathbf{v} \vdash B \Downarrow b, b = \text{false} & \\
& \langle \mathbf{v}, \text{while } B \text{ do } S \text{ od}, w, t \rangle \to \langle \mathbf{v}, \text{skip}, w, t \rangle \\
& w' = \text{pdf}_D(x) \ge 0 & \\
& \langle \mathbf{v}, \text{observe}(x, D), w, t \rangle \to \langle \mathbf{v}, \text{skip}, w \cdot w', \mathbf{t} \rangle \\
& \langle \mathbf{v}, S_1, w, \mathbf{t} \rangle \to \langle \mathbf{v}', S'_1, w', \mathbf{t}' \rangle, S'_1 \ne \text{skip} & \\
& \langle \mathbf{v}, S_1; S_2, w, \mathbf{t} \rangle \to \langle \mathbf{v}', S'_1; S_2, w', \mathbf{t}' \rangle \\
& \langle \mathbf{v}, S_1, w, \mathbf{t} \rangle \to \langle \mathbf{v}', S'_1, w', \mathbf{t}' \rangle, S'_1 = \text{skip} & \\
& \langle \mathbf{v}, S_1; S_2, w, \mathbf{t} \rangle \to \langle \mathbf{v}', S_2, w', \mathbf{t}' \rangle & \\
& \langle \mathbf{v}, \text{return } x, w, \mathbf{t} \rangle \to \langle \mathbf{v}, \text{skip}, w, t \rangle &
\end{array}
$$

<caption>Fig. 10. One-step reduction for probabilistic programs.</caption>

*Posterior Distributions.* Given a probabilistic program *P*, an initial program valuation *v* ∈ Val<sub>*P*</sub> and a measurable set *U* ∈ Σ<sub>ℝ</sub>, we define the set of terminating traces where the value of the return variable falls into *U* as

$$
\mathcal{T}_{P,v,U} := \{ \boldsymbol{t} \in \mathcal{T} \mid \langle \mathbf{v}, P, 1, \boldsymbol{t} \rangle \to^* \langle \mathbf{v}', \text{skip}, w, \, \dots \rangle, \mathbf{v}'[\text{ret}] \in U \}
$$

and the set of all terminating traces as

$$
\mathcal{T}_{P,v} := \{ \boldsymbol{t} \in \mathcal{T} \mid \langle \mathbf{v}, P, 1, \boldsymbol{t} \rangle \to^* \langle \mathbf{v}', \text{skip}, w, \, \dots \rangle \}.
$$

Note that $\mathcal{T}_{P,v} = \mathcal{T}_{P,v,\mathbb{R}}$. Therefore, we can define the *unnormalised density* w.r.t *P*, *v*, *U* as

$$
\llbracket P \rrbracket_{v}(U) := \int_{\mathcal{T}_{P,v,U}} \mathrm{wt}_{P}(\mathbf{v}, t) \mu_{\mathcal{T}}(\mathrm{d}t). \qquad (6)
$$

That is, the integral takes all traces *t* on which *P* starts from **v** and evaluates to a value in *U*, weighting each *t* with the weight wt<sub>P</sub>(**v**, *t*) of the corresponding execution. The *normalising constant* is thus defined by

$$
Z_{P,v} := \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, t) \mu_{\mathcal{T}}(dt). \qquad (7)
$$

Therefore, the normalised *posterior distribution* is defined as $\text{posterior}_P(\mathbf{v}, U) := \frac{\llbracket P \rrbracket_{\mathbf{v}}(U)}{Z_{P,v}}$.

We call a program *P* integrable if its normalised constant is finite, i.e, 0 < Z<sub>P,v</sub> < ∞ for any
v ∈ Val<sub>P</sub>. Given an integrable program, we are interested in deriving lower and upper bounds on
the posterior distribution.
**Definition A.1 (Interval Bounds).** Given an integrable probabilistic program *P*, a program valuation *v* ∈ Val<sub>*P*</sub>, and a measurable set *U* ∈ Σ<sub>*R*</sub>, we call [*l*, *u*] an interval bound of posterior<sub>*P*</sub>(*v*, *U*) if *l* ≤ posterior<sub>*P*</sub>(*v*, *U*) ≤ *u* for two reals 0 ≤ *l* ≤ *u* ≤ 1.

## A.4 Equivalence between Posterior Distributions and Expected Weights

In this subsection, we introduce the equivalence between posterior distributions (in sample-based semantics) and expected weights (in transition-based semantics: WPTS).

**LEMMA A.2.** For all non-negative bounded measurable function *g* : ℝ<sup>|*V*<sub>*P*</sub>|</sup> → ℝ, and a program state Ξ = (ℓ, **v**), and a statement *P* := *S*<sub>1</sub>; *S*<sub>2</sub>, we have that

$$ 
\begin{aligned}
& \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
& = \int_{\mathcal{T}_{S_1,v}} \mathrm{wt}_{S_1}(\mathbf{v}, \mathbf{t}) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) ([\mathrm{val}_{S_1}(\mathbf{v}, \mathbf{t}) = \mathbf{v}']) \\
& \quad \cdot \int_{\mathcal{T}_{S_2,v'}} \mathrm{wt}_{S_2}(\mathbf{v}', \mathbf{t}') \cdot g(\mathrm{val}_{S_2}(\mathbf{v}', \mathbf{t}')) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}')
\end{aligned}
$$

**PROPOSITION A.3.** For all non-negative bounded measurable function *g* : ℝ<sup>|*V*<sub>*P*</sub>|</sup> → ℝ, a probabilistic program *P* and an initial program state Ξ = (ℓ, **v**) ∈ Λ, we have that

$$ \mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)] = \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}). $$

**PROOF.** We prove by induction on the structure of statements.

* Case *P* ≡ “skip”.

$$ 
\begin{aligned}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= \int \hat{w}_T(\omega) \cdot g(\hat{v}_T(\omega)) \mathbb{P}_{\Xi}(d\omega) \\
&= g(\mathbf{v}) \\
&= \int_{\mathcal{T}_{P,v}} [\mathbf{t} = []] \cdot g(\mathbf{v}) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t})
\end{aligned}
$$

* Case *P* ≡ “*x* := *E*”.

$$ 
\begin{aligned}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= \int \hat{w}_T(\omega) \cdot g(\hat{v}_T(\omega)) \mathbb{P}_{\Xi}(d\omega) \\
&= g(\mathbf{v}[x \mapsto [E](\mathbf{v})]) \\
&= \int_{\mathcal{T}_{P,v}} [\mathbf{t} = []] \cdot g(\mathbf{v}[x \mapsto [E](\mathbf{v})]) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t})
\end{aligned}
$$
• Case P ≡ “x := **sample D**”.

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= \int \hat{w}_T(\omega) \cdot g(\hat{v}_T(\omega)) \mathbb{P}_{\Xi}(d\omega) \\
&= \int g(v[x \mapsto r]) \mu_D(dr) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(v, t) \cdot g(\mathrm{val}_P(v, t)) \mu_{\mathcal{T}}(dt)
\end{align*}
$$

• Case P ≡ “**observe**(*x*, *D*)”.

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= \int \hat{w}_T(\omega) \cdot g(\hat{v}_T(\omega)) \mathbb{P}_{\Xi}(d\omega) \\
&= W_{\ell}(\mathbf{v}) \cdot g(\mathbf{v}) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(dt)
\end{align*}
$$

• Case P ≡ “**return x**”.

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= \int \hat{w}_T(\omega) \cdot g(\hat{v}_T(\omega)) \mathbb{P}_{\Xi}(d\omega) \\
&= g(\mathbf{v}) \\
&= \int_{\mathcal{T}_{P,v}} [\mathbf{t} = []] \cdot g(\mathbf{v}) \mu_{\mathcal{T}}(dt) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(dt)
\end{align*}
$$

• Case P ≡ “**if B then S₁ else S₂ fi**”. Assume the next state corresponding to the **then-branch** (resp. **else-branch**) is Ξ₁ = (ℓ₁, v) (resp. Ξ₂ = (ℓ₂, v)). Then we obtain that

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= [[B]](v) = true \cdot \mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)] + [[B]](v) = false \cdot \mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)] \\
&= [[B]](v) = true \cdot \int_{\mathcal{T}_{S_1,v}} \mathrm{wt}_{S_1}(v,t) \cdot g(\mathrm{val}_{S_1}(v,t)) \mu_{\mathcal{T}}(dt) \\
&\quad + [[B]](v) = false \cdot \int_{\mathcal{T}_{S_2,v}} \mathrm{wt}_{S_2}(v,t) \cdot g(\mathrm{val}_{S_2}(v,t)) \mu_{\mathcal{T}}(dt) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(v,t) \cdot g(\mathrm{val}_P(v,t)) \mu_{\mathcal{T}}(dt)
\end{align*}
$$
• Case P ≡ “S₁; S₂”.

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= \int \hat{w}_T(\omega) \mathbb{P}_{\Xi}(d\omega) \int [\omega_T = \Xi'] \cdot \hat{w}_T(\omega') \cdot g(\hat{v}_T(\omega')) \mathbb{P}_{\Xi'}(d\omega') \\
&= \int \hat{w}_T(\omega) \mathbb{P}_{\Xi}(d\omega) \left( [\omega_T = \Xi'] \cdot \int \cdot \hat{w}_T(\omega') \cdot g(\hat{v}_T(\omega')) \mathbb{P}_{\Xi'}(d\omega') \right) \\
&= \int \hat{w}_T(\omega) \mathbb{P}_{\Xi}(d\omega) \left( [\omega_T = \Xi'] \cdot \int_{\mathcal{T}_{S_2, v'}} \mathrm{wt}_{S_2}(\mathbf{v}', \mathbf{t}') \cdot g(\mathrm{val}_{S_2}(\mathbf{v}', \mathbf{t}')) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}') \right) \\
&= \int_{\mathcal{T}_{S_1, v}} \mathrm{wt}_{S_1}(\mathbf{v}, \mathbf{t}) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \left( [\mathrm{val}_{S_1}(\mathbf{v}, \mathbf{t}) = \mathbf{v}'] \cdot \int_{\mathcal{T}_{S_2, v'}} \mathrm{wt}_{S_2}(\mathbf{v}', \mathbf{t}') \cdot g(\mathrm{val}_{S_2}(\mathbf{v}', \mathbf{t}')) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}') \right) \\
&= \int_{\mathcal{T}_{P, v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t})
\end{align*}
$$

Here $\omega_i$ is the i-th element of the sequence $\omega = \{\Xi_n\}_{n \in \mathbb{N}}$, i.e., $\omega_i := \Xi_i$. $\omega_T$ is the last element of $\omega$, and $\Xi' = (\ell', \mathbf{v}')$. The third and fourth equalities follow from the induction hypothesis, and the last equality from Lemma A.2.

• Case P ≡ “while B do S od”. Assume the next state corresponding to the entry of the loop (resp. the exit of the loop) is Ξ₁ = (ℓ₁, v) (resp. Ξ₂ = (ℓ₂, v)). Then we obtain that

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot g(\hat{v}_T)]
&= [[B]](\mathbf{v}) = true \cdot \mathbb{E}_{\Xi_1} [\hat{w}_T \cdot g(\hat{v}_T)] + [[B]](\mathbf{v}) = false \cdot \mathbb{E}_{\Xi_2} [\hat{w}_T \cdot g(\hat{v}_T)] \\
&= [[B]](\mathbf{v}) = true \cdot \int \hat{w}_T(\omega) \cdot g(\hat{v}_T) \mathbb{P}_{\Xi_1}(d\omega) + [[B]](\mathbf{v}) = false \cdot g(\mathbf{v}) \\
&= [[B]](\mathbf{v}) = true \cdot \int_{\mathcal{T}_{S,P,v}} \mathrm{wt}_{S,P}(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_{S,P}(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
&\quad + [[B]](\mathbf{v}) = false \cdot \int_{\mathcal{T}_{skip,v}} [\mathbf{t} = []] \cdot g(\mathrm{val}_{skip}(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot g(\mathrm{val}_P(\mathbf{v}, \mathbf{t})) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t})
\end{align*}
$$

□

THEOREM A.4. Given a probabilistic program P, an initial program state Ξ = (ℓ, v) ∈ Λ and a measurable set U ∈ Σ_R, it holds that E_Ξ [ŵ_T · [v̂_T ∈ U]] = [P]_v(U). Moreover, the expected weight E_Ξ [ŵ_T] is equivalent to the normalising constant Z_P.

PROOF. We instantiate Proposition A.3 with $g(x) = [x \in \mathcal{U}]$. For any initial program state $\Xi = (\ell, \mathbf{v})$, we have that

$$
\begin{align*}
\mathbb{E}_v [\hat{w}_T \cdot [\hat{v}_T \in \mathcal{U}]]
&= \int_{\mathcal{T}_{P,v}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \cdot [\mathrm{val}_P(\mathbf{v}, \mathbf{t}) \in \mathcal{U}] \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
&= \int_{\mathcal{T}_{P,v,U}} \mathrm{wt}_P(\mathbf{v}, \mathbf{t}) \mu_{\mathcal{T}}(\mathrm{d}\mathbf{t}) \\
&= [P]_v(\mathcal{U})
\end{align*}
$$

We prove $\mathbb{E}_v [\hat{w}_T] = Z_P$ by setting $\mathcal{U} = \mathbb{R}$.

□
By Theorem A.4, we show the equivalence between posterior distributions under sampling-based semantics and expected weights under transition-based semantics. In the following, we will focus on developing approaches to infer interval bounds on expected weights.

## A.5 Proof for Proposition 2.8

**PROPOSITION A.5.** Given a WPTS in the form of (†), the interval-bound analysis of $Z_Π$ can be reduced to the interval-bound analysis of the expected weight $[\Pi]_v(\mathbb{R}^{|V_p|})$ for all $v \in \mathcal{V}$.

**PROOF.** Fix a WPTS $Π$ in the form of (†). By Definition 2.7, the normalising constant $Z_Π = [\Pi](\mathbb{R}^{|V_p|}) = \int_V [\Pi]_v(\mathbb{R}^{|V_p|}) \cdot \mu_{init}(dv)$ is the integral of the expected weight function $[\Pi]_v(\mathbb{R}^{|V_p|})$ over the finite set $\mathcal{V} = \text{supp}(\mu_{init})$. Assume there exist two integrable bound functions $Up, Lw$ such that $Lw(v) \le [\Pi]_v(\mathbb{R}^{|V_p|}) \le Up(v)$ for all $v \in \mathcal{V}$. Then the interval bounds $l_Z, u_Z$ for $Z_Π$ can be derived by the integrals of $Up, Lw$ over $\mathcal{V}$, that is,

$$l_Z := \int_V Lw(v) \cdot \mu_{init}(dv) \le [\Pi](\mathbb{R}^{|V_p|}) = \int_V [\Pi]_v(\mathbb{R}^{|V_p|}) \cdot \mu_{init}(dv) \le \int_V Up(v) \cdot \mu_{init}(dv) =: u_Z. \quad \square$$

**PROPOSITION A.6.** Given a WPTS in the form of (†) and any measurable set $\mathcal{U} \subset \mathbb{R}^{|V_p|}$, the interval-bound analysis of $[\Pi](\mathcal{U})$ can be reduced to the interval-bound analysis of the expected weight $[\Pi]_v(\mathbb{R}^{|V_p|})$ for all $v \in \mathcal{V}' \subseteq \mathcal{V}$.

**PROOF.** Fix a WPTS $Π$ in the form of (†) and a measurable set $\mathcal{U} \subset \mathbb{R}^{|V_p|}$. By Definition 2.7, $[\Pi](\mathcal{U}) = \int_V [\Pi]_v(\mathcal{U}) \cdot \mu_{init}(dv)$. By Definition 2.6, $[\Pi]_v(\mathcal{U}) = \mathbb{E}_v [[\hat{v}_T \in \mathcal{U}] \cdot \hat{w}_T]$ for any initial program valuation $v$. As $\hat{v}_T$ is the random (vector) variable of the program valuation at termination and $\mathcal{U} \ne \mathbb{R}^{|V_p|}$, it is not possible to track the valuation of $\hat{v}_T$ for each program run starting from $v$, which makes it challenging to compute interval bounds for $[\Pi]_v(\mathcal{U})$ directly. However, in practice, we only care about the posterior distribution that whether the return variable falls into some target sets of interest.<fnref n="10" /> That is, the set $\mathcal{U}$ is defined such that $[\hat{v}_T \in \mathcal{U}] = 1 \text{ iff } \hat{v}_T[\text{ret}] \in \mathcal{U}(\text{ret})$ where $\mathcal{U}(\text{ret})$ is the element in $\mathcal{U}$ that corresponds to the set w.r.t. the return variable. Below we distinguish the computation into two cases.

Case 1. If the value of $ret$ is determined once and unaffected by loop iterations, then $[\hat{v}_T[\text{ret}] \in \mathcal{U}(\text{ret})]$ is equivalent to $[v_{init}[\text{ret}] \in \mathcal{U}(\text{ret})]$. We can conclude that $[\Pi]_{v_{init}}(\mathcal{U}) = \mathbb{E}_{v_{init}} [[\hat{v}_T \in \mathcal{U}] \cdot \hat{w}_T] = \mathbb{E}_{v_{init}} [\hat{w}_T] = [\Pi]_{v_{init}}(\mathbb{R}^{|V_p|})$ if $v_{init}[\text{ret}] \in \mathcal{U}(\text{ret})$, and 0 otherwise.

Case 2. If the value of $ret$ is affected during loop iterations, then we can construct a WPTS $Π_{qL}$ by adding a conditional branch of the form “if $v_T \notin \mathcal{U}$ then score(0) fi” immediately before the termination of $Π$. And we can prove that $[\Pi]_{v_{init}}(\mathcal{U})$ is equivalent to $[\Pi_{qL}]_{v_{init}}(\mathbb{R}^{|V_p|})$ for all $v_{init} \in \mathcal{V}$.

The correctness of Case 1 above is straightforward to be proved, while the correctness of Case 2 follows from the proposition below. $\square$

**PROPOSITION A.7.** Given a WPTS $Π$ in the form of (†) and a measurable set $\mathcal{U} \subset \mathbb{R}^{|V_p|}$, for any initial program valuation $v_{init} \in \mathcal{V}$, the expected weight $[\Pi]_{v_{init}}(\mathcal{U})$ is equivalent to $[\Pi_{qL}]_{v_{init}}(\mathbb{R}^{|V_p|})$ where $\Pi_{qL}$ is constructed by the method described in Case 2 above.

**PROOF.** Fix a WPTS $Π$ in the form of (†), a measurable set $\mathcal{U} \subset \mathbb{R}^{|V_p|}$ and a new WPTS $Π_{qL}$ constructed by the method described in Case 2. Let the termination location of $Π$ be $\ell_{out}$ and the termination location of $\Pi_{qL}$ be $\ell'_{out}$. Define a predicate $\Phi_{qL}$ over $V_p$ such that $\mathcal{U} = \{v | v \models \Phi_{qL}\}$.

<footnote n="10">¹⁰ For the sake of simplicity, we consider that the WPTS has only one return variable $ret \in V_p$ of interest, but it can be straightforwardly extended to general cases.</footnote>
$\Phi_{\mathcal{U}}$. Then $\Pi_{\mathcal{U}}$ has two additional transitions than $\Pi$, i.e., $\langle \ell_{\text{out}}, \Phi_{\mathcal{U}}, F_{\tau} \rangle$ with $F_{\tau} := \langle \ell'_{\text{out}}, 1, \text{id}, \bar{0} \rangle$ where $\text{id}$ is the identity function and $\bar{0}$ is the constant function that always takes the value 0, and $\langle \ell_{\text{out}}, \neg\Phi_{\mathcal{U}}, F_{\perp} \rangle$ with $F_{\perp} := \langle \ell'_{\text{out}}, 1, \text{id}, \bar{1} \rangle$ where $\bar{1}$ is the constant function that always takes the value 1. Let the termination time of $\Pi$ be $T$ and the termination time of $\Pi_{\mathcal{U}}$ be $T'$. Therefore, for the random variable of the multiplicative likelihood weight at termination, we have that $\hat{w}_{T'} = [\hat{v}_{T'} \in \mathcal{U}] \cdot \hat{w}_{T} + [\hat{v}_{T'} \notin \mathcal{U}] \cdot 0$. By Definition 2.6,

$$ \begin{align*} [[\Pi]]_{v_{\text{init}}} (\mathcal{U}) &= \mathbb{E}_{v_{\text{init}}} [[\hat{v}_T \in \mathcal{U}] \cdot \hat{w}_T] \\ &= \mathbb{E}_{v_{\text{init}}} [[\hat{v}_T \in \mathcal{U}] \cdot \hat{w}_T + [\hat{v}_T \notin \mathcal{U}] \cdot 0] \\ &= \mathbb{E}_{v_{\text{init}}} [[\hat{v}_{T'} \in \mathcal{U}] \cdot \hat{w}_T + [\hat{v}_{T'} \notin \mathcal{U}] \cdot 0] \\ &= \mathbb{E}_{v_{\text{init}}} [\hat{w}_{T'}] \\ &= \mathbb{E}_{v_{\text{init}}} [ [\hat{v}_{T'} \in \mathbb{R}^{|V_{\text{P}}|}] \cdot \hat{w}_{T'} ] \\ &= [[\Pi_{\mathcal{U}}]]_{v_{\text{init}}} (\mathbb{R}^{|V_{\text{P}}|}) \end{align*} $$

where the third equality is derived from the fact that $\hat{v}_{T'} = \hat{v}_T$ and the fourth equality is obtained from the fact that $[\hat{v}_{T'} \in \mathbb{R}^{|V_{\text{P}}|}] = 1$. $\square$

Based on Proposition A.5 and Proposition A.6, in the following (Section 4) we develop approaches to obtain interval bounds for expected weights in the form $[[\Pi]]_v(\mathbb{R}^{|V_{\text{P}}|})$ where $v$ is an initial program valuation.

# B SUPPLEMENTARY MATERIAL FOR SECTION 4

## B.1 Basics of Fixed Point Theory

Given a partial order $\sqsubseteq$ on a set $K$ and a subset $K' \subseteq K$, an *upper bound* of $K'$ is an element $u \in K$ that is no smaller than every element of $K'$, i.e., $\forall k' \in K'$. $k' \sqsubseteq u$. Similarly, a *lower bound* for $K'$ is an element $l$ that is no greater than every element of $K'$, i.e. $\forall k' \in K'$. $l \sqsubseteq k'$. The *supremum* of $K'$, denoted by $\sqsupseteq K'$, is an element $u^* \in K$ such that $u^*$ is an upper-bound of $K'$ and for every upper bound $u$ of $K'$, we have $u^* \sqsubseteq u$. Similarly, the *infimum* $\sqinfloor K'$ is a lower bound $l^*$ of $K'$ such that for every lower-bound $l$ of $K'$, we have $l \sqsubseteq l^*$. We define $\perp := \sqinfloor K$ and $\top := \sqsupseteq K$. In general, suprema and infima may not exist.

A partial order $(K, \sqsubseteq)$ is called a *complete lattice* if every subset $K' \subseteq K$ has a supremum and an infimum. Given a partial order $(K, \sqsubseteq)$, a function $f: K \to K$ is called *monotone* if for every $k_1 \sqsubseteq k_2$ in $K$, we have $f(k_1) \sqsubseteq f(k_2)$.

Given a complete lattice $(K, \sqsubseteq)$, a function $f: K \to K$ is called *continuous* if for every increasing chain $k_0 \sqsubseteq k_1 \sqsubseteq \dots$ in $K$, we have $f(\sqcup\{k_n\}_{n=0}^{\infty}) = \sqcup\{f(k_n)\}_{n=0}^{\infty}$, and *cocontinuous* if for every decreasing chain $k_0 \supseteq k_1 \supseteq \dots$ of elements of $K$, we have $f(\sqprod\{k_n\}_{n=0}^{\infty}) = \sqprod\{f(k_n)\}_{n=0}^{\infty}$.

Given a complete lattice $(K, \sqsubseteq)$ and a function $f: K \to K$, an element $k \in K$ is called a *fixed-point* if $f(k) = k$. Moreover, $k$ is a *pre fixed-point* if $f(k) \sqsubseteq k$ and a *post fixed-point* if $k \sqsubseteq f(k)$. The *least fixed-point* of $f$, denoted by $lfpf$, is the fixed-point that is no greater than every fixed-point under $\sqsubseteq$. Analogously, the *greatest fixed-point* of $f$, denoted by $gfpf$, is the fixed-point that is no smaller than all fixed-points.

**THEOREM B.1 (KLEENE [42])**. Let $(K, \sqsubseteq)$ be a complete lattice and $f: K \to K$ be an continuous function. Then, we have

$$ lfp f = \sqcup_{i \ge 0} \{f^{(i)}(\perp)\}. $$

Analogously, if $f$ is cocontinuous, then we have

$$ gfp f = \sqprod_{i \ge 0} \{f^{(i)}(\top)\}. $$
## B.2 Proofs for Our Fixed-Point Approach

**Theorem 4.5.** The expected-weight function $ew$ is the least fixed point of the expected-weight transformer $ewt$.

**PROOF.** Define the step-bounded weight random variable $\hat{w}^n(\ell, v)$ starting from any program state $\Xi = (\ell, v)$ for a step bound $n \in \mathbb{N}$ by

$$ \hat{w}^n(\ell, v)(\omega) = \begin{cases} \hat{w}(\ell, v)(\omega) & \text{if } T(\omega) \le n \\ 0 & \text{otherwise} \end{cases}. $$

Since we always assume that the underlying WPTS is almost-surely terminating, it follows that the sequence of random variables $\{\hat{w}^n\}_{n \in \mathbb{N}}$ converges non-decreasingly to $W$.

Given any program state $\Xi = (\ell, v)$ with a unique transition $\tau = (\ell, \phi_\tau, f_\tau)$ satisfying $v \models \phi_\tau$, define the step-bounded expected-weight function $ew^n$ by $ew^n(\ell, v) = E_{(\ell, v)}[\hat{w}^n]$. Without loss of generality, we assume there is only one fork $f_\tau$ in this transition. Assume the next sampling valuation from $\Xi$ is $r_0$ and the next program state is $\Xi' = (\ell', v')$, i.e., $v' = f_\tau(v, r_0)$. Following the symbols in Section 2, we denote the probability space of the WPTS $\Pi$ starting from $(\ell, v)$, i.e., the program runs starting from $\Xi = (\ell, v)$ as $(\Omega, \mathcal{F}, \mathbb{P})_\Xi$. By Tonelli-Fubini Theorem, we have that for all $n \ge 0$,

$$ \begin{aligned} ew^{n+1}(\ell, v) &= \int \hat{w}^{n+1}(\ell, v) d\mathbb{P}_\Xi \\ &= \int \hat{w}^{n+1}(\ell, v) d(D_{r_0} \times \mathbb{P}_{\Xi'}) \\ &= \int wt \cdot W^n_{(\ell', v')}(\omega) d(D_{r_0} \times \mathbb{P}_{\Xi'}) \\ &= \int_{r_0} \int_{\omega_{\Xi'}} wt \cdot W^n_{(\ell', v')}(\omega) d\mathbb{P}_{\Xi'} dD_{r_0} \\ &= \int_{r_0} wt \cdot \left( \int_{\omega_{\Xi'}} \hat{w}^n_{(\ell', v')}(\omega) d\mathbb{P}_{\Xi'} \right) dD_{r_0} \\ &= \int_{r_0} wt \cdot ew^n(\ell', v') dD_{r_0} \\ &= E_{r_0} [wt \cdot ew^n(\ell', v')] \\ &= ewt(ew^n)(\ell, v) \end{aligned} $$

By applying MCT to the both sides of the equality above, we have that

$$ ew(\ell, v) = ewt(ew)(\ell, v). $$

This shows that $ew$ is a fixed point of $ewt$. Furthermore, given any fixed point $h$ of $ewt$, since (i) $ew^0 \le h$ and (ii) $ew^n \le h$ implies $ew^{n+1} = ewt(ew^n) \le ewt(h) = h$, one can prove by a straightforward induction on $n$ that $ew^n \le h$ for all $n \ge 0$. It follows from $ew = \lim_{n \to \infty} ew^n$ that $ew$ is the least fixed point of $ewt$. $\square$

In order to show the uniqueness of the fixed point, we first prove that $ewt$ is both continuous and cocontinuous.

**PROPOSITION B.2.** If $M \in [0, \infty)$, then the expected-weight transformer $ewt: K_M \to K_M$ is both continuous and cocontinuous.

**PROOF.** We first prove that $ewt$ is well-defined. Given an arbitrary $h \in K_M$, for any $\Xi = (\ell, v) \in \Lambda$,
* When $\ell = \ell_{\text{out}}$, $\text{ewt}(h)(\ell, \mathbf{v}) = 1$.

* When $\ell \neq \ell_{\text{out}}$, for a unique transition $\tau = (\ell, \phi_{\tau}, f_{\tau}, \ell')$ such that $\mathbf{v} \models \phi_{\tau}$,

$$ 
\begin{aligned}
\text{ewt}(h)(\ell, \mathbf{v}) &= \mathbb{E}_{\mathbf{r}} [h(\ell', f_{\tau}(\mathbf{v}, \mathbf{r})) \cdot W(\ell, \mathbf{v})] \\
&\leq M \cdot \text{maxscore} \\
&< \infty
\end{aligned}
$$

where maxscore is the maximum of W given any state $\Xi$. As W is a non-negative function, we can prove that $\text{ewt}(h)(\ell, \mathbf{v}) \geq 0$. Thus, $\text{ewt}$ is well defined. Next, we prove that $\text{ewt}$ is monotone. Given any two functions $h_1, h_2 \in \mathcal{K}_M$ such that $h_1 \leq h_2$, by case analysis on $(\ell, \mathbf{v})$,

* If $\ell = \ell_{\text{out}}$, $\text{ewt}(h_1)(\ell, \mathbf{v}) = 1 = \text{ewt}(h_2)(\ell, \mathbf{v})$.

* If $\ell \neq \ell_{\text{out}}$, given a unique transition $\tau = (\ell, \phi_{\tau}, f_{\tau}, \ell')$ such that $\mathbf{v} \models \phi_{\tau}$,

$$ 
\begin{aligned}
\text{ewt}(h_1)(\ell, \mathbf{v}) &= \mathbb{E}_{\mathbf{r}} [h_1(\ell', f_{\tau}(\mathbf{v}, \mathbf{r})) \cdot W(\ell, \mathbf{v})] \\
&\leq \mathbb{E}_{\mathbf{r}} [h_2(\ell', f_{\tau}(\mathbf{v}, \mathbf{r})) \cdot W(\ell, \mathbf{v})] \\
&= \text{ewt}(h_2)(\ell, \mathbf{v})
\end{aligned}
$$

Therefore, $\text{ewt}(h_1) \leq \text{ewt}(h_2)$, hence it is monotone. Then we prove upper continuity of $\text{ewt}$. Choose any increasing chain $h_0 \sqsubseteq h_1 \sqsubseteq h_2 \sqsubseteq \dots$ and do another case analysis on $(\ell, \mathbf{v})$:

* If $\ell = \ell_{\text{out}}$, then

$$ \text{ewt}(\bigcup_{n \geq 0} \{h_n\})(\ell, \mathbf{v}) = 1 = \bigcup_{n \geq 0} \{\text{ewt}(h_n)\}(\ell, \mathbf{v}). $$

* Otherwise, for a unique transition $\tau = (\ell, \phi_{\tau}, f_{\tau})$ such that $\mathbf{v} \models \phi_{\tau}$:

$$ 
\begin{aligned}
\text{ewt}(\bigcup_{n \geq 0} \{h_n\})(\ell, \mathbf{v}) &= \mathbb{E}_{\mathbf{r}} \left[ \text{wt}_j(\mathbf{v}, \mathbf{r}) \cdot \left( \bigcup_{n \geq 0} \{h_n\} \right) (\ell', f_{\tau}(\mathbf{v}, \mathbf{r})) \right] \\
&= \mathbb{E}_{\mathbf{r}} \left[ \sup_{n \geq 0} \{h_n(\ell', f_{\tau}(\mathbf{v}, \mathbf{r})) \} \right] \\
&= \mathbb{E}_{\mathbf{r}} \left[ \lim_{n \to \infty} \{h_n(\ell', f_{\tau}(\mathbf{v}, \mathbf{r})) \} \right] \\
\stackrel{\text{MCT}}{=} & \lim_{n \to \infty} \mathbb{E}_{\mathbf{r}} [h_n(\ell', f_{\tau}(\mathbf{v}, \mathbf{r}))] \\
&= \lim_{n \to \infty} \text{ewt}(h_n)(\ell, \mathbf{v}) \\
&= \sup_{n \geq 0} \{\text{ewt}(h_n)(\ell, \mathbf{v})\} \\
&= \bigcup_{n \geq 0} \{\text{ewt}(h_n)\} (\ell, \mathbf{v})
\end{aligned}
$$

The “MCT” above denotes the monotone convergence theorem. A similar argument establishes cocontinuity for integrable $h_0$ and decreasing chains. □

Then the uniqueness follows from Theorem B.1.

**Theorem 4.5.** Let $\Pi$ be a non-score-recursive WPTS whose weights are bounded in $[-M, M]$ for a finite $M \geq 1$. Then the expected-weight function $ew$ is the least fixed point of the expected-weight
transformer ewt in the complete lattice ($\mathcal{K}_M$, $\le$). Furthermore, if the WPTS $\Pi$ is AST, then the function ew is the unique fixed point of the higher-order function ewt in ($\mathcal{K}_M$, $\le$) when $M \ge 1$.

PROOF. The proof follows similar arguments in [50, Theorem 4.4]. By Proposition B.2, we have that for every state $\Xi = (\ell, \mathbf{v})$,

* lfp ewt($\ell, \mathbf{v}$) = $\lim_{n \to \infty}$ ewt$^n(\perp)(\ell, \mathbf{v})$, and

* gfp ewt($\ell, \mathbf{v}$) = $\lim_{n \to \infty}$ ewt$^n(\top)(\ell, \mathbf{v})$.

By the definition of ewt$^n_M$ and Proposition B.2, we have that

* ewt$^n(\perp)(\ell, \mathbf{v}) = \mathbb{E}_{\Xi}[W \cdot [T \le n]] - M \cdot \mathbb{P}(T > n)$, and

* ewt$^n(\top)(\ell, \mathbf{v}) = \mathbb{E}_{\Xi}[W \cdot [T \le n]] + M \cdot \mathbb{P}(T > n)$.

Recall that we assume the underlying PTS to be almost-surely terminating. Hence, $\lim_{n \to \infty} \mathbb{P}(T > n) = \mathbb{P}(T = \infty) = 0$. It follows that lfp ewt($\ell, \mathbf{v}$) = gfp ewt($\ell, \mathbf{v}$), i.e., the fixed point is unique. $\square$

## B.3 Proof for the OST Variant

Theorem 4.8 (OST Variant) Let {$X_n$}$_{n=0}^\infty$ be a supermartingale adapted to a filtration $\mathcal{F} = \{\mathcal{F}_n\}^\infty_{n=0}$, and $\kappa$ be a stopping time w.r.t. the filtration $\mathcal{F}$. Then the following condition (U) is sufficient to ensure that $\mathbb{E}(|X_\kappa|) < \infty$ and $\mathbb{E}(X_\kappa) \le \mathbb{E}(X_0)$:

(U) There exist integers $b_1, b_2 > 0$ and real numbers $c_1 > 0, c_2 > c_3 > 0$ such that (i) $\mathbb{P}(\kappa > n) \le c_1 \cdot e^{-c_2 \cdot n}$ for sufficiently large $n \in \mathbb{N}$, and (ii) for all $n \in \mathbb{N}$, $|X_{n+1} - X_n| \le b_1 \cdot n^{b_2} \cdot e^{c_3 \cdot n}$ holds almost surely.

PROOF. For every $n \in \mathbb{N}_0$,

$$ 
\begin{aligned}
|X_{\kappa \wedge n}| &= \left| X_0 + \sum_{k=0}^{\kappa \wedge n-1} (X_{k+1} - X_k) \right| \\
&= \left| X_0 + \sum_{k=0}^{\infty} (X_{k+1} - X_k) \cdot \mathbf{1}_{\kappa > k \wedge n > k} \right| \\
&\le |X_0| + \sum_{k=0}^{\infty} |(X_{k+1} - X_k) \cdot \mathbf{1}_{\kappa > k \wedge n > k}| \\
&\le |X_0| + \sum_{k=0}^{\infty} |(X_{k+1} - X_k) \cdot \mathbf{1}_{\kappa > k}| .
\end{aligned}
$$
Then

$$
\begin{align*}
& \mathbb{E}\left( |X_0| + \sum_{k=0}^{\infty} |(X_{k+1} - X_k) \cdot \mathbf{1}_{\kappa>k}| \right) \\
&= \text{(By Monotone Convergence Theorem)} \\
& \mathbb{E}(|X_0|) + \sum_{k=0}^{\infty} \mathbb{E}(|(X_{k+1} - X_k) \cdot \mathbf{1}_{\kappa>k}|) \\
&= \mathbb{E}(|X_0|) + \sum_{k=0}^{\infty} \mathbb{E}(|X_{k+1} - X_k| \cdot \mathbf{1}_{\kappa>k}) \\
&\le \mathbb{E}(|X_0|) + \sum_{k=0}^{\infty} \mathbb{E}\left( \lambda \cdot k^d \cdot e^{c_3 \cdot k} \cdot \mathbf{1}_{\kappa>k} \right) \\
&= \mathbb{E}(|X_0|) + \sum_{k=0}^{\infty} \lambda \cdot k^d \cdot e^{c_3 \cdot k} \cdot \mathbb{P}(\kappa > k) \\
&\le \mathbb{E}(|X_0|) + \sum_{k=0}^{\infty} \lambda \cdot k^d \cdot e^{c_3 \cdot k} \cdot c_1 \cdot e^{-c_2 \cdot k} \\
&= \mathbb{E}(|X_0|) + \lambda \cdot c_1 \cdot \sum_{k=0}^{\infty} k^d \cdot e^{-(c_2-c_3) \cdot k} \\
&< \infty .
\end{align*}
$$

where the first inequality is obtained by Condition (ii), and the second inequality is derived from Condition (i).

Therefore, by Dominated Convergence Theorem and the fact that $X_\kappa = \lim_{n \to \infty} X_{\kappa \wedge n}$ a.s.,

$$
\mathbb{E}(X_{\kappa}) = \mathbb{E}\left(\lim_{n \to \infty} X_{\kappa \wedge n}\right) = \lim_{n \to \infty} \mathbb{E}(X_{\kappa \wedge n}) .
$$

Finally, the result follows from properties for the stopped process {$X_{\kappa \wedge n}$}$_{n \in \mathbb{N}_0}$ that

$$
\mathbb{E}(X_{\kappa}) \leq \mathbb{E}(X_0) .
$$

□

B.4 Proofs for Our OST-Based Approach

**Theorem 4.9.** [Upper Bounds on Expected Weights] Let Π be a score-recursive WPTS that has the bounded update property. Suppose that there exist real numbers $c_1 > 0$ and $c_2 > c_3 > 0$ such that

• $\mathbb{P}(\kappa > n) \le c_1 \cdot e^{-c_2 \cdot n}$ for sufficiently large $n \in \mathbb{N}$, and

• for each score function $wt$ in $\Pi$ we have $|wt| \le e^{c_3}$.

Then for any PUWF $h$ over $\Pi$, we have that $[\Pi]_{v_{\text{init}}}(\mathbb{R}^{|V_{\text{P}}|}) \le h(\ell_{\text{init}}, v_{\text{init}})$ for any initial state $(\ell_{\text{init}}, v_{\text{init}})$.

**PROOF.** Define the stochastic process {$X_n$}$_{n=0}^\infty$ as $X_n := h(\ell_n, v_n)$ where $(\ell_n, v_n)$ is the program state at the $n$-th step of a program run. Then construct a stochastic process {$Y_n$}$_{n=0}^\infty$ such that $Y_n := X_n \cdot \prod_{i=0}^{n-1} W_i$ where $W_i$ is the weight at the $i$-th step of the program run. According to Condition (C1), we have that $\mathbb{E}[X_{n+1} \cdot W_n | \mathcal{F}_n] \le X_n$. Therefore, by the “take out what is known” property of conditional expectation (see [54]), it follows that
$$ \mathbb{E}\left[X_{n+1} \cdot \prod_{i=0}^{n} W_i | \mathcal{F}_n\right] \leq X_n \cdot \prod_{i=0}^{n-1} W_i \\ \Leftrightarrow \mathbb{E}[Y_{n+1} | \mathcal{F}_n] \leq \mathbb{E}[Y_n], $$

which means that $\mathbb{E}[Y_{n+1}] \leq \mathbb{E}[Y_n]$ from the basic property of conditional expectation. By an easy induction on $n$, we have that $\mathbb{E}[Y_n] \leq \mathbb{E}[Y_0] < \infty$ for all $n \geq 0$, thus the conditional expectation is also taken in the normal sense as each $Y_n$ is indeed integrable. Hence, $\{Y_n\}_{n=0}^{\infty}$ is a supermartingale. Moreover, we have from the bounded-update property that $|X_{n+1}| \leq \zeta \cdot (n+1)^d$ for a real number $\zeta > 0$. By definition, we obtain that for sufficiently large $n$,

$$ \begin{aligned} |Y_{n+1} - Y_n| &= \left| X_{n+1} \cdot \prod_{i=0}^{n} W_i - X_n \cdot \prod_{i=0}^{n-1} W_i \right| \\ &\leq \left| X_{n+1} \cdot \prod_{i=0}^{n} W_i \right| + \left| X_n \cdot \prod_{i=0}^{n-1} W_i \right| \\ &< e^{c_3 \cdot n} \cdot (|X_{n+1}| + |X_n|) \\ &\leq e^{c_3 \cdot n} \cdot [\zeta \cdot (n+1)^d + \zeta \cdot n^d] \\ &\leq \lambda \cdot n^d \cdot e^{c_3 \cdot n} \end{aligned} $$

where the first inequality is induced by the triangle inequality, and the second inequality is derived from the bounded stepwise weight condition such that each $W_i \in [0, e^{c_3}]$ and the fact $W_0 = 1$. By applying the OST variant (Theorem 4.8), we obtain that $\mathbb{E}[Y_T] \leq \mathbb{E}[Y_0]$. By definition and Condition (C2),

$$ \begin{aligned} Y_T &= h(\ell_T, \mathbf{v}_T) \cdot \prod_{i=0}^{T-1} W_i \\ &= h(\ell_{\text{out}}, \mathbf{v}_T) \cdot \prod_{i=0}^{T-1} W_i \\ &= \prod_{i=0}^{T-1} W_i \end{aligned} $$

Finally, we have that $\mathbb{E}_{\Xi}[W_{\infty}] = \mathbb{E}[\prod_{i=0}^{T-1} W_i] \leq \mathbb{E}[Y_0] = h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$. $\square$

**Theorem 4.9.** [Lower Bounds on Expected Weights] Let $\Pi$ be a score-recursive WPTS that has the bounded update property. Suppose that there exist real numbers $c_1 > 0$ and $c_2 > c_3 > 0$ such that

* $\mathbb{P}(\kappa > n) \leq c_1 \cdot e^{-c_2 \cdot n}$ for sufficiently large $n \in \mathbb{N}$, and

* for each score function $wt$ in $\Pi$ we have $|wt| \leq e^{c_3}$.

Then for any PLWF $h$ over $\Pi$, we have that $\|\Pi\|_{\mathbf{v}_{\text{init}}}(\mathbb{R}^{|V_{\mathbf{p}}|}) \geq h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$ for any initial state $(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$.

**PROOF.** Define the stochastic process $\{X_n\}_{n=0}^{\infty}$ as $X_n := h(\ell_n, \mathbf{v}_n)$ where $(\ell_n, \mathbf{v}_n)$ is the program state at the $n$-th step of a program run. Then construct a stochastic process $\{-Y_n\}_{n=0}^{\infty}$ such that $-Y_n := -X_n \cdot \prod_{i=0}^{n-1} W_i$ where $W_i$ is the weight at the $i$-th step of the program run. According to Condition (C1'), we have that $\mathbb{E}[-X_{n+1} \cdot W_n | \mathcal{F}_n] \leq -X_n$. Therefore, by the “take out what is known” property of conditional expectation (see [54]), it follows that
$$ \mathbb{E}\left[-X_{n+1} \cdot \prod_{i=0}^{n} W_i | \mathcal{F}_n\right] \leq -X_n \cdot \prod_{i=0}^{n-1} W_i \\ \Leftrightarrow \mathbb{E}[-Y_{n+1} | \mathcal{F}_n] \leq \mathbb{E}[-Y_n], $$

which means that $\mathbb{E}[-Y_{n+1}] \leq \mathbb{E}[-Y_n]$ from the basic property of conditional expectation. By an easy induction on $n$, we have that $\mathbb{E}[-Y_n] \leq \mathbb{E}[-Y_0] < \infty$ for all $n \geq 0$, thus the conditional expectation is also taken in the normal sense as each $Y_n$ is indeed integrable. Hence, $\{-Y_n\}_{n=0}^{\infty}$ is a supermartingale. Moreover, we have from the bounded-update property that $|X_{n+1}| \leq \zeta \cdot (n+1)^d$ for a real number $\zeta > 0$. By definition, we obtain that for sufficiently large $n$,

$$ \begin{aligned} |-Y_{n+1} - (-Y_n)| &= \left| X_{n+1} \cdot \prod_{i=0}^{n} W_i - X_n \cdot \prod_{i=0}^{n-1} W_i \right| \\ &\leq \left| X_{n+1} \cdot \prod_{i=0}^{n} W_i \right| + \left| X_n \cdot \prod_{i=0}^{n-1} W_i \right| \\ &< e^{c_3 \cdot n} \cdot (|X_{n+1}| + |X_n|) \\ &\leq e^{c_3 \cdot n} \cdot [\zeta \cdot (n+1)^d + \zeta \cdot n^d] \\ &\leq \lambda \cdot n^d \cdot e^{c_3 \cdot n} \end{aligned} $$

where the first inequality is induced by the triangle inequality, and the second inequality is derived from the bounded stepwise weight condition such that each $W_i \in [0, e^{c_3}]$ and the fact $W_0 = 1$. By applying the variant of Optional Stopping Theorem (Theorem 4.8), we obtain that $\mathbb{E}[-Y_T] \leq \mathbb{E}[-Y_0]$, so $\mathbb{E}[Y_T] \geq \mathbb{E}[Y_0]$. By definition and Condition (C2'),

$$ \begin{aligned} -Y_T &= -h(\ell_T, \mathbf{v}_T) \cdot \prod_{i=0}^{T-1} W_i \\ &= -h(\ell_{\text{out}}, \mathbf{v}_T) \cdot \prod_{i=0}^{T-1} W_i \\ &= -\prod_{i=0}^{T-1} W_i \end{aligned} $$

Finally, we have that $\mathbb{E}_{\Sigma}[W_{\infty}] = \mathbb{E}[\prod_{i=0}^{T-1} W_i] \geq \mathbb{E}[Y_0] = h(\ell_{\text{init}}, \mathbf{v}_{\text{init}})$. $\square$

## B.5 Correctness of Truncation

**Theorem 4.11.** Let $\Pi$ be a WPTS in the form of $(\uparrow)$, $\mathcal{B}$ a truncation function and $\mathcal{M}$ an approximation function. Suppose that the following condition (*) holds:

(*) for each fork $F^{M, \#} = \langle \# , p, \text{upd}, \mathcal{M} \rangle$ in the truncated WPTS $\Pi_{\mathcal{B}, \mathcal{M}}$ that is derived from some fork $F = \langle \ell', p, \text{upd}, \text{wt} \rangle$ with the source location $\ell$ in the original WPTS (see sentence (b) in Definition 4.10), we have that $[\Pi]_{\text{v}}(\mathbb{R}^{|\mathcal{V}_p|}) \leq \mathcal{M}(\text{v})$ for all $\text{v}$ such that the state $(\ell, \text{v})$ is reachable and $\text{v} \nleq \Phi_{\mathcal{B}}$.

Then $[\Pi]_{\text{v}_{\text{init}}}(\mathbb{R}^{|\mathcal{V}_p|}) \leq [\Pi_{\mathcal{B}, \mathcal{M}}]_{\text{v}_{\text{init}}}(\mathbb{R}^{|\mathcal{V}_p|})$ for all initial program valuations $\mathbf{v}_{\text{init}}$. Analogously, if it holds the condition $(\star)$ which is almost the same as (*) except for that "$[\Pi]_{\text{v}}(\mathbb{R}^{|\mathcal{V}_p|}) \leq \mathcal{M}(\text{v})$" is replaced with "$[\Pi]_{\text{v}}(\mathbb{R}^{|\mathcal{V}_p|}) \geq \mathcal{M}(\text{v})$", then we have $[\Pi]_{\text{v}_{\text{init}}}(\mathbb{R}^{|\mathcal{V}_p|}) \geq [\Pi_{\mathcal{B}, \mathcal{M}}]_{\text{v}_{\text{init}}}(\mathbb{R}^{|\mathcal{V}_p|})$ for all initial program valuations $\mathbf{v}_{\text{init}}$.
PROOF. We first prove that when every score function $\mathcal{M}$ in a $F^{\mathcal{M}, \#}$ derived from a transition with source location $\ell$ is equal to the function $ewt(\ell, -)$, we have that $\Pi_{\mathcal{B}, \mathcal{M}}$ is equal to $\Pi$. By Theorem 4.5, the expected weight functions $ew_{\Pi}$, $ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$ are the least fixed point of the higher-order operator $ewt$ defined in Definition 4.3. We prove that both $ew_{\Pi} \le ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$ and $ew_{\Pi} \ge ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$ holds. Note that since we choose the scoring function to be the exact expected weight function of $\Pi$, it holds that $ew_{\Pi}(-, \mathbf{v}) = ew_{\Pi_{\mathcal{B}, \mathcal{M}}}(-, \mathbf{v})$ for all program valuations outside $B$. Thus, the nontrivial part is to consider program valuations inside the truncated range.

* $ew_{\Pi} \le ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$: To show that $ew_{\Pi} \le ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$, it suffices to observe that $ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$ satisfies $ewt_{\Pi}(ew_{\Pi_{\mathcal{B}, \mathcal{M}}}) = ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$. Since $ew_{\Pi}$ is the least fixed point of the higher order equation, we directly obtain that $ew_{\Pi} \le ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$.

* $ew_{\Pi} \ge ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$: To show that $ew_{\Pi} \ge ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$, it suffices to observe that $ew_{\Pi}$ (extended with the $\#)$ location whose score function is 1) satisfies the higher-order equation of $\Pi_{\mathcal{B}, \mathcal{M}}$. Thus, we directly have that $ew_{\Pi} \ge ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$.

Then we prove the theorem. We only prove the upper-bound case, since the lower-bound case can be proved similarly. The proof follows from Theorem 4.5. Denote $\perp$ as the bottom element of the complete lattice ($\mathcal{K}_M, \leq$). Then by Theorem 4.5, we have that $\lim_{n \to \infty} ew_{\Pi}^n(\perp) = ew_{\Pi}$ and $\lim_{n \to \infty} ew_{\Pi_{\mathcal{B}, \mathcal{M}}}^n(\perp) = ew_{\Pi_{\mathcal{B}, \mathcal{M}}}$. Since $ew_{\Pi}(\mathbf{v}) \le \mathcal{M}(\mathbf{v})$ for all $\mathbf{v} \in \text{exit}(\Pi)$ where $\text{exit}(\Pi)$ is the set of all program valuations $\mathbf{v}$ such that the state $(\ell_{\text{out}}, \mathbf{v})$ is reachable, one can perform a straightforward induction on $n$ that $ew_{\Pi_{\mathcal{B}}}^n, ew_{\Pi} \le ew_{\Pi_{\mathcal{B}}}^{n, f}$ for all $n$. $\square$

# C SUPPLEMENTARY MATERIAL FOR SECTION 5

## C.1 Error Analysis for Polynomial Approximation

THEOREM C.1. Let $\Pi$ be a non-score-recursive WPTS with score functions $g_1, \dots, g_k$ on the transitions to the termination location $\ell_{\text{out}}$. Suppose we have a non-negative real number $\epsilon$ and polynomials $g'_1, \dots, g'_k$ such that for all $x \in \text{exit}(\Pi)$ and $1 \le j \le k$, $|g'_j(x) - g_j(x)| \le \epsilon$. Then we have that $|\llbracket\Pi\rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_{\mathbf{P}}|}) - \llbracket\Pi'\rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_{\mathbf{P}}|})| \le \epsilon$ for all initial program valuation $\mathbf{v}$, where $\Pi'$ is obtained from $\Pi$ by replacing each $g_j$ ($1 \le j \le k$) with $g'_j$.

PROOF. By Theorem 4.5, we have that $\lim_{n \to \infty} ewt_{\Pi}^n(\perp) = ew_{\Pi}$ and $\lim_{n \to \infty} ewt_{\Pi'}^n = ew_{\Pi'}$. Since $|g'_j - g_j| \le \epsilon$ for every $j$, one can perform a straightforward induction on $n$ to prove that for all $n \ge 0$, it holds that $|ewt_{\Pi}^n(\perp) - ewt_{\Pi'}^n(\perp)| \le \epsilon$. By the definition of $\llbracket\Pi\rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_{\mathbf{P}}|})$, we have that $|\llbracket\Pi\rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_{\mathbf{P}}|}) - \llbracket\Pi'\rrbracket_{\mathbf{v}}(\mathbb{R}^{|V_{\mathbf{P}}|})| \le \epsilon$. $\square$

In **Step A4** of our algorithm, constraints are established in the form $\forall \mathbf{v} \in P.(g(\mathbf{v}) \ge 0)$ where $P$ is a polyhedron over program variables $V_{\mathbf{P}}$ and are grouped conjunctively. Thus, the key point is how to tackle the each such constraint. In our algorithm, we follow the exact treatment through the application of Putinar's and Handelman's Positivstellensatz in [24, 36]. Below we describe the detailed application.

## C.2 Application of Putinar's Positivstellensatz

We recall Putinar's Positivstellensatz below.

THEOREM C.2 (PUTINAR'S POSITIVSTELLENSATZ [36]). Let $V$ be a finite set of real-valued variables and $g, g_1, \dots, g_m \in \mathbb{R}[V]$ be polynomials over $V$ with real coefficients. Consider the set $S := \{\mathbf{x} \in \mathbb{R}^V \mid g_i(\mathbf{x}) \ge 0 \text{ for all } 1 \le i \le m\}$ which is the set of all real vectors at which every $g_i$ is non-negative.
If (i) there exists some $g_k$ such that the set $\{\mathbf{x} \in \mathbb{R}^V \mid g_k(\mathbf{x}) \ge 0\}$ is compact and (ii) $g(\mathbf{x}) > 0$ for all $\mathbf{x} \in S$, then we have that

$$g = f_0 + \sum_{i=1}^{m} f_i \cdot g_i \qquad (8)$$

for some polynomials $f_0, f_1, \dots, f_m \in \mathbb{R}[V]$ such that each polynomial $f_i$ is the a sum of squares (of polynomials in $\mathbb{R}[V]$), i.e. $f_i = \sum_{j=0}^{k} q_{i,j}^2$ for polynomials $q_{i,j}$'s in $\mathbb{R}[V]$.

In this work, we utilize the sound form in (8) for witnessing a polynomial $g$ to be non-negative over a polyhedron $P$ for each constraint $\forall \mathbf{v} \in P.(g(\mathbf{v}) \ge 0)$ from **Step A5** of our algorithm. Let $\forall \mathbf{v} \in P.(g(\mathbf{v}) \ge 0)$ be such a constraint for which the polyhedron $P$ is defined by the linear inequalities $g_1 \ge 0, \dots, g_m \ge 0$. Let $V_p = \{v_1, v_2, \dots, v_t\}$ be the set of program variables and define $M_d(V_p) = \{m_1, m_2, \dots, m_r\}$ as the set of all monomials of degree at most $d$ over $V_p$, i.e. $M_d(V_p) := \{\prod_{i=1}^t v_i^{\alpha_i} \mid \forall i \alpha_i \in \mathbb{N} \wedge \sum_{i=1}^t \alpha_i \le d\}$. The application of Putinar's to $\forall \mathbf{v} \in P.(g(\mathbf{v}) \ge 0)$ has the following steps.

* First, represent each $f_i$ in Eq. (8) as the positive semidefinite form $f_i = \mathbf{v}^T Q_i \mathbf{v}$ subject to the positive semidefinite constraint where each $Q_i$ is a real matrix whose every entry is an unknown parameter.

* Second, compute an equation in the form (8) whose coefficients are affine expressions in the unknown coefficients from our templates and the unknown entries in the matrices $Q_i$'s.

* Third, establish the affine constraints between the unknown coefficients in the templates and the unknown entries in the matrices $Q_i$'s by matching the coefficients at the LHS and the RHS of the equation obtained from the previous step.

The overall application processes all such constraints from **Step A5** of our algorithm by (i) collecting all the affine and the semidefinite constraints from the first and the third steps above and (ii) solve them by semidefinite programming.

### C.3 Application of Handelman's Positivstellensatz

To present Handelman's Positivstellensatz, we need the notion of monoid as follows. Below we consider an arbitrary finite collection $\Gamma = \{g_1, \dots, g_k\} (k \ge 1)$ of linear functions (i.e., degree-1 polynomials) in the program variables.

**Definition C.3 (Monoid).** The monoid of $\Gamma$ is defined by:

$$\text{Monoid}(\Gamma) := \left\{ \prod_{i=1}^{k} h_i \mid k \in \mathbb{N}_0 \text{ and } h_1, \dots, h_k \in \Gamma \right\}.$$

Then in our context, Handelman's Positivstellensatz can be formulated as follows.

**THEOREM C.4 (HANDELMAN'S POSITIVSTELLENSATZ [24]).** Let $g$ be a polynomial in the program variables such that $g(\mathbf{v}) > 0$ for all program valuations $\mathbf{v} \in P := \{\mathbf{v}' \in \mathbb{R}^{|V_p|} \mid g_1(\mathbf{v}') \ge 0, \dots, g_k(\mathbf{v}') \ge 0\}$. If $P$ is compact, then we have

$$g = \sum_{i=1}^{d} a_i \cdot u_i \qquad (9)$$

for some $d \in \mathbb{N}$, real numbers $a_1, \dots, a_d \ge 0$ and $u_1, \dots, u_d \in \text{Monoid}(\Gamma)$.

To apply Handelman's Positivstellensatz, we consider a natural number which serves as a bound on the number of multiplicands allowed to form an element in $\text{Monoid}(\Gamma)$. Then Eq. (9) results in a system of linear equalities that involves $a_1, \dots, a_d$ and the coefficients of $g$. The application of Handelman's Positivstellensatz to each $\forall \mathbf{v} \in P.(g(\mathbf{v}) \ge 0)$ is simpler than that of Putinar's Positivstellensatz, and is as follows.
* First, compute an equation in the form (9) whose coefficients are affine expressions in the unknown coefficients from our templates and the fresh variables $a_1, \dots, a_d$ from Eq. (9).

* Second, establish the affine constraints between the unknown coefficients in the templates and the fresh variables $a_1, \dots, a_d$ from Eq. (9) by matching the coefficients at the LHS and the RHS of the equation obtained from the previous step.

The overall application processes all such constraints from **Step A5** of our algorithm by (i) collecting all the affine constraints from the second steps above and (ii) solve them by linear programming.

# D SUPPLEMENTARY MATERIALS FOR SECTION 6

## D.1 Possible Approaches for Computing $M_{up}$ and $M_{low}$ of Score-recursive WPTS's

Fix a score-recursive WPTS $\Pi$, and assume it has (1) the concentration property, i.e., $\mathbb{P}(T > n) \le c_1 \cdot e^{-c_2 \cdot n}$ for $c_1, c_2 > 0$, (2) the bounded-update property, and (3) the stepwise weight is bounded by $e^{c_3}$ for $0 < c_3 < c_2$. Then given the bounded ranges $B$ and $B'$ as computed in Section 5, we derive the upper bound $M_{up}$ and the lower bound $M_{low}$ for the expected weight from $B'\setminus B$ as follows.

For any $v \in B'\setminus B$,

$$ 
\begin{aligned}
\llbracket \Pi \rrbracket(v) &= \mathbb{E}_v[w_T] \\
&= \sum_{n=1}^{\infty} \mathbb{P}(T = n) \cdot w_n \\
&\le \sum_{n=0}^{\infty} \mathbb{P}(T > n) \cdot w_n \\
&\le 1 + \sum_{n=1}^{\infty} \mathbb{P}(T > n) \cdot w_n \\
&= 1 + \sum_{n=1}^{n^*-1} \mathbb{P}(T > n) \cdot w_n + \sum_{n=n^*}^{\infty} \mathbb{P}(T > n) \cdot w_n \\
&\le 1 + M + \sum_{n=n^*}^{\infty} \mathbb{P}(T > n) \cdot w_n \\
&= M' + \sum_{n=n^*}^{\infty} \mathbb{P}(T > n) \cdot w_n \\
&\le M' + \sum_{n=n^*}^{\infty} c_1 \cdot e^{-c_2 \cdot n} \cdot w_n \\
&\le M' + \sum_{n=n^*}^{\infty} c_1 \cdot e^{-c_2 \cdot n} \cdot e^{c_3 \cdot n} \\
&= M' + c_1 \cdot \sum_{n=n^*}^{\infty} (e^{c_3-c_2})^n \\
&= M' + c_1 \cdot \frac{a}{1-q} \\
&= M_{up}
\end{aligned}
$$
where

$$ \sum_{n=1}^{n^*-1} \mathbb{P}(T > n) \cdot w_n \le \sum_{n=1}^{n^*-1} w_n \le \sum_{n=1}^{n^*-1} (e^{c_3})^n = \frac{a' \cdot (1 - (q')^{(n^*-1)})}{1 - q'} =: M $$

and $a' = e^{c_3}$, $q' = e^{c_3}$. The first inequality is obtained from the fact that

$$ \mathbb{P}(T > n) = \mathbb{P}(T \ge n + 1) = \mathbb{P}(T = n + 1) + \mathbb{P}(T = n + 2) + \dots, $$

thus,

$$ \mathbb{P}(T = n + 1) \le \mathbb{P}(T > n). $$

The second inequality is derived by the fact that $\mathbb{P}(T > 0) \le 1$ and $w_0 = w_{init} = 1$. The third inequality is obtained by the definition of M above. The fourth inequality is obtained by the concentration property, while the fifth inequality is derived by the bounded stepwise weight condition.

For $M_{low}$, we trivially set $M_{low} = 0$. We can refine it heuristically, e.g., according to the monotonicity of the scoring function.

## D.2 Overapproximation via Polynomial Interpolations

Given a non-polynomial function $f(x)$ over the interval $I = [a, b]$, we aim to approximate $f(x)$ by polynomials $p(x)$'s. The correctness of approximation is based on a classical theorem called Weierstrass' Theorem [26].

**THEOREM D.1 (WEIERSTRASS' THEOREM).** Let $f(x)$ be a continuous function on the (closed) interval $[a, b]$. Then there is a sequence of polynomials $p_n(x)$ (of degree $n$) such that

$$ \lim_{n \to \infty} \|f - p_n\|_{\infty} = 0. $$

We also need the following theorem to measure the derived polynomials. The property of Lipschitz continuity supports the following theorem easily.

**THEOREM D.2.** Suppose $r(x)$ is a continuous and differentiable function on a compact convex set $\Psi \subseteq \mathbb{R}$. Assume that a collection of points $\{x_1, x_2, \dots, x_k\}$ are sampled uniformly from $\Psi$ and $s \in \mathbb{R}_{>0}$ is the sampling spacing. Let $r_0 = \max\{|r(x_1)|, |r(x_2)|, \dots, |r(x_k)|\}$, and $\beta = \sup_{x \in \Psi} ||\nabla r(x)||$, then

$$ |r(x)| \le \beta \cdot s + r_0, \quad \forall x \in \Psi. \qquad (10) $$

Then our scheme is as follows.

* Split the interval $I = [a, b]$ uniformly into $m$ partitions, i.e., $I_1 = [a_1, b_1]$, $I_2 = [a_2, b_2]$, ..., $I_m = [a_m, b_m]$.

* For each partition $I_i = [a_i, b_i]$, define a $n$-degree polynomial $p_n^i(x) := \sum_{j=0}^{n} c_{ij} \cdot x^j$.

(1) Pick a non-negative integer $k > n$ and sample $k$ points uniformly from $f$ over $I_i$. That is,

$$ D = \{(x_1, f(x_1)), (x_2, f(x_2)), \dots, (x_k, f(x_k))\} $$

where $x_l \in I_i$ for all $1 \le l \le k$.

(2) Let $p_n^i(x_l) = f(x_l)$ for all $1 \le l \le k$, then we have a linear system $V \cdot c = f$ where

$$ V = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_k & x_k^2 & \cdots & x_k^n \end{bmatrix}, $$

$$ c = [c_{i0}, c_{i1}, \dots, c_{in}]^T \text{ and } f = [f(x_1), f(x_2), \dots, f(x_k)]^T. $$
(3) By solving the above overdetermined system, we obtain $p_n^i(x)$ as the approximation of $f(x)$ over the interval $I_i = [a_i, b_i]$.

(4) Having $p_n^i(x)$, evaluate an error bound $\gamma_i$ such that

$$ \forall x \in I_i, |f(x) - p_n^i(x)| \le \gamma_i. \qquad (11) $$

Let $r(x) = f(x) - p_n^i(x)$ and $\Psi = I_i$, then we obtain $r_0 = \max\{|r(x_1)|, |r(x_2)|, \dots, |r(x_k)|\}$ by Theorem D.2. To derive the Lipschitz constant $\beta$ of $r(x)$ over the interval $I_i$, we pick a non-negative integer $q = 10k$, and sample $q$ points uniformly from $f$, i.e., we have another collection of points $\{x'_1, x'_2, \dots, x'_q\}$. Let $\beta = \max\{|\nabla r(x'_1)|, \dots, |\nabla r(x'_q)|\}$, then

$$ \gamma_i := \beta \cdot s + r_0 $$

where $s$ is the corresponding sampling spacing of the $q$ points.

* Now we have a set $D_p$ of tuples of intervals, polynomials and error bounds, i.e.,

$$ D_p = \{(I_1, p_n^1(x), \gamma_1), \dots, (I_m, p_n^m(x), \gamma_m)\} \qquad (12) $$

The approximation error bounds $\gamma_i$'s are taken into account when we synthesize the polynomial template $h$. Given a non-polynomial function $f(x)$ such that score($f(x)$) occurs in the program, we obtain a set $D_p$ in the form of (12). For each interval $I_i$, we introduce a new variable $r_i$ and approximate $f(x)$ over $I_i$ as $p_n^i(x) + r_i$ with $r_i \in [-\gamma_i, \gamma_i]$. That is, for $1 \le i \le m$, we have

$$ \forall x \in I_i, f(x) \approx p_n^i(x) + r_i \text{ with } r_i \in [-\gamma_i, \gamma_i]. \qquad (13) $$

For a state $(\ell, \mathbf{v})$ such that $\ell$ is the location before the command score($f(x)$), there is the unique transition $\tau = \langle \ell, \text{true}, F \rangle$ such that $F = \langle \ell', 1, \mathbf{1}, f \rangle$ and $\ell'$ is the location that follows the command score($f(x)$). Then for all valuations $\mathbf{v} \in I(\ell) \wedge \Phi_B$ and $1 \le i \le m$, it should hold that

* for all $\mathbf{v}[x] \in I_i$ and $r_i \in [-\gamma_i, \gamma_i]$, we have that $\text{ewt}(h)(\ell, \mathbf{v}) \le h(\ell, \mathbf{v})$ (for upper bounds) and $\text{ewt}(h)(\ell, \mathbf{v}) \ge h(\ell, \mathbf{v})$ (for lower bounds) where

$$ \text{ewt}(h)(\ell, \mathbf{v}) = (p_n^i(x) + r_i) \cdot h(\ell', \mathbf{v}). $$

## D.3 Other Experimental Results
<figure><img src="image_16.png" /><figcaption>Fig. 11. Part 1: NPD Bounds of Novel Examples</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>(a) Hare Tortoise Race v1</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Random Walk v2</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c) Random Walk v3</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) Random Walk v4</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(e) Pedestrian Multi-branches v4</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(f) Phylogenetic Birth Model</figcaption></figure>

<caption>Fig. 12. Part 2: NPD Bounds of Novel Examples</caption>
<figure><img src="image_23.png" /><figcaption>Fig. 13. Part 1: NPD Bounds of Our Approach and GuBPI</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 14. Part 2: NPD Bounds of Our Approach and GuBPI</figcaption></figure>