1

# This Is the Moment for Probabilistic Loops

Marcel Moosbrugger  0000-0002-2006-3741 TU WienViennaAustria
<EMAIL_ADDRESS>,  Miroslav Stankovič  0000-0001-5978-7475 TU
WienViennaAustria<EMAIL_ADDRESS>,  Ezio Bartocci
0000-0002-8004-6601 TU WienViennaAustria<EMAIL_ADDRESS>and  Laura
Kovács  0000-0002-8299-2714 TU WienViennaAustria<EMAIL_ADDRESS>
(2022)

###### Abstract.

We present a novel static analysis technique to derive higher moments for
program variables for a large class of probabilistic loops with potentially
uncountable state spaces. Our approach is fully automatic, meaning it does not
rely on externally provided invariants or templates. We employ algebraic
techniques based on linear recurrences and introduce program transformations
to simplify probabilistic programs while preserving their statistical
properties. We develop power reduction techniques to further simplify the
polynomial arithmetic of probabilistic programs and define the theory of
moment-computable probabilistic loops for which higher moments can precisely
be computed. Our work has applications towards recovering probability
distributions of random variables and computing tail probabilities. The
empirical evaluation of our results demonstrates the applicability of our work
on many challenging examples.

Probabilistic Programs, Higher Moments, Linear Recurrences, Distribution
Recovery

††journal: PACMPL††journalvolume: 6††journalnumber: OOPSLA2††journalyear:
2022††submissionid: oopslab22main-p505-p††publicationmonth: 10††doi:
10.1145/3563341††copyright: rightsretained††ccs: Mathematics of computing
Markov processes††ccs: Computing methodologies Symbolic and algebraic
algorithms††ccs: Theory of computation Random walks and Markov chains

##  1\. Introduction

Probabilistic programming languages enrich classical imperative or functional
languages with native primitives to draw samples from random distributions,
such as Bernoulli, Uniform, and Normal distributions. The resulting
probabilistic programs (PPs) (Kozen, 1985; Barthe et al., 2020) embed
uncertain quantities, represented by random variables, within standard program
control flows. As such, PPs offer a unifying framework to naturally encode
probabilistic machine learning models (Ghahramani, 2015), for example Bayesian
networks (Kaminski et al., 2016), into programs. Moreover, PPs enable
programmers to handle uncertainty resulting from sensor measurements and
environmental perturbations in cyber-physical systems (Selyunin et al., 2015;
Chou et al., 2020). Other notable examples of PPs include the implementation
of cryptographic (Barthe et al., 2012a) and privacy (Barthe et al., 2012b)
protocols, as well as randomized algorithms (Motwani and Raghavan, 1995) such
as Herman’s self-stabilization protocol (Herman, 1990) for recovering from
faults in a process token ring — see our example in Figure 2.

#### Analysis of PPs.

The random nature of PPs makes their functional analysis very challenging as
one needs to reason about probability distributions of random variables
instead of computing with single variable values (Barthe et al., 2020). A
standard approach towards handling probability distributions associated with
random variables is to estimate such distributions by sampling PPs using Monte
Carlo simulation techniques (Hastings, 1970). While such approaches work well
for statistical model checking (Younes and Simmons, 2006), they are not
suitable for the analysis of PPs with potentially infinite program loops as
simulating infinite-state behavior is not always viable. Moreover, even for
PPs with finitely many states, simulation-based analysis is inherently
approximative.

With the aim of precisely, and not just approximately, handling random
variables, probabilistic model checking (Kwiatkowska et al., 2011; Dehnert et
al., 2017) became a prominent approach in the analysis of PPs with finite
state spaces. For analyzing unbounded PPs, these techniques would however
require non-trivial user guidance, in terms of assertion templates and/or
invariants.

_In this paper, we address the challenge of precisely analyzing, and even
recovering, probability distributions induced by PPs with both countably and
uncountably infinite state spaces._ We do so by extending both expressivity
and automation of the state-of-the-art in PP analysis: We (i) focus on PPs
with probabilistic infinite loops (see Figure 1) and (ii) fully automate the
analysis of such loops by computing exact higher-order statistical moments of
program variables $x$ parameterized by a loop counter $n$.

Functional representations $f(n)$ for a program variable $x$, with $f(n)$
characterizing the $k$th moment ${\mathbb{E}}(x^{k}_{n})$ of $x$ at iteration
$n$, can be interpreted as a quantitative invariant
${\mathbb{E}}(x^{k}_{n})-f(n)=0$, as the equation is true for all loop
iterations $n\in{\mathbb{N}}$. Inferring quantitative invariants is arguably
not novel. On the contrary, it is one of the most challenging aspects of PP
analysis, dating back to the seminal works of (McIver and Morgan, 2005; Katoen
et al., 2010) introducing the weakest pre-expectations calculus. Template-
based approaches to discover invariants or (super-)martingales emerged (Barthe
et al., 2016; Kura et al., 2019) by translating the invariant generation
problem into a constraint solving one. The derived quantitative invariants are
generally provided in terms of expected values (Chakarov and Sankaranarayanan,
2014; Katoen et al., 2010; McIver and Morgan, 2005). Nevertheless, the
expected value alone — also referred to as the _first moment_ — provides only
partial information about the underlying probability distribution. This
motivates the critical importance of higher moments for PP analysis (Kura et
al., 2019; Bartocci et al., 2020b; Wang et al., 2021; Stankovic et al., 2022).

#### Higher Moments for PP Analysis.

Using concentration-of-measure inequalities (Boucheron et al., 2013), we can
utilize higher moments ${\mathbb{E}}(X^{k})$ to obtain upper and lower bounds
on tail probabilities ${\mathbb{P}}(X>t)$, measuring the probability that a
given random variable $X$, corresponding for example to our program variables
$x$ from Figure 1, surpasses some value $t$. In this paper, we also show that
when a program variable $x$ admits only $k<\infty$ many values, we can fully
recover its probability mass function as a closed-form expression in the loop
counter $n$ using the first $k{-}1$ raw moments (see Section 6). Furthermore,
raw moments can be used to compute central moments
${\mathbb{E}}((X-{\mathbb{E}}(X))^{k})$ and thus provide insights on other
important characteristics of the distribution such as the _variance_ ,
_skewness_ and _kurtosis_ (Durrett, 2019). However, computing exact higher
statistical moments for PPs is computationally expensive (Kaminski et al.,
2019), a challenge which we also address in this paper, as illustrated in
Figures 1–2 and described next.

⬇

toggle,sum,x,y,z = 0,$s_{0}$,1,1,1

while $\star$:

toggle = 1-toggle

if toggle == 0:

x = x+1 {1/2} x+2

y = y+z+x**2 {1/3} y-z-x

z = z+y {1/4} z-y

end

l,g = Laplace(x+y, 1),Normal(0,1)

if g < 1/2: sum = sum+x end

end

${\mathbb{E}}(\text{toggle}_{n})=\frac{1}{2}-\frac{(-1)^{n}}{2}$
${\mathbb{E}}(x_{n})=\frac{5}{8}+\frac{3n}{4}+\frac{3(-1)^{n}}{8}$
${\mathbb{E}}(x_{n}^{2})=\frac{15}{32}+\frac{17n}{16}+\frac{9n(-1)^{n}}{16}+\frac{17(-1)^{n}}{32}+\frac{9n^{2}}{16}$
${\mathbb{E}}(l_{n})=\frac{-17}{8}-\frac{15n}{4}+\frac{67\cdot 2^{-n}\cdot
6^{\frac{n}{2}}}{10}+\frac{67\cdot 2^{-n}6^{\frac{1+n}{2}}}{30}-$ $\text{\ \ \
}\frac{37\cdot 3^{-n}6^{\frac{n}{2}}}{10}-\frac{37\cdot
3^{-n}6^{\frac{1+n}{2}}}{20}+\frac{67\cdot 6^{\frac{n}{2}}(-1)^{n}}{10\cdot
2^{n}}-$ $\text{\ \ \ }\frac{15(-1)^{n}}{8}-\frac{67\cdot
6^{\frac{1+n}{2}}(-1)^{n}}{30\cdot 2^{n}}+\frac{37\cdot
6^{\frac{1+n}{2}}(-1)^{n}}{20\cdot 3^{n}}-$ $\text{\ \ \ }\frac{37\cdot
6^{\frac{n}{2}}(-1)^{n}}{10\cdot 3^{n}}$

Figure 1. An example of a multi-path PP loop, with Laplace and Normal
distributions parametrized by program variables. Our work fully automates the
analysis of such and similar PP loops by computing higher moments. Several
moments for program variables in the loop counter $n$ are listed on the right.
Each moment was automatically generated.

#### Computing Higher Moments.

The theory we establish in this paper describes how to compute higher moments
of program variables for a large class of probabilistic loops and how to
utilize these moments to gain more insights into the analyzed programs. We
call this theory the _theory of moment-computable probabilistic loops_
(Section 5). Our approach is fully automatic, meaning it does _not_ rely on
externally provided invariants or templates. Unlike constraint solving over
templates (Barthe et al., 2016; Kura et al., 2019), we employ algebraic
techniques based on systems of _linear recurrences_ with constant coefficients
describing so-called _C-finite sequences_ (Kauers and Paule, 2011). Different
equivalence preserving program transformations (Section 3) and power reduction
of finite valued variables (Section 4) allow us to simplify PPs and represent
their higher moments as linear recurrence systems in the loop counter. Figure
1 shows a PP with many unique features supported by our work towards PP
analysis: it has an uncountable state-space, contains if-statements, symbolic
constants, draws from continuous probability distributions with state-
dependent parameters, and employs polynomial arithmetic as well as circular
variable dependencies. We are not aware of other works automating the
reasoning about such and similar probabilistic loops, in particular for
computing precise higher moments of variables. Figure 1 lists some of the
variables’ moments computed automatically by our work. Further, these moments
can be used to compute tail probability bounds or central moments, such as the
variance, to characterize the distribution of the program variables as the
loop progresses.

Thanks to our power reduction techniques (Section 4), our approach supports
arbitrary polynomial dependencies among finite valued variables. Moreover, our
work can fully recover the value distributions of finite valued program
variables, from finitely many higher moments, as illustrated in Section 6 for
Herman’s self-stabilization algorithm from Figure 2.

#### Theory and Practice in Computing Higher Moments.

In theory, our approach can compute any higher moment for any variable and PP
of our program model, under assumptions stated in Sections 3 and 5. We also
establish the necessity of these assumtions in Section 5.3. In a nutshell, the
completeness theorem (Theorem 3) holds for probabilistic loops for which non-
finite program variables are not polynomially self-dependent and all branching
conditions are over finite valued variables. We strengthen the theory of
(Bartocci et al., 2019) to support if-statements, circular variable
dependencies, state-dependent distribution parameters, simultaneous
assignments, and multiple assignments, and establish the necessity of our
assumptions. Moreover, unlike (Wang et al., 2021), our approach does not rely
on templates and provides exact closed-form representations of higher moments
parameterized by the loop counter.

In practice, our approach is implemented in the Polar tool and compared
against exact as well as approximate methods (Moosbrugger et al., 2022). Our
experiments (Section 7) show that Polar outperforms the state-of-the-art of
moment computation for probabilistic loops in terms of supported programs and
efficiency. Furthermore, Polar is able to compute exact higher moments
magnitudes faster than sampling can establish reasonable confidence intervals.

#### Contributions.

Our main contributions are listed below:

  * •

An automated approach for computing higher moments of program variables for a
large class of probabilistic loops with potentially uncountable state spaces
(Sections 3-5).

  * •

We develop power reduction techniques to reduce the degrees of finite valued
program variables in polynomials (Section 4).

  * •

We prove completeness of our work for computing higher moments (Section 5).

  * •

We fully recover the distributions of finite valued program variables and
approximate distributions for unbounded/continuous program variables from
finitely many moments (Section 6).

  * •

We provide an implementation and empirical evaluation of our work,
outperforming the state-of-the-art in PP analysis in terms of automation and
expressivity (Section 7).

⬇

x1, x2, x3 = 1, 1, 1

t1, t2, t3 = 1, 1, 1

p = 1/2; tokens = t1 + t2 + t3

while $\star$:

x1o, x2o, x3o = x1, x2, x3

if x1o == x3o: x1=Bernoulli(p) else: x1=x3o end

if x2o == x1o: x2=Bernoulli(p) else: x2=x1o end

if x3o == x2o: x3=Bernoulli(p) else: x3=x2o end

if x1 == x3: t1 = 1 else: t1 = 0

if x2 == x1: t2 = 1 else: t2 = 0

if x3 == x2: t3 = 1 else: t3 = 0

tokens = t1 + t2 + t3

end

${\mathbb{E}}(tokens_{n})=1+2\cdot 4^{-n}$ —
${\mathbb{E}}(tokens^{2}_{n})=1+8\cdot 4^{-n}$ —
${\mathbb{E}}(tokens^{3}_{n})=1+26\cdot 4^{-n}$

Figure 2. Herman’s self stabilization algorithm with three nodes encoded as a
probabilistic loop together with three moments of $tokens$.

##  2\. Preliminaries

We use the symbol ${\mathbb{P}}$ for probability measures and ${\mathbb{E}}$
for the expectation operator. The support of a random variable $X$ is denoted
by $\text{supp}(X)$.

###  2.1. Probability Theory

Operationally, a probabilistic program is a Markov chain with potentially
uncountably many states. Let us recall some notions about Markov chains. For
more details on Markov chains and probability theory in general we refer the
reader to (Durrett, 2019).

For a fixed set $S$, a $\sigma$-algebra is a non-empty set of subsets of $S$
closed under complementation and countable unions.

###### Definition 0 (Sequence Space).

Let $(S,\mathcal{S})$ be a measurable space, that is, $S$ is a set with a
$\sigma$-algebra $\mathcal{S}$. Its _sequence space_ is the measurable space
$(S^{\omega},\mathcal{S}^{\omega})$ where
$S^{\omega}:=\\{(s_{1},s_{2},\dots):s_{i}\in S\\}$ and $\mathcal{S}^{\omega}$
is the $\sigma$-algebra generated by the _cylinder sets_
$Cyl[B_{1},\dots,B_{n}]:=\\{\theta:\theta_{i}\in B_{i},1\leq i\leq n\\}$ for
all prefixes $B_{1},\dots,B_{n}\in\mathcal{S}$ and all $n\in{\mathbb{N}}$.

A _Markov kernel_ is, on a high level, a generalization of transition
probabilities between states to uncountable state spaces and is required for
the definition of a _Markov chain_.

###### Definition 0 (Markov Chain).

Let $(S,\mathcal{S},{\mathbb{P}})$ be a probability space and
$p:S\times\mathcal{S}\to[0,1]$ a Markov kernel. A stochastic process $X_{n}$
is a _Markov chain_ with Markov kernel $p$ if

(1) ${\mathbb{P}}(X_{n+1}\in B\mid
X_{0}=x_{0},X_{1}=x_{1},\dots,X_{n}=x_{n})=p(X_{n},B).$

Given a measurable space $(S,\mathcal{S})$ , an _initial distribution_ $\mu$,
a stochastic process $X_{n}$ and a Markov kernel $p$, _Kolmogorov’s Extension
Theorem_ says that there is a unique measure ${\mathbb{P}}$ such that $X_{n}$
is a Markov chain in $(S^{\omega},\mathcal{S^{\omega}},{\mathbb{P}})$.

For a random variable $X$, central moments
${\mathbb{E}}((X-{\mathbb{E}}(X))^{k})$ can be computed from raw moments
${\mathbb{E}}(X^{k})$ and vice versa through the transformation of center:

(2)
${\mathbb{E}}\left((X-b)^{k}\right)={\mathbb{E}}\left(((X-a)+(a-b))^{k}\right)=\sum_{i=0}^{k}\binom{k}{i}{\mathbb{E}}\left((X-a)^{i}\right)(a-b)^{k-i}.$

###  2.2. Linear Recurrences

We briefly recall standard terminology on algebraic sequences and recurrences.
For further details, we refer the reader to (Kauers and Paule, 2011). A
sequence $(a_{n})_{n=0}^{\infty}$ is called _C-finite_ if it obeys a linear
recurrence with constant coefficients, that is, $(a_{n})_{n=0}^{\infty}$
satisfies an equation of the form

$a_{n+l}=c_{l-1}\cdot a_{n-l-1}+c_{l-2}\cdot a_{n-l-2}+\dots+c_{0}\cdot
a_{n},$

for some _order_ $l\in{\mathbb{N}}$, some constants $c_{i}\in{\mathbb{R}}$ and
all $n\in{\mathbb{N}}$.

###### Theorem 3 (Closed-form (Kauers and Paule, 2011)).

Every C-finite sequence $(a_{n})_{n=0}^{\infty}$ can be written as an
_exponential polynomial_ , that is $a_{n}=\sum_{i=1}^{m}n^{d_{i}}u_{i}^{n}$
for some natural numbers $d_{i}\in{\mathbb{N}}$ and complex numbers
$u_{i}\in{\mathbb{C}}$. We refer to $\sum_{i=1}^{m}n^{d_{i}}u_{i}^{n}$ as the
_closed-form_ or the _solution_ of the sequence $(a_{n})_{n=0}^{\infty}$ or
its recurrence.

An important fact is that closed-forms of linear recurrences with constant
coefficients of _any order_ always exist and are computable. This also holds
for all variables in _systems_ of linear recurrences with constant
coefficients.

##  3\. Probabilistic Program Model

In this section we introduce our programming model (Section 3.1) and describe
its semantics in terms of Markov chains (Section 3.2). Moreover, we introduce
transformations (3.3) normalizing a probabilistic program to simplify its
analysis.

###  3.1. Probabilistic Program Syntax

$\mathit{lop}\in\\{\mathit{and},\mathit{or}\\}$,
$\mathit{cop}\in\\{=,\neq,<,>,\geq,\leq\\}$,
$\mathit{Dist}\in\\{\mathit{Bernoulli},\mathit{Normal},\mathit{Uniform},\dots\\}$
¡sym¿ ::= ”a” — ”b” — $\dots$ ¡var¿ ::= ”x” — ”y” — $\dots$

¡const¿ ::= $r\in{\mathbb{R}}$ — ¡sym¿ — ¡const¿ ( ”+” — ”*” — ”/” ) ¡const¿

¡poly¿ ::= ¡const¿ — ¡var¿ — ¡poly¿ (”+” — ”-” — ”*”) ¡poly¿ — ¡poly¿”**n”

¡assign¿ ::= ¡var¿ ”=” ¡assign_right¿ — ¡var¿ ”,” ¡assign¿ ”,” ¡assign_right¿

¡categorical¿ ::= ¡poly¿ (”{”¡const¿”}” ¡poly¿)* [”{”¡const¿”}”]

¡assign_right¿ ::= ¡categorical¿ — Dist”(”¡poly¿∗”)” —
”Exponential(”¡const¿”/”¡poly¿”)”

¡bexpr¿ ::= ”true” ($\star$) — ”false” — ¡poly¿ ¡cop¿ ¡poly¿ — ”not” ¡bexpr¿ —
¡bexpr¿ ¡lop¿ ¡bexpr¿

¡ifstmt¿ ::= ”if” ¡bexpr¿”:” ¡statems¿ (”else if” ¡bexpr¿”:” ¡statems¿)∗
[”else:” ¡statems¿] ”end”

¡statem¿ ::= ¡assign¿ — ¡ifstmt¿ ¡statems¿ ::= ¡statem¿+

¡loop¿ ::= ¡statem¿* ”while” ¡bexpr¿ ”:” ¡statems¿ ”end”

Figure 3. Grammar describing the syntax of probabilistic loops
$\langle\textit{loop}\rangle$.

The syntax defining our program model is given by the grammar in Figure 3.
Throughout the paper, we will use the phrases _(probabilistic) loops_ and
_(probabilistic) programs_ interchangeably for loops adhering to the syntax in
Figure 3. In our work, we infer higher moments ${\mathbb{E}}(x_{n}^{k})$ of
program variables $x$ parameterized by the loop counter $n$. We abstract from
concrete loop guards by defining the guards of programs in our program model
to be _true_ (written as $\star$). Guarded loops while $\phi$: $\dots$ can be
modeled as an infinite loops while $\star$: if $\phi$: $\dots$, with the limit
behaviour giving the moments after termination (cf. Section 5.1).

Our program model defined in Figure 3 contains non-nested while-loops which
are preceded by a loop-free initialization part. The loop-body and
initialization part allow for (nested) if-statements, polynomial arithmetic,
drawing from common probability distributions, and symbolic constants.
Symbolic constants can be used to represent arbitrary real numbers and are
also used for uninitialized program variables. Categorical expressions
(defined by the non-terminal $\langle\textit{categorical}\rangle$ in Figure 3)
are expressions of the form $v_{1}\\{p_{1}\\}\dots v_{l}\\{p_{l}\\}$ such that
$\sum p_{i}=1$. Their intended meaning is that they evaluate to $v_{i}$ with
probability $p_{i}$. The last parameter $p_{l}$ can be omitted and in that
case is set to $p_{l}:=1-\sum_{i=1}^{l-1}p_{i}$. For a program ${\mathcal{P}}$
we denote with $\text{Vars}({\mathcal{P}})$ the set of ${\mathcal{P}}$’s
variables appearing on the left-hand side of an assignment in
${\mathcal{P}}$’s loop-body. The programs of Figures 1-2 are examples of our
program model defined in Figure 3. In comparison to the _probabilistic Guarded
Command Language (pGCL)_ (Barthe et al., 2020), programs of our model contain
exactly one while-loop, no non-determinism 111Non-determinism is different
from probabilistic choice. Demonic (angelic) non-determinism is concerned with
the worst-case (best-case) behavior. For instance, a variable can be assigned
to $0$ or $1$ both with probability $\nicefrac{{1}}{{2}}$. This is different
from assigning $0$ or $1$ non-deterministically, where the probability is not
specified. but support continuous distributions and simultaneous assignments.

###  3.2. Program Semantics

In what follows we define the semantics of probabilistic programs in terms of
Markov chains on a measurable space. We then introduce the notion of
_normalized_ probabilistic loops by means of so-called _${\mathcal{P}}$
-preserving program transformations_ (Section 3.3).

###### Definition 0 (State & Run Space).

Let ${\mathcal{P}}$ be a probabilistic program with $m$ variables. We denote
by _$\text{ND}({\mathcal{P}})$ the non-probabilistic program obtained from
${\mathcal{P}}$_ by replacing every probabilistic choice $C$ in
${\mathcal{P}}$ by a non-deterministic choice over $\text{supp}(C)$. Let
$\text{States}_{\mathcal{P}}\subseteq{\mathbb{R}}^{m}$ be the set of program
states of $\text{ND}({\mathcal{P}})$ reachable from any initial state. The
_state space_ of ${\mathcal{P}}$ is the measurable space
$(\text{States}_{\mathcal{P}},\mathcal{S}_{\mathcal{P}})$, where
$\mathcal{S}_{\mathcal{P}}$ is the Borel $\sigma$-algebra on
${\mathbb{R}}^{m}$ restricted to $\text{States}_{\mathcal{P}}$. The _run
space_ of ${\mathcal{P}}$ is the sequence space
$(\text{States}_{\mathcal{P}}^{\omega},\mathcal{S}_{\mathcal{P}}^{\omega})=:(\text{Runs}_{\mathcal{P}},\mathcal{R}_{\mathcal{P}})$.

In what follows, we omit the subscript ${\mathcal{P}}$ whenever the program
${\mathcal{P}}$ is irrelevant or clear from the context. Executions/runs of a
probabilistic program ${\mathcal{P}}$ define a stochastic process, as follows.

###### Definition 0 (Run Process).

Let ${\mathcal{P}}$ be a probabilistic program with $m$ variables. The _run
process_ $\Phi_{n}:\text{Runs}\to\text{States}$ is a stochastic process in the
run space mapping a program run to its $n$th state, that means,
$\Phi_{n}(run):=run_{n}$.

For program variable $x$ with index $i\geq 1$, we denote by $x_{n}$ _the
projection of $\Phi_{n}$ to its $i$th component_ $\Phi_{n}(\cdot)(i)$. Given
an arithmetic expression $A$ over ${\mathcal{P}}$’s variables, we write
$A_{n}$ for the stochastic process where every program variable $x$ in $A$ is
replaced by $x_{n}$.

###### Remark 0.

Given an initial distribution of program states $\mu$ and a Markov kernel $p$
defined according to the _standard meaning_ of the program statements, by
_Kolmogorov’s Extension Theorem_ we conclude that there is a unique
probability measure ${\mathbb{P}}_{\mathcal{P}}$ on
$(\text{Runs}_{\mathcal{P}},\mathcal{R}_{\mathcal{P}})$ such that the run
process is a Markov chain.
$(\text{Runs}_{\mathcal{P}},\mathcal{R}_{\mathcal{P}},{\mathbb{P}}_{\mathcal{P}})$
is the probability space associated to program ${\mathcal{P}}$. Distributions
and (higher) moments of ${\mathcal{P}}$’s variables are to be understood with
respect to this probability space.

For probabilistic loops according to the syntax in Figure 3, the initial
distribution $\mu$ of values of loop variables is the distribution of states
after the statements $\langle\textit{statem}\rangle^{*}$ just before the
while-loop. Moreover, the loop body in Figure 3 is considered to be atomic,
meaning the Markov kernel $p$ describes the transition between full iterations
in contrast to single statements.

###  3.3. ${\mathcal{P}}$-Preserving Transformations

Our probabilistic programs defined by the grammar in Figure 3 support rich
arithmetic and complex probabilistic behavior/distributions. Such an
expressivity of Figure 3 comes at the cost of turning the analysis of programs
defined by Figure 3 cumbersome. In this section, we address this difficulty
and introduce a number of program transformations that allow us to simplify
our probabilistic programs to a so-called _normal form_ while preserving the
joint distribution of program variables. Normal forms allow us to extract
recursive properties from the program, which we will later use to compute
moments for program variables (Section 5).

#### Schemas and Unification.

The program transformations we introduce in this section build on the notion
of _schemas_ and program parts. A _program part_ is an empty word or any word
resulting from any non-terminal of the grammar in Figure 3. For our purposes,
a schema $S$ is a program part with some subtrees in the program part’s syntax
tree being replaced by placeholder symbols $\dot{s_{1}},\dots,\dot{s_{l}}$. A
_substitution_ is a finite mapping $\sigma=\\{\dot{s_{1}}\mapsto
p_{1},\dots,\dot{s_{l}}\mapsto p_{l}\\}$ where $p_{1},\dots,p_{l}$ are program
parts. We denote by $S[\sigma]$ the program part resulting from $S$ by
replacing every $\dot{s_{i}}$ by $p_{i}$, assuming $S[\sigma]$ is well-formed.
For two schemas $S_{1}$ and $S_{2}$ a substitution $u$ such that
$S_{1}[u]=S_{2}[u]$ is called a _unifier_ (with respect to $S_{1}$ and
$S_{2}$). In this case $S_{1}$ and $S_{2}$ are called _unifiable_ (by $u$).

#### Transformations.

In what follows, we consider ${\mathcal{P}}$ to be a fixed probabilistic
program defined by Figure 3 and give all definitions relative to
${\mathcal{P}}$. A _transformation_ $T$ is a mapping from program parts to
program parts with respect to a schema Old. $T$ is _applicable_ to a
subprogram $S$ of ${\mathcal{P}}$ if $S$ and Old are unifiable by the unifier
$u$. Then, the transformed subprogram is defined as $T(S):=\textit{New}[u]$
where New is a schema depending on Old and $u$. A transformation is fully
specified by defining how New results from Old and $u$. We write
$T({\mathcal{P}},S)$ for the program resulting from ${\mathcal{P}}$ by
replacing the subprogram $S$ of ${\mathcal{P}}$ by $T(S)$.

The first transformation we consider removes simultaneous assignments from
${\mathcal{P}}$. For this, we store a copy of each assignment in an auxiliary
variable to preserve the values used for simultaneous assignments, in case an
assigned variable appears in an assignment expression. Variables are then
assigned their intended value.

###### Definition 0 (Simultaneous Assignment Transformation).

A _simultaneous assignment transformation_ is the transformation defined by

$\dot{x_{1}}$, $\dots$, $\dot{x_{l}}$ = $\dot{v_{1}}$, $\dots$, $\dot{v_{l}}$
$\mapsto$ $t_{1}$=$\dot{v_{1}}$; $\dots$; $t_{l}$=$\dot{v_{l}}$;
$\dot{x_{1}}$=$t_{1}$; $\dots$; $\dot{x_{l}}$=$t_{l}$,
where $t_{1},\dots,t_{l}$ are fresh variables.

In what follows, we assume that parameters of common distributions used in
programs are constant. Nevertheless, the following transformation enables the
use of some non-constant distribution parameters.

###### Definition 0 (Distribution Transformation).

A _distribution transformation_ is a transformation defined by either of the
mappings

  * •

$\dot{x}$ = Normal($\dot{p}$,$\dot{v}$) 2.75cm$\text{\ }\mapsto\text{\ }$ t =
Normal($0$,$\dot{v}$); $\dot{x}$=$\dot{p}$+t

  * •

$\dot{x}$ = Uniform($\dot{p_{1}}$,$\dot{p_{2}}$) 2.75cm$\text{\
}\mapsto\text{\ }$ t = Uniform($0$,$1$);
$\dot{x}$=$\dot{p_{1}}$+($\dot{p_{2}}$-$\dot{p_{1}}$)*t

  * •

$\dot{x}$ = Laplace($\dot{p}$,$\dot{b}$) 2.75cm$\text{\ }\mapsto\text{\ }$ t =
Laplace($0$,$\dot{b}$); $\dot{x}$=$\dot{p}$+t

  * •

$\dot{x}$ = Exponential($\dot{c}$/$\dot{p}$) 2.75cm$\text{\ }\mapsto\text{\ }$
t = Exponential($\dot{c}$); $\dot{x}$=$\dot{p}$*t

where, for every mapping, $t$ is a fresh variable.

###### Example 0.

Consider Figure 1. The simultaneous assignment l,g =
Laplace(x+y,$1$),Normal($0$,$1$) can be transformed using the transformation
rules from Definitions 3-4 as follows:

$\overset{\scriptscriptstyle\mathrm{(sim)}}{\text{\ }\mapsto\text{\ }}$

⬇

t1=Laplace(x+y,1)

t2=Normal(0,1)

x=t1

y=t2

$\overset{\scriptscriptstyle\mathrm{(dist)}}{\text{\ }\mapsto\text{\ }}$

⬇

t3=Laplace(0,1)

t1=x+y+t3

t2=Normal(0,1)

x=t1

y=t2

To simplify the structure of probabilistic loops, we assume else if branches
to be syntactic sugar for nested if else statements. We remove else by
splitting it into if-statements (if $C$ and if not $C$). Since variables in
$C$ could be changed within the first branch, we store their original values
in auxiliary variables and use those for the condition $C^{\prime}$ of the
second if statement. We capture this transformation in the following
definition.

###### Definition 0 (Else Transformation).

An _else transformation_ is the transformation

⬇

if $\dot{C}$: $\dddot{Branch_{1}}$

else: $\dddot{Branch_{2}}$ end

, u $\text{\ }\mapsto\text{\ }$

⬇

$t_{1}$=$x_{1}$;$\dots$;$t_{l}$=$x_{l}$

if $\dot{C}$: $\dddot{Branch_{1}}$ end

if not $C^{\prime}$: $\dddot{Branch_{2}}$ end

where $x_{1},\dots,x_{l}$ are all variables appearing in $\dot{C}[u]$ which
are also being assigned in $\dddot{Branch_{1}}[u]$. Every $t_{i}$ is a fresh
variable and $C^{\prime}$ results from $\dot{C}[u]$ by substituting every
$x_{i}$ with $t_{i}$.

To further simplify the loop body into a flattened list of assignments, we
equip every assignment $a$ of form “x = value”with a condition $C_{a}$
(initialized to true $\top$) and a default variable $d_{a}$ (initialized to
$x$), written as “x = value [$C_{a}$] $d_{a}$”. The semantics of the
conditioned assignment is that $x$ is assigned $value$ if $C_{a}$ holds just
before the assignment and $d_{a}$ otherwise. With conditioned assignments, the
loop body’s structure can be flattened using the following transformation.

###### Definition 0 (If Transformation).

An _if transformation_ is the transformation defined by

⬇

if $\dot{C_{1}}$:

$\dot{x}$ = $\dot{v}$ [$\dot{C_{2}}$] $\dot{x}$

$\dddot{Rest}$ end

, u $\text{\ }\mapsto\text{\ }$

⬇

t = $\dot{x}$

$\dot{x}$ = $\dot{v}$ [$\dot{C_{1}}$ and $\dot{C_{2}}$] $\dot{x}$

if $C$: $\dddot{Rest}$ end

where $t$ is a fresh variable and $C$ results from $\dot{C_{1}}[u]$ by
substituting $\dot{x}[u]$ by $t$. If $\dddot{Rest}[u]$ is empty, the line if
$C$: $\dddot{Rest}$ end is omitted from the result. If $\dot{x}[u]$ does not
appear in $\dot{C_{1}}[u]$ the line t = $\dot{x}$ is dropped.

###### Example 0.

The following program containing nested if-statements can be flattened as
follows:

⬇

if x == 1:

x = Bernoulli(1/2)

if x == 0: y = 1 end

end

$\overset{\scriptscriptstyle\mathrm{(if)}}{\text{\ }\mapsto\text{\ }}$

⬇

if x == 1:

x = Bernoulli(1/2)

y = 1 [x == 0] y

end

$\overset{\scriptscriptstyle\mathrm{(if)}}{\text{\ }\mapsto\text{\ }}$

⬇

t = x

x = Bernoulli(1/2) [t == 1] x

if t == 1:

y = 1 [x == 0] y

end

$\overset{\scriptscriptstyle\mathrm{(if)}}{\text{\ }\mapsto\text{\ }}$

⬇

t = x

x = Bernoulli(1/2) [t == 1] x

y = 1 [t == 1 $\land$ x == 0] y

To bring further simplicity to our program ${\mathcal{P}}$, we ensure for each
variable to be modified only once within the loop body. To remove duplicate
assignments we introduce new variables $x_{1},\dots,x_{l-1}$ to store
intermediate states. Assignments to other variables, in between the updates of
$x$, will be adjusted to refer to the latest $x_{i}$ instead of $x$.

###### Definition 0 (Multi-Assignment Transformation).

A _multi-assignment transformation_ is the transformation defined by

⬇

$\dot{x}$ = $\dot{v_{1}}$ [$\dot{C_{1}}$] $\dot{x}$;$\dddot{Rest_{1}}$;

$\dot{x}$ = $\dot{v_{2}}$ [$\dot{C_{2}}$] $\dot{x}$;$\dddot{Rest_{2}}$;

$\dots$; $\dot{x}$ = $\dot{v_{l}}$ [$\dot{C_{l}}$] $\dot{x}$;

, u $\text{\ }\mapsto\text{\ }$

⬇

$x_{1}$ = $\dot{v_{1}}$ [$\dot{C_{1}}$] $\dot{x}$;$Rest_{1}$;

$x_{2}$ = $v_{2}$ [$C_{2}$] $x_{1}$;$Rest_{2}$;

$\dots$; $\dot{x}$ = $v_{l}$ [$C_{l}$] $x_{l-1}$;

where $x_{1},\dots,x_{l-1}$ are fresh variables. For $i\geq 2$, $v_{i}$,
$C_{i}$ and $Rest_{i}$ result from $\dot{v_{i}}[u]$, $\dot{C_{i}}[u]$ and
$\dddot{Rest}_{i}[u]$, respectively, by replacing $\dot{x}[u]$ by $x_{i-1}$.

###### Example 0.

In the program of Figure 2, program line if x1 == x3: t1 = 1 else: t1 = 0 can
be transformed using transformation rules from Definitions 6-9 as follows:

$\overset{\scriptscriptstyle\mathrm{(else)}}{\text{\ }\mapsto\text{\ }}$

⬇

if x1 == x3: t1 = 1

if x1 != x3: t1 = 0

$\overset{\scriptscriptstyle\mathrm{(if)}}{\text{\ }\mapsto\text{\ }}$

⬇

t1 = 1 [x1 == x3] t1

t1 = 0 [x1 != x3] t1

$\overset{\scriptscriptstyle\mathrm{(multi)}}{\text{\ }\mapsto\text{\ }}$

⬇

t11 = 1 [x1 == x3] t1

t1 = 0 [x1 != x3] t11

With program transformations defined, we can turn our attention to program
properties. In particular, we show that our transformations of ${\mathcal{P}}$
do not change the joint distribution of ${\mathcal{P}}$’s variables. Since our
transformations may introduce new variables, we consider program equivalence
with respect to program variables in order to ensure that the distribution of
${\mathcal{P}}$ is maintained/preserved by our transformations.

###### Definition 0 (Program Equivalence).

Let ${\mathcal{P}}_{1}$ and ${\mathcal{P}}_{2}$ be two probabilistic programs.
We define ${\mathcal{P}}_{1}$ and ${\mathcal{P}}_{2}$ to be _equivalent with
respect to a set of program variables $X$_, in symbols
${\mathcal{P}}_{1}\equiv^{X}{\mathcal{P}}_{2}$, if:

  1. (1)

$X\subseteq\text{Vars}({\mathcal{P}}_{1})\cap\text{Vars}({\mathcal{P}}_{2})$,
and

  2. (2)

the joint distributions of $X$ arising from ${\mathcal{P}}_{1}$ and
${\mathcal{P}}_{2}$ are equal.

To relate a program ${\mathcal{P}}$ to its transformed version
$T({\mathcal{P}},S)$ we consider the distribution of variables of the original
program ${\mathcal{P}}$. If ${\mathcal{P}}$ retains the joint distribution of
its variables after applying transformation $T$, we say that $T$ is
${\mathcal{P}}$-preserving.

###### Definition 0 (${\mathcal{P}}$-Preserving Transformation).

We say that a transformation $T$ is _${\mathcal{P}}$ -preserving_ if
${\mathcal{P}}\equiv^{\text{Vars}({\mathcal{P}})}T({\mathcal{P}},S)$ for all
subprograms $S$ of ${\mathcal{P}}$ which are unifiable with $Old$.

It is not hard to argue that the transformations defined above are
${\mathcal{P}}$-preserving, yielding the following result.

###### Lemma 13.

The transformations from Definitions 3-9 are ${\mathcal{P}}$-preserving.

By exhaustively applying the ${\mathcal{P}}$-preserving transformations of
Definitions 3-9 over ${\mathcal{P}}$, we obtain a so-called normalized program
${\mathcal{P}}_{N}$, as defined below. The normalized ${\mathcal{P}}_{N}$ will
then further be used in computing higher moments of ${\mathcal{P}}$ in
Sections 5, as the ${\mathcal{P}}_{N}$ preserves the moments of
${\mathcal{P}}$ (Theorem 15).

###### Definition 0 (Normal Form).

A program ${\mathcal{P}}$ is in _normal form_ or a _normalized_ program if
none of the transformations from Definitions 3-9 are applicable to
${\mathcal{P}}$.

###### Theorem 15 (Normal Form).

For every probabilistic program ${\mathcal{P}}$ there is a
${\mathcal{P}}_{\mathcal{N}}$ in normal form such that
${\mathcal{P}}\equiv^{\text{Vars}({\mathcal{P}})}{\mathcal{P}}_{\mathcal{N}}$.
Moreover, ${\mathcal{P}}_{\mathcal{N}}$ can be effectively computed from
${\mathcal{P}}$ by exhaustively applying transformations from Definitions 3-9.

###### Proof.

There are two claims in the theorem, which we need to address: (i)
exhaustively applying transformation rules terminates (_termination_), and
(ii) it preserves statistical properties of the (original) program variables
(_correctness_).

For termination, we show that programs become smaller, in some sense, after
every transformation. In particular, we consider the program size to be given
by a tuple (Sim, Dist, Else, If, $\textit{Multi}_{B}$, $\textit{Multi}_{I}$),
representing the number of simultaneous assignments, non-trivial
distributions, else statements, assignments within if branches (weighted for
nested ifs), and number of variables with multiple assignments in the loop
body and initialization part, respectively. With respect to the lexicographic
order, each transformation reduces the size of the program which is lower-
bounded by $0$.

Correctness can be shown by treating each transformation separately and
showing that it does not alter the variables’ distributions after a single
application (Lemma 13). This is true for all transformations from Definitions
3-9. Auxiliary variables are used to store the original value to prevent
intervening variable modifications. For Definition 0 (Multi-Assignment
Transformation). (Definition 9), we also revise the rest of the assignments to
reflect the change of the original variable. The Definition 0 (Distribution
Transformation). (Definition 4) uses statistical properties of well-known
distributions. ∎

#### Properties of Normalized Programs.

Figure 4 shows a normal form for the program from Figure 1. Normalized
programs have the following important properties: (1) all distribution
parameters are constant; (2) the loop body is a sequence of guarded
assignments; (3) every program variable is only assigned once in the loop
body. Moreover for every guarded assignment v =
$\textit{assign}_{\textit{true}}$ [C] $\textit{assign}_{\textit{false}}$, the
guard $C$ is a boolean condition and $\textit{assign}_{\textit{false}}$ is a
single variable which is assigned to v if $C$ evaluates to _false_. If $C$
evaluates to _true_ , the variable v is assigned
$\textit{assign}_{\textit{true}}$. The expression
$\textit{assign}_{\textit{true}}$ is either a distribution with constant
parameters or a probabilistic choice of polynomials as illustrated in Figure
4.

###### Remark 0.

Based on the order in which transformations are applied to a program
${\mathcal{P}}$ and the names used for auxiliary variables, several different
normalized programs can be achieved for ${\mathcal{P}}$. In this work, only
the existence of a normal form is relevant. Moreover, from the definitions of
the transformation, it is apparent that exhaustively applying them leads to a
normal form whose size is linear in the size of the original program.

⬇

toggle = 0; sum = s0

x = 1; y = 1; z = 1

while $\star$:

toggle = 1-toggle

x = 1+x {1/2} 2+x [toggle==0] x

y = y+z+x**2 {1/3} -x+y-z [toggle==0] y

z = z+y {1/4} z-y [toggle==0] z

t1 = Laplace(0, 1)

l = t1+x+y

g = Normal(0, 1)

sum = sum+x [g < 1/2] sum

end

Figure 4. A normal form for the program in Figure 1.

##  4\. Finite Types in Probabilistic Programs

Given a probabilistic program ${\mathcal{P}}$ in our programming model, the
transformations of Section 3.3 simplify ${\mathcal{P}}$ by computing its
normalized form while maintaining the distribution (and hence also moments) of
${\mathcal{P}}$. Nevertheless, the normalized form of ${\mathcal{P}}$ contains
computationally expensive polynomial arithmetic, potentially hindering the
automated analysis of ${\mathcal{P}}$ in Section 5 due to a computational
blowup. Therefore, we introduce further simplifications for ${\mathcal{P}}$ by
means of power reduction techniques.

###### Example 0.

Consider Figure 1 and assume we are interested in the $k$th power of variable
toggle and deriving the raw moment ${\mathbb{E}}(\textit{toggle}_{n}^{k})$.
Our analysis relies on replacing variables with their assignments (see Section
5), leading to the expression $(1-\textit{toggle}_{n-1})^{k}$. When expanded,
this is a polynomial in $\textit{toggle}_{n-1}$ with $k{+}1$ monomials:

$(1-\textit{toggle}_{n-1})^{k}=\sum_{i=0}^{k}\binom{k}{i}(-1)^{i}\textit{toggle}_{n-1}^{i}$

Higher moments, together with the aforementioned replacements, may lead to
blowups of the number of monomials to consider. However, observing that the
variable toggle is binary, we have
$\textit{toggle}_{n}^{k}=\textit{toggle}_{n}=1-\textit{toggle}_{n-1}$ for any
$k\geq 0$. Arbitrary powers of the finite variable toggle with $2$ possible
values can be written in terms of powers smaller than $2$. In the rest of this
section, we show that this phenomenon generalizes from binary variables to
arbitrary finite valued variables, thus simplifying the analysis of higher
moments of finite valued program variables.

As defined in Definition 2, for an arithmetic expression $X$ over program
variables, $X_{n}$ denotes the stochastic process mapping a program run to the
value of $X$ after iteration $n$.

###### Definition 0 (Finite Expression).

Let ${\mathcal{P}}$ be a probabilistic program and $X$ an arithmetic
expression over the variables of ${\mathcal{P}}$. We say that $X$ is _finite_
if there exist $a_{1},\dots,a_{m}\in{\mathbb{R}}$ such that for all
$n\in{\mathbb{N}}$ : $X_{n}\in\\{a_{1},\dots,a_{m}\\}$.

###  4.1. Power Reduction for Finite Types

As established in (Bartocci et al., 2020a, Lemma 1), high powers $k$ of a
random variable $X$ over a finite set can be reduced. We adapt their result to
our setting as follows.

###### Theorem 3 (Finite Power Reduction).

Let $m,k\in{\mathbb{Z}}$ and $X$ be a discrete random variable over
$A=\\{a_{1},\dots,a_{m}\\}$. Then we can rewrite $X^{k}$ as a linear
combination of $1,X,X^{2},\cdots,X^{m-1}$. Furthermore,

(3) $X^{k}=\overline{a^{k}}M^{-1}\overline{X},$

where $\overline{a^{k}}=(a_{1}^{k},\dots,a_{m}^{k})$, $M$ is an $m\times m$
matrix with $M_{ij}=a_{j}^{i-1}$ (with $0^{0}:=1$), and
$\overline{X}=(X^{0},\dots,X^{m-1})^{T}$.

In other words, Theorem 3 implies that any higher moment of $X$, can be
computed from just its first $m{-}1$ moments. Furthermore, we build on Theorem
3 and establish the inverse of matrix $M$ explicitly ($M$ is explicit in
Theorem 3 but its inverse is implicit).

###### Theorem 4 (Reduction Formula).

Recall that the $k$th elementary symmetric polynomial with respect to a set
$V=\\{v_{1},\dots,v_{n}\\}$ is $\displaystyle e_{k}(V)=\sum_{1\leq
j_{1}<\cdots<j_{k}\leq n}v_{j_{1}}\cdots v_{j_{k}}$ and let
$A_{-j}=A\setminus\\{a_{j}\\})$. Then the inverse of $M$ in (3) is given by

(4) $M^{-1}_{ij}=-\frac{(-1)^{j}e_{m-j}(A_{-i})}{\prod_{a\in
A_{-i}}(a-a_{i})}.$

###### Proof.

Let $MN=B$ for $M$ as of (3) and $N$ as of (4). We show that $B=I$ by showing
that $B_{ij}=1$ if $i=j$ and $B_{ij}=0$ otherwise. We have

$\displaystyle\begin{split}B_{ij}=\sum_{1\leq k\leq
m}N_{ik}M_{kj}&=\sum_{1\leq k\leq
m}-\frac{(-1)^{k}e_{m-k}(A_{-i})}{\prod_{a\in A_{-i}}(a-a_{i})}a_{j}^{k-1}\\\
&=\frac{1}{\prod_{a\in A_{-i}}(a-a_{i})}\sum_{1\leq k\leq
m}-(-1)^{k}a_{j}^{k-1}e_{m-k}(A_{-i})\\\ &=\frac{1}{\prod_{a\in
A_{-i}}(a-a_{i})}\sum_{0\leq k\leq m-1}(-a_{j})^{k}e_{m-k-1}(A_{-i})\\\
&=\frac{1}{\prod_{a\in A_{-i}}(a-a_{i})}\prod_{a\in
A_{-i}}(a-a_{j}),\end{split}$

where the last equation comes from the expansion of product $\prod_{a\in
A_{-i}}(a-a_{j})$ and grouping by the exponent of $a_{j}$. We can clearly see
that the last expression is $1$ if $i=j$ and $0$ otherwise. ∎

###### Example 0.

Let $X$ be a random variable over $A:=\\{-2,0,1,3\\}$. Using Theorem 3-4, we
obtain the $10$th power of $X$ as:

$\displaystyle\begin{split}X^{10}&=\begin{pmatrix}(-2)^{10}&0^{10}&1^{10}&3^{10}\end{pmatrix}\begin{pmatrix}0&\nicefrac{{-1}}{{10}}&\nicefrac{{2}}{{5}}&\nicefrac{{-1}}{{30}}\\\
1&\nicefrac{{-5}}{{6}}&\nicefrac{{-1}}{{3}}&\nicefrac{{1}}{{6}}\\\
0&1&\nicefrac{{1}}{{6}}&\nicefrac{{-1}}{{6}}\\\
0&\nicefrac{{-1}}{{15}}&\nicefrac{{1}}{{30}}&\nicefrac{{1}}{{30}}\end{pmatrix}\begin{pmatrix}X^{0}\\\
X^{1}\\\ X^{2}\\\ X^{3}\end{pmatrix}\\\
&=1934X^{3}+2105X^{2}-4038X.\end{split}$

##  5\. Computing Higher Moments of Probabilistic Programs

We now bring together the results from Sections 3-4 to develop the _theory of
moment-computability_ for probabilistic loops. We establish the technical
details leading to sufficient conditions that ensure moment-computability,
culminating in the proof of Theorem 3. The main ideas of our method are
illustrated on the probabilistic loop from Figure 1 in Example 6 at the end of
this section.

###### Definition 0 (Moment-Computability).

A probabilistic loop ${\mathcal{P}}$ is _moment-computable_ if a closed-form
(according to Theorem 3) of ${\mathbb{E}}(x_{n}^{k})$ exists and is computable
for all $x\in\text{Vars}({\mathcal{P}})$ and $k\in{\mathbb{N}}$.

We will describe the class of moment-computable probabilistic loops through
the properties of the dependencies between program variables.

###### Definition 0 (Variable Dependency).

Let ${\mathcal{P}}$ be a probabilistic loop and
$x,y\in\text{Vars}({\mathcal{P}})$. We define:

  * •

$y$ _depends conditionally on_ $x$, if there is an assignment of $y$ within an
if-else-statement and $x$ appears in the if-condition.

  * •

$y$ _depends finitely on_ $x$, if $x$ is finite and appears in an assignment
of $y$.

  * •

$y$ _depends linearly on_ $x$, if $x$ appears only linearly in every
assignment of $y$.

  * •

$y$ _depends polynomially on_ $x$, if there is an assignment of $y$ in which
$x$ appears non-linearly and $x$ is not finite (motivated by Section 4.1).

  * •

$y$ _depends on_ $x$ if it depends on $x$ conditionally, finitely, linearly,
or polynomially.

Furthermore, we consider the transitive closure for variable dependency as
follows: If $z$ depends on $y$ and $y$ depends on $x$, then $z$ depends on
$x$. If one of the two dependencies is polynomial, then $z$ depends
polynomially on $x$.

A crucial point to highlight in Definition 2 is that due to transitivity,
variables can depend on themselves. For instance, if variable $x$ depends on
$y$ and $y$ on $x$, then $x$ is self-dependent. Moreover, if either of the
dependencies between $x$ and $y$ is non-linear, $x$ depends, by Definition 2,
_polynomially on itself._ The absence of such polynomial self-dependencies is
a central condition for our notion of moment-computable loops.

###### Theorem 3 (Moment-Computability).

A probabilistic loop ${\mathcal{P}}$ is moment-computable if _(1) none of its
non-finite variables depends on itself polynomially, and (2) if the variables
in all if-conditions are finite._

Note that none of the program transformations from Section 3.3 can introduce a
polynomial (self-)dependence. We capture this in the following lemma:

###### Lemma 4 (Non-Dependency Preservation).

If a variable $x\in\text{Vars}({\mathcal{P}})$ does not depend on itself
polynomially, neither does $x\in\text{Vars}({\mathcal{P}}_{\mathcal{N}})$.

Before we prove Theorem 3, let us first show its validity for programs in
normal form (as defined in Section 3.3). Recall that a normalized program’s
loop body is a flat list of (guarded) assignments, one for every (possibly
auxiliary) program variable.

###### Lemma 5 (Normal Moment-Computability).

The Moment-Computability Theorem (Theorem 3) holds for loops in normal form.

###### Proof.

We have to show that for an arbitrary normalized program ${\mathcal{P}}$
satisfying the conditions of Theorem 3, all $x\in\text{Vars}({\mathcal{P}})$
and all $k\in{\mathbb{N}}$, the $k$th moment of $x$ (that is
${\mathbb{E}}(x_{n}^{k})$) admits a closed-form as an exponential polynomial.
${\mathbb{E}}(x_{n}^{k})$ admits a closed-form as an exponential polynomial in
$n$ if it satisfies a linear recurrence. We show a slightly more general
statement. That is, we show that for _any monomial of program variables_ $M$
(and hence also for $x^{k}$) ${\mathbb{E}}(M)$ satisfies a linear recurrence.
The main idea of the proof is to show that ${\mathbb{E}}(M)$ only depends on a
finite set of monomials, each (in some sense) _not larger_ than $M$ itself.
Intuitively, the finite set of monomials on which ${\mathbb{E}}(M)$ depends on
are all monomials of program variables such that their expected values
determine ${\mathbb{E}}(M)$. We will show that this set of monomials exists,
is finite, and leads to a system of linear recurrences containing
${\mathbb{E}}(M)$, implying a computable exponential polynomial closed-form
for ${\mathbb{E}}(M)$ by Theorem 3.

Let ${\mathcal{P}}$ be a normalized program satisfying the conditions of
Theorem 3, $x\in\text{Vars}({\mathcal{P}})$, $k,n\in{\mathbb{N}}$ arbitrary,
and $\mathcal{M}$ the set of all monomials over $\text{Vars}({\mathcal{P}})$
with the powers of every finite variable $d$ bounded by the number of possible
values of $d$ (higher powers can be reduced as of Theorem 3).

#### Recurrences over Moments.

Given the syntax of probabilistic programs and properties of expectation, for
any monomial $M\in\mathcal{M}$ there is a natural way to express the expected
value of $M$ in iteration $n{+}1$, that is ${\mathbb{E}}(M_{n+1})$, as a
linear combination of expectations of monomials in iteration $n$:

(5) ${\mathbb{E}}(M_{n+1})=\sum_{N\in M^{*}}c_{N}{\mathbb{E}}(N_{n}),$

for some finite set $M^{*}\subset\mathcal{M}$, and non-zero constants $c_{N}$.
Equation (5) is called _the recurrence_ of ${\mathbb{E}}(M)$. The set $M^{*}$
is the set of monomials that appear in the recurrence of ${\mathbb{E}}(M)$. We
define the ∗ operator to give such a set for any monomial and extend the
definition to sets by

$S^{*}=\bigcup_{M\in S}M^{*}.$

The exact recurrence can be computed from ${\mathbb{E}}(M_{n+1})$ by replacing
variables appearing in $M$ by their assignments and using the linearity of
${\mathbb{E}}$ to convert expected values of polynomials to linear
combinations of expected monomials. Recall that for a program in normal form,
every program variable is only assigned once, all distribution parameters are
constant, and the loop body is a flat list of guarded assignments (Section
3.3). The guarded assignments are of the form

$x=a_{0}\\{p_{0}\\}\dots\\{p_{i-1}\\}a_{i}[C_{x}]d_{x},$

for polynomials of program variables $a_{0},\dots,a_{i}$ and constant
probabilities $p_{0},\dots,p_{i-1}$, or

$x=Dist[C_{x}]d_{x}$

for some admissible distribution $Dist$. The guard $[C_{x}]$ is a boolean
condition and for normalized programs, $d_{x}$ is always a program variable.
Assume, that the variable $x$ appears in the monomial $M$. Hence,
$M=M^{\prime}\cdot x^{k}_{n+1}$ for some monomial $M^{\prime}$ not containing
$x$. If the single assignment of $x$ is a (guarded) probabilistic choice of
polynomials we rewrite ${\mathbb{E}}(M_{n+1})$ to

(6) ${\mathbb{E}}(M_{n+1})={\mathbb{E}}(M^{\prime}\cdot
x^{k}_{n+1})={\mathbb{E}}\left(M^{\prime}\left(d_{x}[\lnot C_{x}]+\sum
p_{i}a_{i}^{k}[C_{x}]\right)\right).$

If the single assignment of $x$ is a (guarded) draw from a distribution we
rewrite ${\mathbb{E}}(M_{n+1})$ to

(7) ${\mathbb{E}}(M_{n+1})={\mathbb{E}}(M^{\prime}\cdot
x^{k}_{n+1})={\mathbb{E}}\left(M^{\prime}d_{x}[\lnot
C_{x}]\right)+{\mathbb{E}}\left(M^{\prime}[C_{x}]\right){\mathbb{E}}\left(Dist^{k}\right).$

Variables in conditions $[C_{x}]$ are all finite since they come from branch
conditions. We further simplify the expressions of Equations (6)-(7) by
replacing the logical conditions $[C_{x}]$ by polynomials that evaluate to $1$
whenever variables satisfy the condition $[C_{x}]$ and to $0$ otherwise. It is
possible to write any logical condition over finitely valued variables as such
a polynomial ((Stankovic et al., 2022)), with
$[x=c]:=\prod_{d\in\omega(x)\setminus\\{c\\}}\frac{x-d}{c-d}$, $[\lnot
C]=1-[C]$, and $[C_{1}\land C_{2}]=[C_{1}]\cdot[C_{2}]$, where $\omega(x)$ is
the set of possible values of $x$222Because negation and conjunction are
functionally complete for propositional logic, we can also handle disjunctions
using De Morgan’s laws: $[P\lor Q]=[\lnot(\lnot P\land\lnot
Q)]=1-(1-[P])\cdot(1-[Q])$. Inequalities can then be transformed into a
disjunction of equalities since we assume that all variables appearing in if-
conditions are finitely valued.. Converting conditions to polynomials in the
equation above leads to polynomials over moments of program variables. We can
compute the recurrence of ${\mathbb{E}}(M)$ in Equation (5) by replacing all
variables in $M$ of iteration $n{+}1$ from last to first (by their appearance
in $\mathcal{P}$’s loop body). Throughout, the linearity of ${\mathbb{E}}$ is
used to convert expected values of polynomials to linear combinations of
expected monomials. In Example 6, we illustrate this computation on the
program from Figure 1 with its normal form from Figure 4.

#### Ordering.

Now that we can compute the recurrences, we need to introduce the order for
monomials, such that the monomials appearing in the recurrence for $M$ are not
larger than $M$. We will need this to show that computing the recurrences as
described above is, indeed, a finite process. Intuitively, a variable $y$ is
larger than (or equal to) $x$ if it depends on $x$. Mutually dependent
variables will form an equivalence class. We then extend the order to
monomials based on their degrees with respect to the variables’ equivalence
classes.

More formally, let $\preceq$ be a smallest total preorder on variables such
that $x\preceq y$ whenever $y$ depends on $x$. We write $x\prec y$ iff
$x\preceq y$ and $y\not\preceq x$, and $x\sim y$ iff $x\preceq y$ and
$y\preceq x$. Let $\widetilde{x}$ be the equivalence class of $x$ induced by
$\sim$. Note, that all variables in an equivalence class are mutually
dependent. However, because of our assumptions on the program ${\mathcal{P}}$,
the mutual dependencies among non-finite variables are all linear (as there
are no polynomial self-dependencies).

We extend $\preceq$ to a preorder on the set of monomials $\mathcal{M}$. For
every monomial $M$ and non-finite variable $x$, we consider the degree
$deg(\widetilde{x},M)$ of $M$ in the equivalence class of
$\widetilde{x}$333$deg(\widetilde{x},M):=\sum_{x\in\widetilde{x}}deg(x,M)$,
where $deg(x,M)$ is the degree of $x$ in $M$.. We associate $M$ with the
sequence of $deg(\widetilde{x},M)$ for equivalence classes of all non-finite
variables, ordered reverse-lexicographically with respect to $\preceq$. Then
the relations $\prec$, $\sim$, and equivalence classes $\widetilde{(-)}$
follow naturally from $\preceq$.

By Theorem 3 and the definition of $\preceq$, the equivalence class
$\widetilde{M}$ is finite for each $M\in\mathcal{M}$. Let $\mathcal{M}_{\sim}$
be the set of equivalence classes of $\sim$. The preorder $\preceq$ induces a
partial order $\preceq_{\sim}$ on $\mathcal{M}_{\sim}$. Note that monomials
only contain non-negative powers and a finite number of variables. These facts
together with $\preceq$ being total imply that $\preceq_{\sim}$ is a well-
order. We will write $\preceq$ instead of $\preceq_{\sim}$ when the meaning is
clear from the context.

With these orders defined, we have $N\preceq M$ for any $N\in M^{*}$ and
$M\in\mathcal{M}$. Intuitively, this means that for every monomial $M$, the
monomials occurring in the recurrence of $M$ are not larger than $M$ itself.
This is true because the order on variables is defined according to variable
dependencies, the order on monomials is a reverse-lexicographic extension of
the order on variables, and the fact that polynomial self-dependencies are not
allowed by assumption.

#### Finite $\mathcal{S}^{(Q)}$

We show, that for any monomial $Q$, there is a finite set of monomials
$\mathcal{S}^{(Q)}\subset\mathcal{M}$ containing $Q$ such that
$M^{*}\subset\mathcal{S}^{(Q)}$ for any $M\in\mathcal{S}^{(Q)}$. Because
$M^{*}$ is the set of all monomials in the recurrence of $M$, this means that
$\mathcal{S}^{(Q)}$ contains all monomials necessary to construct a system of
linear recurrences containing $E(Q)$ (if it exists and is finite).

Let

(8)
$\mathcal{S}^{(Q)}:=\widetilde{Q}\cup\bigcup_{\begin{subarray}{c}A\in\widetilde{Q}^{*}\\\
A\prec Q\end{subarray}}\mathcal{S}^{(A)}.$

Clearly $Q\in\mathcal{S}^{(Q)}$ and $M^{*}\subset\mathcal{S}^{(Q)}$ for any
$M\in\mathcal{S}^{(Q)}$ by construction. We are left to show that
$\mathcal{S}^{(Q)}$ is finite for all
$Q\in\widetilde{Q}\in\mathcal{M}_{\sim}$, which we can do by transfinite
induction over $\mathcal{M}_{\sim}$.

For the base case, we have to show that $\mathcal{S}^{(Q)}$ is finite for all
$Q$ in $\widetilde{1}$ (the trivial monomial). Since
$\widetilde{1}=\\{\prod_{d\in F}d^{\lambda_{d}}\mid\lambda_{d}\leq m_{d}\\}$,
where $F$ is the set of finite program variables and $m_{d}$ are upper bounds
on their powers as of Theorem 3, $\mathcal{S}^{(Q)}=\widetilde{1}$ is finite
for all $Q\in\widetilde{1}$. Suppose $\mathcal{S}^{(A)}$ is finite for all
$A\prec Q$. Since $\widetilde{Q}^{*}$ is finite, so is the union in (8) and,
as a result, $\mathcal{S}^{(Q)}$.

#### Moments.

A system of linear recurrences with constant coefficients can be constructed
for monomials in $\mathcal{S}^{(Q)}$. Therefore, the closed-form of any
$E(M)\in\mathcal{M}$ exists and is computable. ∎

We now turn back to Theorem 3 and establish its validity. The crux of our
proof below relies on the fact that our transformations computing normal forms
of ${\mathcal{P}}$ (see Section 3.3) are ${\mathcal{P}}$-preserving.

###### Proof (of Theorem 3).

By Theorem 15 (Theorem 15 (Normal Form). Termination), ${\mathcal{P}}$ can be
transformed to a normalized loop ${\mathcal{P}}_{\mathcal{N}}$. By Lemma 4
(Lemma 4 (Non-Dependency Preservation).), ${\mathcal{P}}_{\mathcal{N}}$
satisfies all conditions from Lemma 5 (Lemma 5 (Normal Moment-
Computability).). Thus, ${\mathcal{P}}_{\mathcal{N}}$ is moment-computable and
the moments are equivalent to those of ${\mathcal{P}}$ for
$x\in\text{Vars}({\mathcal{P}})$ by Theorem 15 (Theorem 15 (Normal Form).
Correctness). ∎

The proofs of Theorem 3 and Lemma 5 are constructive and describe a procedure
to compute (higher) moments of program variables by (1) transforming a
probabilistic loop into a normal form according to Theorem 15, (2)
constructing a system of linear recurrences as in the proof of Lemma 5 and (3)
solving the system of linear recurrences with constant coefficients. In the
following example, we illustrate the whole procedure on the running example
from Figure 1.

###### Example 0.

We return to the probabilistic loop ${\mathcal{P}}$ from Figure 1. A normal
form ${\mathcal{P}}_{\mathcal{N}}$ for ${\mathcal{P}}$ was given in Figure 4.
To compute a closed-form of the expected value of the program variable $z$, we
will model ${\mathbb{E}}(z_{n})$ as a system of linear recurrences according
to the proof of Lemma 5. For a cleaner presentation we will refer to the
variable _toggle_ by $t$. Note that $t$ is binary. To construct the recurrence
for ${\mathbb{E}}(z_{n})$, we start with ${\mathbb{E}}(z_{n+1})$ at the
assignment of $z$ in ${\mathcal{P}}_{\mathcal{N}}$ and repeatedly replace
variables by the right-hand side of their assignment starting from $z$’s
assignment and stopping at the top of the loop body. Throughout the process we
use the linearity of expectation to convert expected values of polynomials to
linear combinations of expected monomials (as required by Equation 5):

$\displaystyle{\mathbb{E}}(z_{n+1})$
$\displaystyle\quad\downarrow\text{assignment of $z$}$
$\displaystyle{\mathbb{E}}(([t_{n+1}=0](z_{n}+y_{n+1})\
\\{\nicefrac{{1}}{{4}}\\}\ [t_{n+1}=0](z_{n}-y_{n+1}))+[t_{n+1}\neq 0]z_{n})$
$\displaystyle\quad\downarrow\text{replace (in)equalities by polynomials}$
$\displaystyle{\mathbb{E}}(((1-t_{n+1})(z_{n}+y_{n+1})\
\\{\nicefrac{{1}}{{4}}\\}\ (1-t_{n+1})(z_{n}-y_{n+1}))+t_{n+1}z_{n})$
$\displaystyle\quad\downarrow\text{${\mathbb{E}}$ on categorical}$
$\displaystyle\frac{1}{4}{\mathbb{E}}((1-t_{n+1})(z_{n}+y_{n+1}))+\frac{3}{4}{\mathbb{E}}((1-t_{n+1})(z_{n}-y_{n+1}))+{\mathbb{E}}(t_{n+1}z_{n})$
$\displaystyle\quad\downarrow\text{simplify; ${\mathbb{E}}$ linearity}$
$\displaystyle{\mathbb{E}}(z_{n})+\frac{1}{2}{\mathbb{E}}(t_{n+1}y_{n+1})-\frac{1}{2}{\mathbb{E}}(y_{n+1})$
$\displaystyle\quad\downarrow_{*}\text{similarly, replace
$y_{n+1},x_{n+1},t_{n+1}$}$
$\displaystyle{\mathbb{E}}(z_{n})-\frac{1}{6}{\mathbb{E}}(t_{n}x_{n})-\frac{1}{2}{\mathbb{E}}(t_{n}y_{n})-\frac{1}{6}{\mathbb{E}}(t_{n}x_{n}^{2})+\frac{1}{12}{\mathbb{E}}(t_{n})+\frac{1}{6}{\mathbb{E}}(t_{n}z_{n}).$

The last line of the calculation represents the recurrence equation of the
expected value of $z$. For every monomial in the recurrence equation of
${\mathbb{E}}(z_{n})$ (that is: $z^{*}=$ $\\{tx,ty,tx^{2},t,tz\\}$), we
compute its respective recurrence equation and recursively repeat this
procedure which eventually terminates according to Lemma 5.

The ordering of variables and monomials, which is essential for proving
termination in Lemma 5 comes from an ordering of variables $t,x,y,z,l,g,sum$.
Based on the program assignments, we need $t\preceq x\prec y\preceq l;y\preceq
z;t\preceq z\preceq y;x\preceq sum;g\preceq sum$. Variables y and z form an
equivalence class, since $y\preceq z$ and $z\preceq y$. Notice that some
variables are not ordered, for example $t$ and $g$. This gives some freedom in
choosing the total preorder, any will work. Let us have $t\prec g\prec x\prec
sum\prec y\sim z\prec l$. This gives an order on the equivalence classes
$\\{t\\}\prec\\{g\\}\prec\\{x\\}\prec\\{sum\\}\prec\\{y,z\\}\prec\\{l\\}$. The
ordering on monomials then considers the classes of non-finite program
variables, i.e. all equivalence classes except $\\{t\\}$. A monomial is then
assigned a sequence of degrees with respect to each equivalence class, e.g.
$z\to(0,0,0,1,0)$, $tz\to(0,0,0,1,0)$, and $x^{2}\to(0,0,2,0,0)$, with
monomials ordered reverse-lexicographically with respect to their sequence.
For the monomials in this example: $1\sim t\prec x\sim tx\prec x^{2}\sim
tx^{2}\prec y\sim z\sim ty\sim tz$.

After computing all necessary recurrence equations, we are faced with a system
of linear recurrences of the expected values of the monomials $z$, $tx$, $ty$,
$x$, $y$, $tx^{2}$, $t$, $tz$, $x^{2}$ and $1$ with recurrence matrix

$\footnotesize\left[\begin{matrix}1&-\frac{1}{6}&-\frac{1}{2}&0&0&-\frac{1}{6}&\frac{1}{12}&\frac{1}{6}&0&0\\\
0&-1&0&1&0&0&0&0&0&0\\\ 0&0&-1&0&1&0&0&0&0&0\\\
0&0&0&1&0&0&\frac{3}{2}&0&0&0\\\
0&\frac{1}{3}&0&0&1&\frac{1}{3}&-\frac{1}{6}&-\frac{1}{3}&0&0\\\
0&0&0&0&0&-1&0&0&1&0\\\ 0&0&0&0&0&0&-1&0&0&1\\\ 1&0&0&0&0&0&0&-1&0&0\\\
0&3&0&0&0&0&\frac{5}{2}&0&1&0\\\ 0&0&0&0&0&0&0&0&0&1\end{matrix}\right].$

Using any computer algebra system, such as sympy444https://www.sympy.org, we
arrive at the closed-form solution for the expected value of $z$,
${\mathbb{E}}(z_{n})$, parameterized by the loop counter $n$:

$\displaystyle\footnotesize\begin{split}{\mathbb{E}}(z_{n})&=\frac{883}{32}+\frac{29n}{16}-\frac{201}{20}2^{-n}6^{\frac{n}{2}}-\frac{67}{20}2^{-n}6^{\frac{1+n}{2}}-\frac{37}{10}3^{-n}6^{\frac{n}{2}}\\\
&-\frac{37}{20}3^{-n}6^{\frac{1+n}{2}}-\frac{201}{20}6^{\frac{n}{2}}\left(\frac{-1}{2}\right)^{n}-\frac{37}{10}6^{\frac{n}{2}}\left(\frac{-1}{3}\right)^{n}+\frac{67}{20}6^{\frac{1+n}{2}}\left(\frac{-1}{2}\right)^{n}\\\
&+\frac{37}{20}6^{\frac{1+n}{2}}\left(\frac{-1}{3}\right)^{n}+\frac{9}{16}n\left(-1\right)^{n}+\frac{29}{32}\left(-1\right)^{n}+\frac{9}{16}n^{2}.\end{split}$

This example highlights the strength of algebraic techniques for probabilistic
program analysis in comparison to constraint-based methods employing
templates. We are not aware of any template-based method able to handle
functions of the complexity of ${\mathbb{E}}(z_{n})$. Our tool Polar is able
to find the closed-form of ${\mathbb{E}}(z_{n})$ in under one second (see
Section 7).

###  5.1. Guarded Loops

When we model a guarded probabilistic loop while $\phi$: $\dots$ as an
infinite loop while $\star$: if $\phi$: $\dots$, we impose the same
restrictions on $\phi$ as on if-conditions (that means $\phi$ only contains
finite variables) to guarantee computability and correctness. The $k$th moment
of $x$ after termination is then given by

(9)
$\lim_{n\to\infty}{\mathbb{E}}(x_{n}^{k}\mid\lnot\phi_{n})=\lim_{n\to\infty}\frac{{\mathbb{E}}(x_{n}^{k}\cdot[\lnot\phi_{n}])}{{\mathbb{E}}([\lnot\phi_{n}])}.$

If the limit exists, we can use standard methods from computer algebra to
compute it, as the (higher) moments our approach computes are given as
exponential polynomials (Gruntz, 1996).

###### Example 0.

Consider the following loop, in which $x$ after termination is geometrically
distributed with parameter $\nicefrac{{1}}{{2}}$:

⬇

x, stop = 0, 0

while stop == 0:

stop=Bernoulli(1/2)

x=x+1

end

With Equation 9 and the techniques from this section we get:

$\displaystyle\footnotesize{\mathbb{E}}(x)=\lim_{n\to\infty}{\mathbb{E}}(x_{n}\mid
stop_{n}=1)=\lim_{n\to\infty}\frac{{\mathbb{E}}(x_{n}\cdot
stop_{n})}{{\mathbb{E}}(stop_{n})}=\lim_{n\to\infty}\frac{-n2^{-n}-2^{1-n}+2}{1-2^{-n}}=2.$

Moreover, whenever the variables in the loop guard are not probabilistic,
traditional techniques can be applied to determine the number of loop
iterations n which can then be plugged into the (higher) moments computed by
our work. Apart from the guarded loops, many systems show the type of infinite
behavior naturally modeled with infinite loops, such as probabilistic
protocols or dynamical systems.

###  5.2. Infinite If-Conditions

Theorem 3 on moment-computability requires the variables in all if-conditions
to be finite. Nevertheless, in some cases, if-conditions containing infinite
variables can be handled by our approach. Let ${\mathcal{P}}$ be a
probabilistic loop containing an if-statement with condition $F$ and $I$ where
$F$ contains only finite variables and $I$ contains infinite variables.
Without loss of generality, no variable in $I$ is assigned in or after the if-
statement. Let the transformation removing $I$ be defined by  ⬇ if $F$ and
$I$:$Branch$ $\text{\ }\mapsto\text{\ }$ ⬇ t = Bernoulli($p$) if $F$ and t ==
1: $Branch$ end where _t_ is a fresh variable and $p:={\mathbb{P}}(I)$
(potentially symbolic). Then, the transformation preserves the distributions
of all $x\in\text{Vars}({\mathcal{P}})$ under the following assumptions:

  1. (1)

$I$ is iteration independent, meaning for every variable $x$ in $I$ neither
$x$ nor any variable $x$ depends on (as of Definition 2) has a self-
dependency.

  2. (2)

$I$ is statistically independent from $F$ and all conditions in $Branch$.

  3. (3)

For every assignment $A$ in $Branch$ and every variable $x$ in $A$ which has
been assigned before A, it holds that $I$ and $x$ are statistically
independent.

Assumption 1 ensures that ${\mathbb{P}}(I)={\mathbb{E}}([I])$ is constant.
Assumption 2 and 3 further ensure that ${\mathbb{E}}([I])$ can always be
“pulled out” (that means
${\mathbb{E}}([I]x)={\mathbb{E}}([I]){\mathbb{E}}(x)$) in the construction of
the recurrences. Assumptions 1-3 can often be checked automatically.

###### Example 0.

Consider the statement if g < 1/2: sum=sum+x of the program from Figure 1,
where the value of $g$ is drawn from a standard normal distribution. In this
case, the transformation’s parameter $p$ represents
${\mathbb{P}}({\text{Normal}}(0,1)<\nicefrac{{1}}{{2}})$, but is left symbolic
for the moment computation. The integral
${\mathbb{P}}({\text{Normal}}(0,1)<\nicefrac{{1}}{{2}})$ can be solved
separately and the result be substituted for $p$.

###  5.3. On the Necessity of the Conditions Ensuring Moment-Computability

Theorem 3 states two conditions that are sufficient to ensure that the closed-
forms of the program variables’ higher moments always exist and are
computable. _Condition 1 enforces that there is no variable with potentially
infinite values with a polynomial self-dependency._ _Condition 2 demands that
all variables appearing in if-conditions are finite._ Our approach for
computing the moments of variables of probabilistic loops can handle precisely
the programs that satisfy these two conditions. We argue that both conditions
are necessary, in the sense that if either of the conditions does not hold,
the existence or computability of the variable moments’ closed-forms as
exponential polynomials cannot be guaranteed for all programs from our program
model when one or both conditions are removed from Theorem 3.

#### Condition 1.

Relaxing condition 1 of Theorem 3 means that we allow for polynomial self-
dependencies of non-finite variables. The _logistic map_ (May, 1976) is a
quadratic first-order recurrence defined by $x_{n+1}=r\cdot x_{n}(1-x_{n})$
and well-known for its chaotic behavior. A famous fact about the logistic map
is that it does not have an analytical solution for most values of $r$
(Maritz, 2020). By neglecting condition 1, we can easily devise a loop
modeling the logistic map:

⬇

while $\star$:

x = r$\cdot$x(1-x)

end

The value of the program variable $x$ after iteration $n$ is equal to the
$n$th term of the logistic map. This means, for most values of $r$ and initial
values of $x$, there does not exist an analytical closed-form solution for the
program variable $x$. Moreover, our counter-example illustrates that condition
1 is necessary already for programs with a single variable and without
stochasticity and if-statements.

#### Condition 2.

Loosening condition 2 of Theorem 3 and allowing for non-finite variables in
if-conditions renders our programming model Turing-complete. Intuitively, one
can model a Turing machine’s tape with two variables l and r such that the
binary representation of l represents the tape’s content left of the read-
write-head. The binary representation of r represents the tape’s content at
the position of the read-write-head and towards the right. The least
significant bit of r is the current symbol the Turing machine is reading. We
can extract the least significant bit of r in our programming model (and
neglecting condition 2) by introducing a variable lsb and using a single if-
statement involving non-finite variables: whenever the loop changes the value
of r, we set $\texttt{lsb}:=\texttt{r}$. The while-loop’s body is of the form
“if lsb > 1: lsb=lsb-2 else $\mathit{transitions}$ end”. The Turing machine’s
transition table can be encoded using if-statements. Writing and shifting can
be accommodated for by multiplying by $2$ or $\nicefrac{{1}}{{2}}$ and using
addition and subtraction. The Turing machine’s state can be modelled by a
single finite variable. Therefore, by dropping condition 2, being able to
model the program variables by linear recurrences would give rise to a
decision procedure for the Halting problem: assume we introduce a variable
terminated which is initialized to $0$ before the loop and set to $1$ whenever
the Turing-machine terminates. If terminated can be modelled by a linear
recurrence of order $k$, it suffices to check the first $k$ values of the
recurrence to determine whether or not terminated is always $0$ (Kauers and
Paule, 2011) and the Turing-machine does not terminate. As the Halting problem
is well-known to be undecidable, condition 2 is necessary to guarantee that
the program variables can be modelled by linear recurrences, even without
stochasticity and polynomial arithmetic.

###### Remark 0 (Sequential & Nested Loops).

Our program model consists of single non-nested loops. Sequential loops can be
analyzed one by one with the same techniques as presented in this section. For
nested loops, one could design a program transformation transforming a nested
loop into a non-nested loop and then apply the techniques presented in this
section. Alternatively, we conjecture that the approach presented for guarded
loops (Section 5.1) could be used to first compute the moments of the most
inner loops and then use the obtained information to compute the moments of
the outer loops. The main challenge lies in ensuring the moment-computability
conditions for the outer loops (Theorem 3) once the inner loops have been
analyzed.

##  6\. Use-Cases of Higher Moments

For probabilistic loops, computing closed-forms of the variables’ (higher)
moments poses a technique for synthesizing quantitative invariants: Given a
program variable $x$ and a closed-form $f(n)$ of its $k$th moment, the
equation ${\mathbb{E}}(x_{n}^{k})-f(n)=0$ is an invariant. Moreover, closed-
forms of raw moments can be converted into closed-forms of _central_ moments,
such as variance, skewness or kurtosis (cf. Section 2.1). In addition, this
section provides hints on two further use-cases of higher moments of
probabilistic loops: (i) deriving tail probabilities (Section 6.1) and (ii)
inferring distributions of random variables from their moments (Section 6.2).

###  6.1. From Moments to Tail Probabilities

Tail probabilities measure the probability that a random variable surpasses
some value. The mathematical literature contains several inequalities
providing upper- or lower bounds on tail probabilities given (higher) moments
(Boucheron et al., 2013). Two examples are _Markov’s inequality_ for upper and
the _Paley-Zygmund inequality_ for lower bounds.

###### Theorem 1 (Markov’s Inequality).

Let $X$ be a non-negative random variable, and $t\geq 0$, then

${\mathbb{P}}(X\geq t)\leq\frac{{\mathbb{E}}(X^{k})}{t^{k}}.$

###### Theorem 2 (Paley-Zygmund Inequality).

Let $X$ be a random variable with $X\geq t$ almost-surely. Then

${\mathbb{P}}(X>t)\geq\frac{({\mathbb{E}}(X)-t)^{2}}{{\mathbb{E}}(X^{2})-2t{\mathbb{E}}(X)+t^{2}}.$

###### Example 0.

For _Herman’s Self-Stabilization_ program from Figure 2 almost-surely
$tokens\in\\{0,1,2,3\\}$. With the techniques from previous sections, we can
compute the first two moments ${\mathbb{E}}(tokens_{n})=1+2\cdot 4^{-n}$ and
${\mathbb{E}}(tokens^{2}_{n})=1+8\cdot 4^{-n}$. Markov’s inequality (Theorem
1) gives us the upper bound ${\mathbb{P}}(tokens_{n}\geq
2)\leq\nicefrac{{1}}{{2}}+4^{-n}$ using the first moment and
${\mathbb{P}}(tokens_{n}\geq 2)\leq\nicefrac{{1}}{{4}}+2\cdot 4^{-n}$
utilizing the second moment.

For the Paley-Zygmund inequality (Theorem 2) both the first and the second
moment are required, yielding the lower bound ${\mathbb{P}}(tokens_{n}\geq
2)={\mathbb{P}}(tokens_{n}>1)\geq 4^{-n}$. The theorem’s precondition that
almost-surely $tokens\geq 1$ might not be apparent at first sight. We take
this for granted for now and will clarify this fact in Example 5.

Markov’s inequality and the Paley-Zygmund inequality are just two examples
showing that our technique for moment computation can be leveraged for further
program analysis using known results from probability theory. Our approach
computes the _exact moments_ of variables in probabilistic loops instead of
just approximations or bounds on moments. This enables our technique to be
readily combined with results from probability theory that require exact
moments.

###  6.2. From Moments to Distributions

For finite random variables, their full distribution can be recovered from
finitely many moments. More precisely, given a random variable $X$ with $m$
possible values, the distribution of $X$ can be recovered from its first
$m{-}1$ moments, as the following theorem states:

###### Theorem 4.

Let $X$ be a random variable over $\\{a_{1},\dots,a_{m}\\}$ and
$p_{i}:=P(X=a_{i})$. The values $p_{i}$ are the solutions of the system of
linear equations given by $\sum_{i=1}^{m}p_{i}a_{i}^{j}={\mathbb{E}}(X^{j})$
for $0\leq j<m$.

###### Proof.

There are $m$ unknowns $p_{i}$ for $1\leq p_{i}\leq m$. Note that all $a_{i}$
are constant and that the first $m{-}1$ raw moments of $X$ are fixed. Using
the definition of raw moments, we get $m$ linear equations
$\sum_{i=1}^{m}p_{i}a_{i}^{j}={\mathbb{E}}(X^{j})$ for $0\leq j<m$. The
solutions of the system of $m$ linear equations are the values $p_{i}$ that
determine the distribution of $X$. ∎

###### Example 0.

Consider _Herman’s Self-Stabilization_ program from Figure 2. In Example 3 we
obtained upper and lower bounds for tail probabilities of the $tokens$
variable using the first one or two moments. With the first _three_ moments we
can fully recover the distribution of the $tokens$ variable. We have that
$tokens\in\\{0,1,2,3\\}$. Let $p_{i}:={\mathbb{P}}(tokens_{n}=i)$ for $0\leq
i\leq 3$. By Theorem 4, we get the following system:

$\displaystyle\begin{split}p_{0}+p_{1}+p_{2}+p_{3}&={\mathbb{E}}(tokens_{n}^{0})=1,\\\
p_{1}+2p_{2}+3p_{3}&={\mathbb{E}}(tokens_{n})=1+2\cdot 4^{-n},\\\
p_{1}+4p_{2}+9p_{3}&={\mathbb{E}}(tokens_{n}^{2})=1+8\cdot 4^{-n},\\\
p_{1}+8p_{2}+27p_{3}&={\mathbb{E}}(tokens_{n}^{3})=1+26\cdot
4^{-n}.\end{split}$

The solution can be obtained using standard techniques and tools, yielding
$p_{0}=0;p_{1}=1-4^{-n};p_{2}=0;p_{3}=4^{-n}$.

Note that probabilities are given as functions of the loop iteration $n$.
Moreover, the solution shows that almost-surely $tokens\geq 1$, which we
assumed to be true in Example 3.

Figure 5.  The empirical density of program variable $x$ for the benchmark
_Bimodal_ (cf. Table 1) and loop iterations $10$ (left) and $100$ (right)
obtained by $10^{5}$ samples, together with two approximations using the Gram-
Charlier A Series with $6$ (red dashed lines) and $12$ (blue solid lines)
moments.

The distributions of program variables with potentially infinitely many
values, including continuous variables, cannot be, in general, fully
reconstructed from finitely many moments. However, expansions such as the
_Gram–Charlier A Series_ (Kolassa, 2006) can be used to approximate a
probability density function using finitely many moments. Figure 5 illustrates
how exact moments computed by our approach can be used to approximate unknown
probability density functions of program variables. While Figure 5 shows
approximations for specific loop iterations, we emphasize that the symbolic
nature of our approach allows for approximating the densities of program
variables for all – potentially infinitely many – loop iterations
simultaneously. Therefore, our technique is constant in the number of loop
iterations, whereas sampling has linear complexity. We compute the
approximations from Figure 5 for all infinitely many loop iterations in
$\sim$$22$ seconds, with the experimental setup from Section 7. In comparison,
sampling the loop $10^{5}$ times takes $\sim$$3.6$ minutes for loop iteration
$10$ and $\sim$$33$ minutes for loop iteration $100$.

##  7\. Implementation and Evaluation

#### Implementation.

The program transformations (Section 3) and (higher) moment computation
(Section 5) are implemented in the new tool Polar (Moosbrugger et al., 2022).
The experiments can be reproduced using the corresponding
artifact555https://doi.org/10.5281/zenodo.7055030. For automatically inferring
finiteness of program variables, we use a standard approach based on _abstract
interpretation_. Polar is implemented in python3, consisting of $\sim$$3300$
LoC, and uses the packages sympy666https://www.sympy.org and
symengine777https://github.com/symengine for symbolic manipulation of
mathematical expressions. Together with all our benchmarks, Polar is publicly
available at https://github.com/probing-lab/polar.

#### Experimental Setting and Evaluation.

The evaluation of our work is split into three parts. First, we evaluate Polar
on the ability of computing higher moments for $15$ probabilistic programs
exhibiting different characteristics (Section 7.1). Second, we compare Polar
to the exact tool Mora (Bartocci et al., 2020b) which computes so-called
_moment-based invariants_ for a subset of our programming model (Section 7.2).
Third, we compare our tool to approximate methods estimating program variable
moments by confidence intervals through sampling (Section 7.3). All
experiments have been run on a machine with a
$2.6\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ Intel i7 (Gen 10) processor and
$32\text{\,}\mathrm{G}\mathrm{B}$ of RAM. Runtime measurements are averaged
over $10$ executions.

Table 1. Evaluation of Polar on $15$ benchmarks. All times are in seconds. #V = number of variables in benchmark; C = benchmark contains circular dependencies; If = benchmark contains if-statements; S = benchmark contains symbolic constants; INF = benchmark’s states space is infinite; CONT = benchmark’s state space is continuous; Moment = Moment to compute; RT = Total runtime. Benchmark | #V | C/If/S/INF/CONT | Moment | RT
---|---|---|---|---
Running-Example (Fig. 1) | 7 | ✓/ ✓/ ✓/ ✓/ ✓ | ${\mathbb{E}}(z)$ | 0.67
Herman-3 | 10 | ✓/ ✓/ ✗/ ✗/ ✗ | ${\mathbb{E}}(\text{tokens}^{3})$ | 0.58
Las-Vegas-Search | 3 | ✗/ ✓/ ✗/ ✓/ ✗ | ${\mathbb{E}}(\text{found}^{20})$ | 0.36
Pi-Approximation | 4 | ✗/ ✓/ ✗/ ✓/ ✓ | ${\mathbb{E}}(\text{count}^{3})$ | 0.47
50-Coin-Flips | 101 | ✗/ ✓/ ✗/ ✗/ ✗ | ${\mathbb{E}}(\text{total})$ | 0.91
Gambler-Ruin-Momentum | 4 | ✓/ ✗/ ✓/ ✓/ ✗ | ${\mathbb{E}}(x^{3})$ | 2.89
Hawk-Dove-Symbolic | 5 | ✗/ ✓/ ✓/ ✓/ ✗ | ${\mathbb{E}}(\text{p1bal}^{4})$ | 2.00
Variable-Swap | 4 | ✓/ ✗/ ✗/ ✓/ ✓ | ${\mathbb{E}}(x^{30})$ | 2.42
Retransmission-Protocol | 4 | ✗/ ✓/ ✓/ ✓/ ✗ | ${\mathbb{E}}(\text{fail}^{3})$ | 1.62
Randomized-Response | 7 | ✗/ ✓/ ✓/ ✓/ ✗ | ${\mathbb{E}}(\text{p1}^{3})$ | 0.59
Duelling-Cowboys | 4 | ✓/ ✓/ ✓/ ✗/ ✗ | ${\mathbb{E}}(\text{ahit})$ | 1.14
Martingale-Bet | 4 | ✗/ ✓/ ✓/ ✓/ ✗ | ${\mathbb{E}}(\text{capital}^{3})$ | 8.44
Bimodal | 5 | ✗/ ✓/ ✗/ ✓/ ✓ | ${\mathbb{E}}(\text{x}^{10})$ | 4.50
DBN-Umbrella | 2 | ✗/ ✓/ ✓/ ✗/ ✗ | ${\mathbb{E}}(\text{umbrella}^{5})$ | 0.77
DBN-Component-Health | 3 | ✗/ ✓/ ✗/ ✗/ ✗ | ${\mathbb{E}}(\text{obs}^{5})$ | 0.26

###  7.1. Experimental Results with Higher Moments

Table 1 shows the evaluation of Polar on the program from Figure 1 and $14$
benchmarks which are either from the literature on probabilistic programming
(Kwiatkowska et al., 2012) (_Herman-3_), (McIver and Morgan, 2005) (_Duelling-
Cowboys_), (Gretz et al., 2013) (_Martingale-Bet_), (Chakarov and
Sankaranarayanan, 2014) (_Hawk-Dove-Symbolic_ , _Variable-Swap_), (Barthe et
al., 2016) (_Gambler-Ruin-Momentum_), (Batz et al., 2021) (_Retransmission-
Protocol_), differential privacy schemes (Warner, 1965) (_Randomized-
Response_), Dynamic Bayesian Networks (_DBN-Umbrella_ , _DBN-Component-
Health_) or well-known stochastic processes (_Las-Vegas-Search_ , _Pi-
Approximation_ , _Bimodal_). The benchmarks _Retransmission-Protocol_ and
_Hawk-Dove-Symbolic_ were further generalized from their original definition
by replacing concrete numbers with symbolic constants. This makes these
benchmarks only harder as solutions to the generalized versions are solutions
for the concretizations. Table 1 illustrates that Polar can compute higher
moments for various probabilistic programs exhibiting different features, like
circular variable dependencies, if-statements, and symbolic constants with
finite, infinite, continuous, and discrete state spaces. Moreover, the table
shows that the number of program variables is _not_ the primary factor for the
complexity of computing moments. For instance, the benchmarks _50-Coin-Flips_
and _Duelling-Cowboys_ have $101$ and $4$ program variables respectively.
Nevertheless, the runtimes for computing first moments for the two benchmarks
only differ by $0.23\text{\,}\mathrm{s}$. The complexity of computing moments
lies in the complexity of the resulting systems of recurrences which depend on
the concrete features present in the benchmarks like specific variable
dependencies, symbolic constants, or degrees of polynomials.

###  7.2. Experimental Comparison to Exact Methods

Table 2. Comparison of Polar to Mora. The runtimes are in seconds per tool,
benchmark, and moment. For Polar the comparison contains in brackets the
seconds spent on parsing, normalizing, and type inference.

Benchmark | Mora | Polar
---|---|---
COUPON  | ${\mathbb{E}}(c)$
---
${\mathbb{E}}(c^{2})$
${\mathbb{E}}(c^{3})$
| $0.25$
---
$0.27$
$0.29$
| $0.29$
---
$0.29$
$0.29$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
COUPON4  | ${\mathbb{E}}(c)$
---
${\mathbb{E}}(c^{2})$
${\mathbb{E}}(c^{3})$
| $0.71$
---
$0.87$
$1.22$
| $0.36$
---
$0.36$
$0.36$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.08)$
RANDOM_WALK_1D  | ${\mathbb{E}}(x)$
---
${\mathbb{E}}(x^{2})$
${\mathbb{E}}(x^{3})$
| $0.07$
---
$0.11$
$0.10$
| $0.12$
---
$0.23$
$0.24$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
SUM_RND_SERIES  | ${\mathbb{E}}(x)$
---
${\mathbb{E}}(x^{2})$
${\mathbb{E}}(x^{3})$
| $0.27$
---
$0.97$
$2.48$
| $0.27$
---
$0.43$
$0.79$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
PRODUCT_DEP_VAR  | ${\mathbb{E}}(p)$
---
${\mathbb{E}}(p^{2})$
${\mathbb{E}}(p^{3})$
| $0.37$
---
$1.41$
$4.03$
| $0.28$
---
$0.46$
$0.97$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.08)$
RANDOM_WALK_2D  | ${\mathbb{E}}(x)$
---
${\mathbb{E}}(x^{2})$
${\mathbb{E}}(x^{3})$
| $0.10$
---
$0.21$
$0.17$
| $0.12$
---
$0.24$
$0.24$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.08)$
BINOMIAL(p)  | ${\mathbb{E}}(x)$
---
${\mathbb{E}}(x^{2})$
${\mathbb{E}}(x^{3})$
| $0.12$
---
$0.32$
$0.79$
| $0.25$
---
$0.29$
$0.44$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$

Benchmark | Mora | Polar
---|---|---
STUTTERING_A  | ${\mathbb{E}}(s)$
---
${\mathbb{E}}(s^{2})$
${\mathbb{E}}(s^{3})$
| $0.29$
---
$1.13$
$3.32$
| $0.26$
---
$0.43$
$0.97$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
STUTTERING_B  | ${\mathbb{E}}(s)$
---
${\mathbb{E}}(s^{2})$
${\mathbb{E}}(s^{3})$
| $0.26$
---
$0.94$
$2.26$
| $0.27$
---
$0.39$
$0.79$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
STUTTERING_C  | ${\mathbb{E}}(s)$
---
${\mathbb{E}}(s^{2})$
${\mathbb{E}}(s^{3})$
| $0.76$
---
$12.43$
$74.83$
| $0.40$
---
$1.94$
$8.19$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.08)$
STUTTERING_D  | ${\mathbb{E}}(s)$
---
${\mathbb{E}}(s^{2})$
${\mathbb{E}}(s^{3})$
| $0.76$
---
$8.19$
$25.67$
| $0.43$
---
$1.34$
$4.33$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
STUTTERING_P  | ${\mathbb{E}}(s)$
---
${\mathbb{E}}(s^{2})$
${\mathbb{E}}(s^{3})$
| $0.26$
---
$1.17$
$3.53$
| $0.28$
---
$0.49$
$1.17$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$
SQUARE  | ${\mathbb{E}}(y)$
---
${\mathbb{E}}(y^{2})$
${\mathbb{E}}(y^{3})$
| $0.30$
---
$0.88$
$1.98$
| $0.29$
---
$0.45$
$0.65$
$\color[rgb]{.25,.25,.25}\definecolor[named]{pgfstrokecolor}{rgb}{.25,.25,.25}\pgfsys@color@gray@stroke{.25}\pgfsys@color@gray@fill{.25}(0.07)$

To the best of our knowledge, Mora is the only other tool capable of computing
higher moments for variables of probabilistic loops without templates – as
described in (Bartocci et al., 2019). Mora operates on so-called _Prob-
solvable loops_ which form a strict subset of our program model (Section 3).
Prob-solvable loops do not admit circular variable dependencies, if-
statements, or state-dependent distribution parameters. We compare Polar
against Mora on the Mora benchmarks taken from (Kura et al., 2019; Chen et
al., 2015; Chakarov and Sankaranarayanan, 2014; Katoen et al., 2010). Details
can be found in Table 2. The experiments illustrate that Polar can handle all
programs and moments that Mora can. Mora, however, cannot compute any moment
for any program in Table 1. On simple benchmarks, Polar is slightly slower
than Mora due to the constant overhead of the program transformations and type
inference to identify finite valued variables. On complex benchmarks Polar
provides a significant speedup compared to Mora. For instance, for the
_STUTTERING_C_ benchmark Polar computes the moment ${\mathbb{E}}(s^{3})$ in
about $8$ seconds, whereas Mora needs over one minute.

Table 3. Comparison of Polar to approximation through sampling. Polar = the tools runtime to compute the precise moment; CI N = an approximated 0.95-CI-interval from N samples; $T_{100.000}$ = the runtime for CI 100.000. The symbolic constant $p$ in _Retransmission-Protocol_ is set to $0.9$ and for _Hawk-Dove-Symbolic_ we set $v=4$ and $c=8$. Benchmark | Moment |  | CI 100
---
CI 1.000
CI 100.000
$T_{100.000}$ | ${\tt Polar}$
Running-Example (Fig. 1) | ${\mathbb{E}}(z_{10})$ |  | $({-}55.6,6.55)$
---
$({-}58.5,{-}39.3)$
$({-}45.6,{-}43.7)$
545.6s | 0.67s
Retransmission-Protocol | ${\mathbb{E}}(\text{fail}_{10})$ |  | $(0.09,0.25)$
---
$(0.11,0.16)$
$(0.109,0.114)$
146.2s | 0.30s
Variable-Swap | ${\mathbb{E}}(y_{10})$ |  | $(4.47,5.75)$
---
$(5.06,5.45)$
$(5.50,5.54)$
245.8s | 0.13s
Hawk-Dove-Symbolic | ${\mathbb{E}}(p1bal_{10})$ |  | $(8.28,12.3)$
---
$(9.07,10.5)$
$(9.86,10.0)$
347.4s | 0.27s

###  7.3. Experimental Comparison with Sampling

For a probabilistic loop with program variable $x$ the moment
${\mathbb{E}}(x^{k}_{n})$ can be approximated for fixed $k$ and $n$ by
sampling $x^{k}_{n}$ and calculating the sample average or confidence
intervals. Table 3 compares Polar to computing confidence intervals by
sampling for $k=1$ and $n=10$. The table shows that our tool is able to
compute precise moments in a fraction of the time needed to sample programs to
achieve satisfactory confidence intervals. An advantage of sampling is that it
is applicable for any probabilistic loop. However, by its nature, sampling
fails to give any formal guarantees or hard bounds on the approximated
moments. This is critical if the loop body contains branches that are executed
with low probability. If applicable, Polar can provide _exact_ moments for
symbolic $n$ (and involving other symbolic constants) faster than sampling can
establish acceptable approximations. Moreover, even if the sampling of the
loops is sped up by using a more efficient implementation, Polar enjoys
complexity theoretical advantages. The complexity of sampling is linear in
both the number of samples and the number of loop iterations. In contrast,
Polar does not need to take multiple samples for higher precision as it
symbolically computes the exact moments. Additionally, our method is
_constant_ in the number of loop iterations. With Polar, computing the moment
for a specific loop iteration, say $10^{5}$, just amounts to evaluate the
closed-form at $10^{5}$.

Figure 6 illustrates the importance of higher moments for probabilistic loops.
The first moment provides a center of mass but contains no information on how
the mass is distributed around this center. For this purpose higher moments
are essential.

(a) Running-Example (Fig. 1)

(b) Variable-Swap

(c) Hawk-Dove-Symbolic

Figure 6. Samples obtained by simulation plotted together with precise moments
computed by Polar. In each benchmark, the thin gray lines are $200$ samples
over $30$ iterations. The thick red line is the precise expected value. The
dotted red lines are the expected values $\pm$ twice the standard deviation
given by the precise first two moments. Figure 6(a) is symmetric log scale.

###  7.4. Evaluation Summary

Our experimental evaluation demonstrates that: (1) Polar can compute higher
moments for a rich class of probabilistic loops with various characteristics,
(2) Polar outperforms the state-of-the-art of moment computation for
probabilistic loops in terms of supported programs and efficiency, and (3)
Polar computes exact moments magnitudes faster than sampling can establish
reasonable approximations.

##  8\. Related Work

Using recurrence equations to extract closed-forms for variables and
quantitative invariants of loops is a well-studied technique for non-
probabilistic programs (Farzan and Kincaid, 2015; Breck et al., 2020; Kincaid
et al., 2019, 2018; de Oliveira et al., 2016; Humenberger et al., 2017, 2018;
Kovács, 2008; Rodríguez-carbonell and Kapur, 2004). Because a classical
program is a special case of a probabilistic program, our technique presented
in Section 5 is a generalization of the closed-form computation for classical
programs to probabilistic programs. Moreover, the generalization to
probabilistic programs is not trivial. One reason for this is that for
classical programs the closed form for $x^{p}$ is just the closed-form for x
to the power $p$. However, this fails for moments, as in general $E(x^{p})$ is
not equal to $E(x)^{p}$.

A common approach to quantitatively and exactly analyze probabilistic programs
is to employ probabilistic model checking techniques (Baier and Katoen, 2008;
Kwiatkowska et al., 2011; Dehnert et al., 2017; Katoen et al., 2011; Holtzen
et al., 2021).

Exact inference for computing precise posterior distributions for
probabilistic programs has been studied in (Gehr et al., 2016; Holtzen et al.,
2020; Narayanan et al., 2016; Claret et al., 2013; Saad et al., 2021). An
interesting direction for future research is using our techniques to assist
probabilistic inference in the presence of loops.

A different approach to characterize the distributions of program variables
are statistical methods such as Monte Carlo and hypothesis testing (Younes and
Simmons, 2006). Simulations are however performed on a chosen finite number of
program steps and do not provide guarantees over a potentially infinite
execution, such as unbounded loops, limiting thus their use (if at all) for
invariant generation.

In (McIver and Morgan, 2005), a deductive approach, the _weakest pre-
expectation calculus_ , for reasoning about PPs with discrete program
variables is introduced. Based on the weakest pre-expectation calculus,
(Katoen et al., 2010) presents the first template-based approach for
generating linear quantitative invariants for PPs. Other works (Feng et al.,
2017; Chen et al., 2015) also address the synthesis of non-linear invariants
or employ _martingale_ expressions (Barthe et al., 2016). All of these works
target a slightly different problem and, unlike our approach, rely on
templates. The first data-driven technique for invariant generation for PPs is
presented in (Bao et al., 2022).

Another line of related work comes with computing bounds over expected values
(Bouissou et al., 2016; Karp, 1994; Chatterjee et al., 2020) and higher
moments (Kura et al., 2019; Wang et al., 2021). The approach in (Bouissou et
al., 2016) can provide bounds for higher moments and can handle non-linear
terms at the price of producing more conservative bounds. In contrast, our
approach natively supports probabilistic polynomial assignments and provides a
precise symbolic expression for higher moments.

The technique presented in (Bartocci et al., 2019) automates the generation of
so-called moment-based invariants for a subclass of PPs with polynomial
probabilistic updates and sets the basis for fully automatic exact higher
moment computation. Relative to our approach, (Bartocci et al., 2019) supports
neither if-statements (thus also no guarded loops), state-dependent
distribution parameters, nor circular variable dependencies. Our work
establishes stronger theoretical foundations.

##  9\. Conclusion

We describe a fully automated approach for inferring exact higher moments for
program variables of a large class of probabilistic loops with complex control
flow, polynomial assignments, symbolic constants, circular dependencies among
variables, and potentially uncountable state spaces. Our work uses program
transformations to normalize and simplify probabilistic programs while
preserving the joint distribution of program variables. We propose a power
reduction technique for finite program variables to ease the complex
polynomial arithmetic of probabilistic programs. We prove soundness and
completeness of our approach, by establishing the theory of moment-computable
probabilistic loops. We demonstrate use cases of exact higher moments in the
context of computing tail probabilities and recovering distributions from
moments. Our experimental evaluation illustrates the applicability of our
work, solving several examples whose automation so far was not yet supported
by the state-of-the-art in probabilistic program analysis.

## Data Availability Statement

The tool Polar together with all benchmarks and scripts necessary to reproduce
the results reported in this paper are available through an openly accessible
artifact (Moosbrugger et al., 2022).

###### Acknowledgements.

This research was supported by the WWTF ICT19-018 grant ProbInG, the ERC
Consolidator Grant ARTIST 101002685, the Austrian FWF project W1255-N23, and
the SecInt Doctoral College funded by TU Wien. We thank the anonymous
reviewers for their outstanding and detailed feedback.

## References

  * (1)
  * Baier and Katoen (2008) Christel Baier and Joost-Pieter Katoen. 2008.  _Principles of model checking_.  MIT Press.  ISBN 978-0-262-02649-9.
  * Bao et al. (2022) Jialu Bao, Nitesh Trivedi, Drashti Pathak, Justin Hsu, and Subhajit Roy. 2022.  Data-Driven Invariant Learning for Probabilistic Programs. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-031-13185-1_3
  * Barthe et al. (2016) Gilles Barthe, Thomas Espitau, Luis María Ferrer Fioriti, and Justin Hsu. 2016.  Synthesizing Probabilistic Invariants via Doob’s Decomposition. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-319-41528-4_3
  * Barthe et al. (2012a) Gilles Barthe, Benjamin Grégoire, and Santiago Zanella Béguelin. 2012a.  Probabilistic Relational Hoare Logics for Computer-Aided Security Proofs. In _Proc. of MPC_.  https://doi.org/10.1007/978-3-642-31113-0
  * Barthe et al. (2020) Gilles Barthe, Joost-Pieter Katoen, and Alexandra Silva. 2020.  _Foundations of Probabilistic Programming_.  Cambridge University Press.  https://doi.org/10.1017/9781108770750
  * Barthe et al. (2012b) Gilles Barthe, Boris Köpf, Federico Olmedo, and Santiago Zanella Béguelin. 2012b.  Probabilistic Relational Reasoning for Differential Privacy. In _Proc. of POPL_.  https://doi.org/10.1145/2103656.2103670
  * Bartocci et al. (2019) Ezio Bartocci, Laura Kovács, and Miroslav Stankovic. 2019\.  Automatic Generation of Moment-Based Invariants for Prob-Solvable Loops. In _Proc. of ATVA_.  https://doi.org/10.1007/978-3-030-31784-3_15
  * Bartocci et al. (2020a) Ezio Bartocci, Laura Kovács, and Miroslav Stankovic. 2020a.  Analysis of Bayesian Networks via Prob-Solvable Loops. In _Proc. of ICTAC_.  https://doi.org/10.1007/978-3-030-64276-1
  * Bartocci et al. (2020b) Ezio Bartocci, Laura Kovács, and Miroslav Stankovic. 2020b.  Mora - Automatic Generation of Moment-Based Invariants. In _Proc. of TACAS_.  https://doi.org/10.1007/978-3-030-45190-5_28
  * Batz et al. (2021) Kevin Batz, Mingshuai Chen, Benjamin Lucien Kaminski, Joost-Pieter Katoen, Christoph Matheja, and Philipp Schröer. 2021\.  Latticed k-Induction with an Application to Probabilistic Programs. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-030-81688-9_25
  * Boucheron et al. (2013) Stéphane Boucheron, Gábor Lugosi, and Pascal Massart. 2013.  _Concentration Inequalities - A Nonasymptotic Theory of Independence_.  Oxford University Press.  https://doi.org/10.1093/acprof:oso/9780199535255.001.0001
  * Bouissou et al. (2016) Olivier Bouissou, Eric Goubault, Sylvie Putot, Aleksandar Chakarov, and Sriram Sankaranarayanan. 2016\.  Uncertainty Propagation Using Probabilistic Affine Forms and Concentration of Measure Inequalities. In _Proc. of TACAS_.  https://doi.org/10.1007/978-3-662-49674-9_13
  * Breck et al. (2020) Jason Breck, John Cyphert, Zachary Kincaid, and Thomas W. Reps. 2020.  Templates and Recurrences: Better Together. In _Proc. of PLDI_.  https://doi.org/10.1145/3385412.3386035
  * Chakarov and Sankaranarayanan (2014) Aleksandar Chakarov and Sriram Sankaranarayanan. 2014.  Expectation Invariants for Probabilistic Program Loops as Fixed Points. In _Proc. of SAS_.  https://doi.org/10.1007/978-3-319-10936-7_6
  * Chatterjee et al. (2020) Krishnendu Chatterjee, Hongfei Fu, Amir Kafshdar Goharshady, and Ehsan Kafshdar Goharshady. 2020\.  Polynomial Invariant Generation for Non-deterministic Recursive Programs. In _Proc. of PLDI_.  https://doi.org/10.1145/3385412
  * Chen et al. (2015) Yu-Fang Chen, Chih-Duo Hong, Bow-Yaw Wang, and Lijun Zhang. 2015.  Counterexample-Guided Polynomial Loop Invariant Generation by Lagrange Interpolation. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-319-21690-4
  * Chou et al. (2020) Yi Chou, Hansol Yoon, and Sriram Sankaranarayanan. 2020\.  Predictive Runtime Monitoring of Vehicle Models Using Bayesian Estimation and Reachability Analysis. In _Proc. of IROS_.  https://doi.org/10.1109/IROS45743.2020.9340755
  * Claret et al. (2013) Guillaume Claret, Sriram K. Rajamani, Aditya V. Nori, Andrew D. Gordon, and Johannes Borgström. 2013.  Bayesian Inference Using Data Flow Analysis. In _Proc. of ESEC/FSE_.  https://doi.org/10.1145/2491411.2491423
  * de Oliveira et al. (2016) Steven de Oliveira, Saddek Bensalem, and Virgile Prevosto. 2016\.  Polynomial Invariants by Linear Algebra. In _Proc. of ATVA_.  https://doi.org/10.1007/978-3-319-46520-3_30
  * Dehnert et al. (2017) Christian Dehnert, Sebastian Junges, Joost-Pieter Katoen, and Matthias Volk. 2017.  A Storm is Coming: A Modern Probabilistic Model Checker. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-319-63390-9
  * Durrett (2019) Rick Durrett. 2019\.  _Probability: Theory and Examples_.  Cambridge University Press.  https://doi.org/10.1017/9781108591034
  * Farzan and Kincaid (2015) Azadeh Farzan and Zachary Kincaid. 2015.  Compositional Recurrence Analysis. In _Proc. of FMCAD_.  https://doi.org/10.1109/FMCAD.2015.7542253
  * Feng et al. (2017) Yijun Feng, Lijun Zhang, David N. Jansen, Naijun Zhan, and Bican Xia. 2017.  Finding Polynomial Loop Invariants for Probabilistic Programs. In _Proc. of ATVA_.  https://doi.org/10.1007/978-3-319-68167-2_26
  * Gehr et al. (2016) Timon Gehr, Sasa Misailovic, and Martin T. Vechev. 2016\.  PSI: Exact Symbolic Inference for Probabilistic Programs. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-319-41528-4_4
  * Ghahramani (2015) Zoubin Ghahramani. 2015\.  Probabilistic Machine Learning and Artificial Intelligence.  _Nature_ (2015).  https://doi.org/10.1038/nature14541
  * Gretz et al. (2013) Friedrich Gretz, Joost-Pieter Katoen, and Annabelle McIver. 2013.  Prinsys - On a Quest for Probabilistic Loop Invariants. In _Proc. of QEST_.  https://doi.org/10.1007/978-3-642-40196-1_17
  * Gruntz (1996) Dominik Gruntz. 1996\.  _On Computing Limits in a Symbolic Manipulation System_.  Ph. D. Dissertation. ETH Zürich.  https://doi.org/10.3929/ETHZ-A-001631582
  * Hastings (1970) Wilfred K. Hastings. 1970\.  Monte Carlo Sampling Methods Using Markov Chains and Their Applications.  _Biometrika_ (1970).  https://doi.org/10.2307/2334940
  * Herman (1990) Ted Herman. 1990\.  Probabilistic Self-Stabilization.  _Inform. Process. Lett._ (1990).  https://doi.org/10.1016/0020-0190(90)90107-9
  * Holtzen et al. (2020) Steven Holtzen, Guy Van den Broeck, and Todd D. Millstein. 2020\.  Scaling Exact Inference for Discrete Probabilistic Programs. In _Proc. of OOPSLA_.  https://doi.org/10.1145/3428208
  * Holtzen et al. (2021) Steven Holtzen, Sebastian Junges, Marcell Vazquez-Chanlatte, Todd D. Millstein, Sanjit A. Seshia, and Guy Van den Broeck. 2021.  Model Checking Finite-Horizon Markov Chains with Probabilistic Inference. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-030-81688-9_27
  * Humenberger et al. (2017) Andreas Humenberger, Maximilian Jaroschek, and Laura Kovács. 2017.  Automated Generation of Non-Linear Loop Invariants Utilizing Hypergeometric Sequences. In _Proc. of ISSAC_.  https://doi.org/10.1145/3087604.3087623
  * Humenberger et al. (2018) Andreas Humenberger, Maximilian Jaroschek, and Laura Kovács. 2018.  Invariant Generation for Multi-Path Loops with Polynomial Assignments. In _Proc. of VMCAI_.  https://doi.org/10.1007/978-3-319-73721-8_11
  * Kaminski et al. (2019) Benjamin Lucien Kaminski, Joost-Pieter Katoen, and Christoph Matheja. 2019.  On the Hardness of Analyzing Probabilistic Programs.  _Acta Informatica_ (2019).  https://doi.org/10.1007/s00236-018-0321-1
  * Kaminski et al. (2016) Benjamin Lucien Kaminski, Joost-Pieter Katoen, Christoph Matheja, and Federico Olmedo. 2016.  Weakest Precondition Reasoning for Expected Run-Times of Probabilistic Programs. In _Proc. of ESOP_.  https://doi.org/10.1007/978-3-662-49498-1
  * Karp (1994) Richard M. Karp. 1994\.  Probabilistic Recurrence Relations.  _J. ACM_ (1994).  https://doi.org/10.1145/195613.195632
  * Katoen et al. (2010) Joost-Pieter Katoen, Annabelle McIver, Larissa Meinicke, and Carroll C. Morgan. 2010.  Linear-Invariant Generation for Probabilistic Programs: Automated Support for Proof-Based Methods. In _Proc. of SAS_.  https://doi.org/10.1007/978-3-642-15769-1_24
  * Katoen et al. (2011) Joost-Pieter Katoen, Ivan S. Zapreev, Ernst Moritz Hahn, Holger Hermanns, and David N. Jansen. 2011\.  The Ins and Outs of the Probabilistic Model Checker MRMC.  _Perform. Eval._ (2011).  https://doi.org/10.1016/j.peva.2010.04.001
  * Kauers and Paule (2011) Manuel Kauers and Peter Paule. 2011.  _The Concrete Tetrahedron - Symbolic Sums, Recurrence Equations, Generating Functions, Asymptotic Estimates_.  Springer.  https://doi.org/10.1007/978-3-7091-0445-3
  * Kincaid et al. (2019) Zachary Kincaid, Jason Breck, John Cyphert, and T. Reps. 2019\.  Closed Forms for Numerical Loops. In _Proc. of POPL_.  https://doi.org/10.1145/3290368
  * Kincaid et al. (2018) Zachary Kincaid, John Cyphert, Jason Breck, and T. Reps. 2018\.  Non-Linear Reasoning for Invariant Synthesis. In _Proc. of POPL_.  https://doi.org/10.1145/3158142
  * Kolassa (2006) John E. Kolassa. 2006\.  _Series Approximation Methods in Statistics_.  Springer.  https://doi.org/10.1007/0-387-32227-2
  * Kovács (2008) Laura Kovács. 2008\.  Reasoning Algebraically About P-Solvable Loops. In _Proc. of TACAS_.  https://doi.org/10.1007/978-3-540-78800-3_18
  * Kozen (1985) Dexter Kozen. 1985\.  A Probabilistic PDL.  _J. Comput. System Sci._ (1985).  https://doi.org/10.1016/0022-0000(85)90012-1
  * Kura et al. (2019) Satoshi Kura, Natsuki Urabe, and Ichiro Hasuo. 2019\.  Tail Probabilities for Randomized Program Runtimes via Martingales for Higher Moments. In _Proc. of TACAS_.  https://doi.org/10.1007/978-3-030-17465-1_8
  * Kwiatkowska et al. (2011) Marta Z. Kwiatkowska, Gethin Norman, and David Parker. 2011\.  PRISM 4.0: Verification of Probabilistic Real-Time Systems. In _Proc. of CAV_.  https://doi.org/10.1007/978-3-642-22110-1
  * Kwiatkowska et al. (2012) Marta Z. Kwiatkowska, Gethin Norman, and David Parker. 2012\.  Probabilistic Verification of Herman’s Self-Stabilisation Algorithm.  _Formal Aspects of Computing_ (2012).  https://doi.org/10.1007/s00165-012-0227-6
  * Maritz (2020) Milton F. Maritz. 2020\.  A Note on Exact Solutions of the Logistic Map.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ (2020).  https://doi.org/10.1063/1.5125097
  * May (1976) Robert M. May. 1976\.  Simple Mathematical Models With Very Complicated Dynamics.  _Nature_ (1976).  https://doi.org/10.1038/261459a0
  * McIver and Morgan (2005) Annabelle McIver and Carroll Morgan. 2005.  _Abstraction, Refinement and Proof for Probabilistic Systems_.  Springer.  https://doi.org/10.1007/b138392
  * Moosbrugger et al. (2022) Marcel Moosbrugger, Miroslav Stankovič, Ezio Bartocci, and Laura Kovács. 2022.  _This is the Moment for Probabilistic Loops \- Artifact (Polar)_.  https://doi.org/10.5281/zenodo.7055030
  * Motwani and Raghavan (1995) Rajeev Motwani and Prabhakar Raghavan. 1995.  _Randomized Algorithms_.  Cambridge University Press.  https://doi.org/10.1017/cbo9780511814075
  * Narayanan et al. (2016) Praveena Narayanan, Jacques Carette, Wren Romano, Chung chieh Shan, and Robert Zinkov. 2016\.  Probabilistic Inference by Program Transformation in Hakaru (System Description). In _Proc. of FLOPS_.  https://doi.org/10.1007/978-3-319-29604-3_5
  * Rodríguez-carbonell and Kapur (2004) Enric Rodríguez-carbonell and Deepak Kapur. 2004.  Automatic Generation of Polynomial Loop Invariants: Algebraic Foundations. In _Proc. of ISSAC_.  https://doi.org/10.1145/1005285.1005324
  * Saad et al. (2021) Feras A. Saad, Martin C. Rinard, and Vikash K. Mansinghka. 2021\.  SPPL: Probabilistic Programming With Fast Exact Symbolic Inference. In _Proc. of PLDI_.  https://doi.org/10.1145/3453483.3454078
  * Selyunin et al. (2015) Konstantin Selyunin, Denise Ratasich, Ezio Bartocci, Md. Ariful Islam, Scott A. Smolka, and Radu Grosu. 2015.  Neural Programming: Towards Adaptive Control in Cyber-Physical Systems. In _Proc. of CDC_.  https://doi.org/10.1109/CDC.2015.7403319
  * Stankovic et al. (2022) Miroslav Stankovic, Ezio Bartocci, and Laura Kovács. 2022\.  Moment-Based Analysis of Bayesian Network Properties.  _Theor. Comput. Sci._ (2022).  https://doi.org/10.1016/j.tcs.2021.12.021
  * Wang et al. (2021) Di Wang, Jan Hoffmann, and Thomas Reps. 2021.  Central Moment Analysis for Cost Accumulators in Probabilistic Programs. In _Proc. of PLDI_.  https://doi.org/10.1145/3453483.3454062
  * Warner (1965) Stanley L. Warner. 1965\.  Randomized Response: A Survey Technique for Eliminating Evasive Answer Bias.  _J. Am. Stat. Assoc._ (1965).  https://doi.org/10.1080/01621459.1965.10480775
  * Younes and Simmons (2006) Håkan L. S. Younes and Reid G. Simmons. 2006.  Statistical probabilistic model checking with a focus on time-bounded properties.  _Inf. Comput._ (2006).  https://doi.org/10.1016/j.ic.2006.05.002