may not terminate Other semantics often use sub-probability measures to express non-termination [4,23]. Thus, these semantics would say that Listing 5 results in 0 with probability 1 2 (and nothing else). We propose to track the probability of non-termination explicitly by an additional state -, just as we track the probability of observation failure (-).

*Partial Functions.* Many functions that are practically useful are only partial (meaning they are not defined for some inputs). Examples include **uniform**(a, b) (undefined for b<a) and √x (undefined for x < 0). Listing 6 shows an example program using √x. Usually, semantics do not explicitly address partial functions [23,24,28,33] or use

```
if flip( 1
         2 ) {
  x:=0;
  observe(flip( 1
                 2 ));
}else{
  x:=1;
  observe(flip( 1
                 4 ));
}
Listing 3. The need for
```
tracking -

```
n:=0;
while !flip( 1
              2 ) {
  n=n+1;
}
return n;
```
Listing 4. Geometric distribution

```
x := 5;
while x>0 {
  if x<10 {
    x+=2*flip( 1
                2 )-1;
  }
}
return x;
```
Listing 5. Program that

```
x:=uniform(-1,1);
x=√x;
return x;
Listing 6. Using par-
tial functions
```
Fig. 1. Visual comparison of the exception handling capabilities of different semantics. For example, is filled in [34] because its semantics can handle non-termination. However, the intersection between and is not filled because [34] cannot distinguish non-termination from observation failure.

partial functions without dealing with failure (e.g. [19] use **Bernoulli**(p) without stating what happens if p /∈ [0, 1]). Most of these languages could use a subprobability distribution that misses weight in the presence of errors (in these languages, this results in conflating errors with non-termination and observation failures).

We introduce a third exception state ⊥ that can be produced when partial functions are evaluated outside of their domain. Thus, Listing 6 results in ⊥ with probability 1 2 and returns a value from [0, 1] with probability 1 2 (larger values are more likely). Some previous work uses an error state to capture failing computations, but does not propagate this failure implicitly [34,35]. In particular, if an early expression in a long program may fail evaluating √−4, every expression in the program that depends on this failing computation has to check whether an exception has occurred. While it may seem possible to skip the rest of the function in case of a failing computation (by applying the pattern **if** (x = ⊥) {**return** ⊥} **else** {rest of function}), this is non-modular and does not address the result of the function being used in other parts of a program.

Although our semantics treat ⊥ and similarly, there is an important distinction between the two: ⊥ means the program terminated due to an error, while means that according to observed evidence, the program did not actually run.

#### 2.2 Interaction of Exception States

Next, we illustrate the interaction of different exception states. We explain how our semantics handles these interactions when compared to existing semantics. Fig. 1 gives an overview of which existing semantics can handle which (interactions of) exceptions. We note that our semantics could easily distinguish more kinds of exceptions, such as division by zero or out of bounds accesses to arrays.

*Non-termination and Observation Failure.* Listing 7 shows a program that has been investigated in [22]. Based on the observations, it only admits a single behavior, namely always sampling x = 0 in the third line. This behavior results in non-termination, but it occurs with probability 0. Hence, the program fails an observation (ending up in state -) with probability 1. If we try to

x:=0; **while** x=0 { x=**flip**( 1 2 ); **observe**(x=0); }

Listing 7. Mixing loops and observations condition on not failing any observation (by rescaling appropriately), this results in a division by 0, because the probability of not failing any observation is 0.

The semantics of Listing 7 thus only has weight on -, and does not allow conditioning on not failing any observation. This is also the solution that [22] proposes, but in our case, we can formally back up this claim with our semantics.

Other languages handle both non-termination and observation failure by subprobability distributions, which makes it impossible to conclude that the missing weight is due to observation failure (and not due to non-termination) [4,24,34]. The semantics in [28] cannot directly express that the missing weight is due to observation failure (rather, the semantics are undefined due to a division by zero). However, the semantics enables a careful reader to determine that the missing weight is due to observation failure (by investigating the conditional weakest precondition and the conditional weakest liberal precondition). Some other languages can express neither while loops nor observations [23,33,35].

*Assertions and Non-termination.* For some programs, it is useful to check assumptions explicitly. For example, the implementation of the factorial function in Listing 8 explicitly checks whether x is a valid argument to the factorial function. If x /∈ **N**, the program should run into an error (i.e. only have weight on ⊥). If x ∈ **N**, the program should return x! (i.e. only have weight on x!). This example illustrates that earlier exceptions (like failing an assertion) should *bypass* later exceptions (like nontermination, which occurs for x /∈ **N** if the programmer

```
assert(x≥0);
assert(x=x);
fac:=1;
while x=0 {
  fac=fac*x;
  x=x-1;
}
return fac;
```
Listing 8. Explicitly checking assumptions

forgets the first two assertions). This is not surprising, given that this is also the semantics of exceptions in most deterministic languages. Most existing semantics either cannot express Listing 8 ([23,34] have no assertions, [35] has no iteration) or cannot distinguish failing an assertion from non-termination [24,28,33]. The consequence of the latter is that removing the first two assertions from Listing 8 does not affect the semantics. Handling assertion failure by sum types (as e.g. in [34]) could be a solution, but would force the programmer to deal with assertion failure explicitly. Only the semantics in [4] has the expressiveness to implicitly handle assertion errors in Listing 8 without conflating those errors with non-termination.

Listing 9 shows a different interaction between nontermination and failing assertions. Here, even though the loop condition is always true, the first iteration of the loop will run into an exception. Thus, Listing 9 results in ⊥ with probability 1. Again, this behavior should not be surprising given the behavior of deterministic languages. For

x:=0; **while** 1 { x=x/x; }

Listing 9. Guaranteed failure

Listing 9, conflating errors with non-termination means the program semantics cannot express that the missing weight is due to an error and not due to nontermination.

*Observation Failure and Assertion Failure.* In our PPL, earlier exceptions bypass later exceptions, as illustrated in Listing 8. However, because we are operating in a probabilistic language, exceptions can occur probabilistically. Listing 10 shows a program that may run into

**observe**(**flip**( 1 2 )); **assert**(**flip**( 1 2 ));

Listing 10. Observation or assertion failure

an observation failure, or into an assertion failure, or neither. If it runs into an observation failure (with probability 1 2 ), it bypasses the rest of the program, resulting in with probability 1 2 and in ⊥ with probability 1 4 . Conditioning on the absence of observation failures, the probability of ⊥ is 1 2 .

An important observation is that reordering the two statements of Listing 10 will result in a different behavior. This is the case, even though there is no obvious data-flow between the two statements. This is in sharp contrast to the semantics in [34], which guarantee (in the absence of exceptions) that only data flow is relevant and that expressions can be reordered. Our semantics illustrate that even if there is no explicit data-dependency, some seemingly obvious properties (like commutativity) may not hold in the presence of exceptions. Some languages either cannot express Listing 10 ([23,33] lack observations), cannot distinguish observation failure from assertion failure [24] or cannot handle exceptions implicitly [34,35].

*Summary.* In this section, we showed examples of probabilistic programs that exhibit non-termination, observation failures and errors. Then, we provided examples that show how these exceptions can interact, and explained how existing semantics handle these interactions.

### 3 Preliminaries

In this section, we provide the necessary theory. Most of the material is standard, however, our treatment of exception states is interesting and important for providing semantics to probabilistic programs in the presence of exceptions. All key lemmas (together with additional definitions and examples) are proven in Appendix A.

*Natural Numbers,* [n]*, Iverson Brackets, Restriction of Functions.* We include 0 in the natural numbers, so that **N** := {0, 1,... }. For n ∈ **N**, [n] := {1,...,n}. The *Iverson brackets* [·] are defined by [b]=1 if b is true and [b]=0 if b is false. A particular application of the Iverson brackets is to characterize the indicator function of a specific set S by [x ∈ S]. For a function f : X → Y and a subset of the domain S ⊆ X, f restricted to S is denoted by f|S : S → Y .

*Set of Variables, Generating Tuples, Preservation of Properties, Singleton Set.* Let Vars be a set of admissible variable names. We refer to the elements of Vars by x, y, z and xi, yi, zi, vi, wi, for i ∈ **N**. For v ∈ A and n ∈ **N**, v!n := (v,..., v) ∈ An denotes the tuple containing n copies of v. A function f : An → A *preserves a property* if whenever a1,...,an ∈ A have that property, f(a1,...,an) ∈ A has that property. Let **1** denote the set which only contains the empty tuple (), i.e. **1** := {()}. For sets of tuples S ⊆ n i=1 Ai, there is an isomorphism S × **1 1** × S S. This isomorphism is intuitive and we sometimes silently apply it.

*Exception States, Lifting Functions to Exception States.* We allow the extension of sets with some symbols that stand for the occurrence of special events in a program. This is important because it allows us to capture the event that a given program runs into specific exceptions. Let X := {⊥, -, -} be a (countable) set of exception states.  Intuitively, ⊥ corresponds to assertion failures, corresponds to observation failures and corresponds to non-termination. For a function f : A → B, *f lifted to exception states*, denoted by f : A → B is defined by f(a) = a if a ∈ X and f(a) = f(a) if a /∈ X . For a function f : n i=1 Ai → B, *f lifted to exception states*, denoted by f : n i=1 Ai → B, propagates the first exception in its arguments, or evaluates f if none of its arguments are exceptions. Formally, it is defined by f(a1,...,an) = a1 if a1 ∈ X , f(a1,...,an) = a2 if a1 ∈ X / and a2 ∈ X , and so on. Only if a1,...,an ∈ X / , we have f(a1,...,an) = f(a1,...,an). Thus, f(-, a, ⊥) =-. In particular, we write (a, b) for lifting the tupling function, resulting in for example (-, -) = -. To remove notation clutter, we do not distinguish the two different liftings f : A → B and f : n i=1 Ai → B notationally. Whenever we write f, it will be clear from the context which lifting we mean. We write S×T for {(s, t) | s ∈ S, t ∈ T}.

*Records.* A *record* is a special type of tuple indexed by variable names. For sets (Si)i∈[n], a record r ∈ n i=1(xi : Si) has the form r = {x1 → v1,...,xn → vn}, where vi ∈ Si, with the convenient shorthand r = {xi → vi}i∈[n]. We can access the elements of a record by their name: r[xi] = vi.

In what follows, we provide the measure theoretic background necessary to express our semantics.

σ*-algebra, Measurable Set,* σ*-algebra Generated by a Set, Measurable Space, Measurable Functions.* Let A be some set. A set ΣA ⊆ P (A) is called a σ*-algebra on A* if it satisfies three conditions: A ∈ ΣA, ΣA is closed under complements (S ∈ ΣA implies A\S ∈ ΣA) and ΣA is closed under countable unions (for any collection {Si}i∈**N** with Si ∈ ΣA, we have i∈**N** Si ∈ ΣA). The elements of ΣA are called *measurable sets*. For any set A, a trivial σ-algebra on A is its power set P (A). Unfortunately, the power set often contains sets that do not behave well. To come up with a σ-algebra on A whose sets do behave well, we often start with a set S ⊆ P (A) that is not a σ-algebra and extend it until we get a σ-algebra. For this purpose, let A be some set and S ⊆ P (A) a collection of subsets of A. The σ*-algebra generated by S* denoted by σ(S) is the smallest σ-algebra that contains S. Formally, σ(S) is the intersection of all σ-algebras on A containing S. For a set A and a σ-algebra ΣA on A, (A, ΣA) is called a *measurable space*. We often leave ΣA implicit; whenever it is not mentioned explicitly, it is clear from the context. Table 1 provides the implicit σ-algebras for some common sets. As an example, some elements of Σ**R** include [0, 1] ∪ {⊥} and {1, 3, π}. For measurable spaces (A, ΣA) and (B,ΣB), a function f : A → B is called *measurable*,

Table 1. Implicit σ-algebras on common sets, for measurable spaces (A, Σ*A*), (A*i*, Σ*A*i )

if ∀S ∈ ΣB : f −1(S) ∈ ΣA. Here, f −1(S) := {a ∈ A: f(a) ∈ S}. If one is familiar with the notion of Lebesgue measurable functions, note that our definition does not include all Lebesgue measurable functions. As a motivation to why we need measurable functions, consider the following scenario. We know the distribution of some variable x, and want to know the distribution of y = f(x). To figure out how likely it is that y ∈ S for a measurable set S, we can determine how likely it is that x ∈ f −1(S), because f −1(S) is guaranteed to be a measurable set.

*Measures, Examples of Measures.* For a measurable space (A, ΣA), a function μ: ΣA → [0,∞] is called a *measure on A* if it satisfies two properties: null empty set (μ(∅)=0) and countable additivity (for any countable collection {Si}i∈I of pairwise disjoint sets Si ∈ ΣA, we have μ i∈I Si = - i∈I μ(Si)). Measures allow us to quantify the probability that a certain result lies in a measurable set. For example, μ([1, 2]) can be interpreted as the probability that the outcome of a process is between 1 and 2.

The *Lebesgue measure* λ: B → [0,∞] is the (unique) measure that satisfies λ([a, b]) = b − a for all a, b ∈ **R** with a ≤ b. The *zero measure* **0**: ΣA → [0,∞] is defined by **0**(S)=0 for all S ∈ ΣA. For a measurable space (A, ΣA) and some a ∈ A, the *Dirac measure* δa : ΣA → [0,∞] is defined by δa(S)=[a ∈ S].

Unfortunately, there are measures that do not satisfy some important properties (for example, they may not satisfy Fubini's theorem, which we discuss later on). The usual way to deal with this is to restrict our attention to σ-finite measures, which are well-known and were studied in great detail. However, σ-finite measures are too restrictive for our purposes. In particular, the s-finite kernels that we introduce later on can induce measures that are not σ-finite. This is why in the following, we work with s-finite measures. Table 2 gives an overview of the different kinds of measures that are important for understanding our work. The expression 1/2 · δ1 stands for the pointwise multiplication of the measure δ1 by 1/2: 1/2 · δ1 = λS. 1/2 · δ1(S). Here, the λ refers to λ-abstraction and not to the Lebesgue measure. To distinguish the two λs, we always write "λx." (with a dot) when we refer to λ-abstraction. For more details on the definitions and for proofs about the provided examples, see Appendix A.1.

Table 2. Definition and comparison of different measures μ: Σ*A* → [0, ∞] on measurable spaces (A, Σ*A*). Reading the table top-down, we get from the most restrictive definition to the most permissive definition. For example, any sub-probability measure is also a σ-finite measure. We also provide an example for each type of measure that is not an example of the more restrictive type of measure. For example, the Lebesgue measure λ is σ-finite but not s-finite.

*Product of Measures, Product of Measures in the Presence of Exception States.* For s-finite measures μ: ΣA → [0,∞] and μ : ΣB → [0,∞], we denote the *product of measures* by μ × μ : ΣA×B → [0,∞], and define it by

$$(\mu \times \mu')(S) = \int\_{a \in A} \int\_{b \in B} [(a, b) \in S] \mu'(db) \mu(da)$$

For s-finite measures μ: ΣA → [0,∞] and μ : ΣB → [0,∞], we denote the *lifted product of measures* by μ×μ : ΣA×B → [0,∞] and define it using the lifted tupling function: (μ×μ )(S) = a∈A b∈B[(a, b) ∈ S]μ (db)μ(da). While the product of measures μ × μ is well known for combining two measures to a joint measure, the concept of a lifted product of measures μ×μ is required to do the same for combining measures that have weight on exception states. Because the formal semantics of our probabilistic programming language makes use of exception states, we always use × to combine measures, appropriately handling exception states implicitly.

Lemma 1. *For measures* μ: ΣA → [0,∞]*,* μ : ΣB → [0,∞]*, let* S ∈ ΣA *and* T ∈ ΣB*. Then,* (μ × μ )(S × T) = μ(S) · μ (T)*.*

For μ: ΣA → [0,∞], μ : ΣB → [0,∞] and S ∈ ΣA, T ∈ ΣB, in general we have (μ×μ )(S × T) = μ(S) · μ (T), due to interactions of exception states.

Lemma 2. × *and* × *for s-finite measures are associative, left- and right-distributive and preserve (sub-)probability and s-finite measures.*

*Lebesgue Integrals, Fubini's Theorem for s-finite Measures.* Our definition of the Lebesgue integral is based on [31]. It allows integrating functions that sometimes evaluate to ∞, and Lebesgue integrals evaluating to ∞.

Here, (A, ΣA) and (B,ΣB) are measurable spaces and μ: ΣA → [0,∞] and μ : ΣB → [0,∞] are measures on A and B, respectively. Also, E ∈ ΣA and F ∈ ΣB. Let s: A → [0,∞) be a measurable function. s is a *simple function* if s(x) = n i=1 αi[x ∈ Ai] for Ai ∈ ΣA and αi ∈ **R**. For any simple function s, the Lebesgue integral of s over E with respect to μ, denoted by a∈E s(a)μ(da), is defined by n i=1 αi · μ(Ai ∩ E), making use of the convention 0 · ∞ = 0. Let f : A → [0,∞] be measurable but not necessarily simple. Then, the *Lebesgue integral* of f over E with respect to μ is defined by

$$\int\_{a \in E} f(a)\mu(da) := \sup \left\{ \int\_{a \in E} s(a)\mu(da) \, \middle| \, s: A \to [0, \infty) \text{ is simple}, 0 \le s \le f \right\}$$

Here, the inequalities on functions are pointwise. Appendix A.2 lists some useful properties of the Lebesgue integral. Here, we only mention Fubini's theorem, which is important because it entails a commutativity-like property of the product of measures: (μ × μ )(S)=(μ × μ)(swap(S)), where swap switches the dimensions of S: swap(S) = {(b, a) | (a, b) ∈ S}. The proof of this property is straightforward, by expanding the definition of the product of measures and applying Fubini's theorem. As we show in Sect. 5, this property is crucial for the commutativity of expressions. In the presence of exceptions, it does not hold: (μ×μ )(S) = (μ ×μ)(swap(S)) in general.

Theorem 1 (Fubini's theorem). *For s-finite measures* μ: ΣA → [0,∞] *and* μ : ΣB → [0,∞] *and any measurable function* f : A × B → [0,∞]*,*

$$\int\_{a \in A} \int\_{b \in B} f(a, b) \mu'(db) \mu(da) = \int\_{b \in B} \int\_{a \in A} f(a, b) \mu(da) \mu'(db)$$

*For s-finite measures* μ: ΣA → [0,∞] *and* μ : ΣB → [0,∞] *and any measurable function* f : A × B → [0,∞]*,*

$$\int\_{a \in \overline{A}} \int\_{b \in \overline{B}} \overline{f}(a, b) \mu'(db) \mu(da) = \int\_{b \in \overline{B}} \int\_{a \in \overline{A}} \overline{f}(a, b) \mu(da) \mu'(db)$$

*(Sub-)probability Kernels, s-finite Kernels, Dirac Delta, Lebesgue Kernel, Motivation for s-finite Kernels.* In the following, let (A, ΣA) and (B,ΣB) be measurable spaces. A *(sub-)probability kernel with source A and target B* is a function κ: A × ΣB → [0,∞) such that for all a ∈ A: κ(a, ·): ΣB → [0,∞) is a (sub-)probability measure, and ∀S ∈ ΣB : κ(·, S): A → [0,∞) is measurable. κ: A × ΣB → [0,∞] is an *s-finite kernel with source A and target B* if κ is a pointwise sum of sub-probability kernels κi : A × ΣB → [0,∞), meaning κ = - i∈**N** κi. We denote the set of s-finite kernels with source A and target B by A → B ⊆ A × ΣB → [0,∞]. Because we only ever deal with s-finite kernels, we often refer to them simply as kernels.

We can understand the Dirac measure as a probability kernel. For a measurable space (A, ΣA), the *Dirac delta* δ : A → A is defined by δ(a, S)=[a ∈ S]. Note that for any a, δ(a, ·): ΣA → [0,∞] is the Dirac measure. We often write δ(a)(S) or δa(S) for δ(a, S). Note that we can also interpret δ : A → A as an s-finite kernel from A → B for A ⊆ B. The *Lebesgue kernel* λ∗ : A → **R** is defined by λ∗(a)(S) = λ(S), where λ is the Lebesgue measure. The definition of s-finite kernels is a lifting of the notion of s-finite measures. Note that for an s-finite kernel κ, κ(a, ·) is an s-finite measure for all a ∈ A. In the context of probabilistic programming, s-finite kernels have been used before [34].

Working in the space of sub-probability kernels is inconvenient, because, for example, λ∗ : **R** → **R** is not a sub-probability kernel. Even though λ∗(x) is σfinite measure for all x ∈ **R**, not all s-finite kernels induce σ-finite measures in this sense. As an example, (λ∗;λ∗)(x) is not a σ-finite measure for any x ∈ **R** (see Lemma 15 in Appendix A.1). We introduce (;) shortly in Definition 1.

Working in the space of s-finite kernels is convenient because s-finite kernels have many nice properties. In particular, the set of s-finite kernels A → B is the smallest set that contains all sub-probability kernels with source A and target B and is closed under countable sums.

*Lifting Kernels to Exception States, Removing Weight from Exception States.* For kernels κ: A → B or kernels κ: A → B, κ *lifted to exception states* κ: A → B is defined by κ(a) = κ(a) if a ∈ A and κ(a) = δ(a) if a /∈ A. When transforming κ into κ, we preserve (sub-)probability and s-finite kernels.

*Composing kernels, composing kernels in the presence of exception states.*

Definition 1. *Let* (;): (A → B) → (B → C) → (A → C) *be defined by* (f;g)(a)(S) = b∈B g(b)(S) f(a)(db)*.*

Note that f;g intuitively corresponds to first applying f and then g. Throughout this paper, we mostly use >=> instead of (;), but we introduce (;) because it is well-known and it is instructive to show how our definition of >=> relates to (;).

Lemma 3. (;) *is associative, left- and right-distributive, has neutral element*2 δ *and preserves (sub-)probability and s-finite kernels.*

Definition 2. *Let* (>=>): (A → B) → (B → C) → (A → C) *be defined by* (f >=> g)(a)(S) = b∈B g(b)(S) f(a)(db)*.*

We sometimes write f(a) = g for (f >=> g)(a).

Lemma 4. *For* f : A → B *and* g : B → C*,* a ∈ A *and* S ∈ ΣC *,*

$$(f \succ \succ g)(a)(S) = (f;g)(a)(S) + \sum\_{x \in \mathcal{X}} \delta(x)(S)f(a)(\{x\})$$

Lemma 4 shows how >=> relates to (;), by splitting f >=> g into nonexceptional behavior of f (handled by (;)) and exceptional behavior of f (handled by a sum). Intuitively, if f produces an exception state ∈ X , then g is not even evaluated. Instead, this exception is directly passed on, as indicated by δ(x)(S).

2 δ is a neutral element of (;) if (δ;κ)=(κ;δ) = κ for all kernels κ.

If f(a)(X )=0 for all a ∈ A, or if S ∩ X = ∅, then the definitions are equivalent in the sense that (f;g)(a)(S)=(f >=> g)(a)(S). The difference between >=> and (;) is the treatment of exception states produced by f. Note that technically, the target B of f : A → B does not match the source B of g : B → C. Therefore, to formally interpret f;g, we silently restrict the domain of f to A × ΣB.

Lemma 5. >=> *is associative, left-distributive (but not right-distributive), has neutral element* δ *and preserves (sub-)probability and s-finite kernels.*

*Product of Kernels, Product of Kernels in the Presence of Exception States.* For s-finite kernels κ: A → B, κ : A → C, we define the *product of kernels*, denoted by κ × κ : A → B × C, as (κ × κ )(a)(S)=(κ(a) × κ (a))(S). For s-finite kernels κ: A → B and κ : A → C, we define the *lifted product of kernels*, denoted by κ×κ : A → B × C, as (κ×κ )(a)(S)=(κ(a)×κ (a))(S). × and × allow us to combine kernels to a joint kernel. Essentially, this definition reduces the product of kernels to the product of measures.

Lemma 6. × *and* × *for kernels preserve (sub-)probability and s-finite kernels, are associative, left- and right-distributive.*

*Binding Conventions.* To avoid too many parentheses, we make use of some binding conventions, ordering (in decreasing binding strength) ×, ×, ;, >=>, +.

*Summary.* The most important concepts introduced in this section are exception states, records, Lebesgue integration, Fubini's theorem and (s-finite) kernels.

### 4 A Probabilistic Language and Its Semantics

We now describe our probabilistic programming language, the typing rules and the denotational semantics of our language.

### 4.1 Syntax

Let **V** := **Q** ∪ {π, e} ⊆ **R** be a (countable) set of constants expressible in our programs. Let i, n ∈ **N**, r ∈ **V**, x ∈ Vars, a generic unary operator (e.g., − inverts the sign of a value, ! is logical negation mapping 0 to 1 and all other numbers to 0, · and · round down and up respectively), ⊕ a generic binary operator (e.g., +, −, ∗, /, ∧ for addition, subtraction, multiplication, division and exponentiation, &&, || for logical conjunction and disjunction, =, =, <, ≤, >, ≥ to compare values). Let f : A → **R** → [0,∞) be a measurable function that maps a ∈ A to a probability density function. We check if f is measurable by uncurrying f to f : A × **R** → [0,∞). Fig. 2 shows the syntax of our language.

Our expressions capture () (the only element of **1**), r (real numbers), x (variables), (e1,...,en) (tuples), e[i] (accessing elements of tuples for i ∈ **N**), e (unary operators), e1 ⊕ e2 (binary operators), e1[e2] (accessing array elements), e1[e2 → e3] (updating array elements), **array**(e1, e2) (creating array of length e1

$$e \implies ()\mid r \mid x \mid (e\_1, \ldots, e\_n) \mid e[i] \mid \in e \mid e\_1 \oplus e\_2 \mid e\_1[e\_2] \mid \tag{\textbf{Expressions}}$$

$$\mathtt{assert}(e\_1, e\_2) \mid e\_1[e\_2 \mapsto e\_3] \mid F(e)$$

$$F \implies \lambda x.\{P; \mathtt{return} \; e;\} \mid \mathtt{`\mathtt{List} \; \mathtt{//} \mathtt{while} \, \mathtt{m} \, \mathtt{//} \mathtt{same} \, \mathtt{IoFron} \} \tag{\textbf{Functions}}$$

$$P \implies \mathtt{skip} \mid x := e \mid x = e \mid P\_1; P\_2 \mid \mathtt{if} \; e \; \{P\_1\} \; \mathtt{else} \; \{P\_2\} \mid \{P\} \mid \quad \text{(\textbf{Statements})}$$

$$\mathtt{assert}(e) \mid \mathtt{observe}(e) \mid \mathtt{while} \; e \; \{P\}$$

Fig. 2. The syntax of our probabilistic language.

containing e2 at every index) and F(e) (evaluating function F on argument e). To handle functions F(e1,...,en) with multiple arguments, we interpret (e1,...,en) as a tuple and apply F to that tuple.

Our functions express λx.{P; **return** e; } (function taking argument x running P on x and returning e), **flip**(e) (random choice from {0, 1}, 1 with probability e), **uniform**(e1, e2) (continuous uniform distribution between e1 and e2) and **sampleFrom**f (e) (sample value distributed according to probability density function f(e)). An example for f is the density of the exponential distribution, indexed with rate λ. Formally, f : (0,∞) → **R** → [0,∞) is defined by f(λ)(x) = λe−λx if x ≥ 0 and f(λ)(x)=0 otherwise. Often, f is partial (e.g., λ ≤ 0 is not allowed). Intuitively, arguments outside the allowed range of f produce the error state ⊥.

Our statements express **skip** (no operation), x := e (assigning to a fresh variable), x = e (assigning to an existing variable), P1; P2 (sequential composition of programs), **if** e {P1} **else** {P2} (if-then-else), {P} (static scoping), **assert**(e) (asserting that an expression evaluates to true, assertion failure results in ⊥), **observe**(e) (observing that an expression evaluates to true, observation failure results in -) and **while** e {P} (while loops, non-termination results in -). We additionally introduce syntactic sugar e1[e2] = e3 for e1 = e1[e2 → e3], **if** (e) {P} for **if** e {P} **else** {**skip**} and func(e2) for λx.{P; **return** e1; }(e2) (using the name func for the function with argument x and body {P; **return** e1}).

### 4.2 Typing Judgments

Let n ∈ **N**. We define types by the following grammar in BNF, where τ [] denotes arrays over type τ . We sometimes write n i=1 τi for the product type τ1×···×τn.

$$\tau ::= \mathbb{1} \mid \mathbb{R} \mid \tau \lbrack \ | \ \tau\_1 \times \dots \times \tau\_n \rbrack$$

Note that we also use the type τ1 → τ2 of kernels with source τ1 and target τ2, but we do not list it here to avoid higher-order functions (discussed in Sect. 4.5).

Formally, a *context* Γ is a set {xi : τi}i∈[n] that assigns a type τi to each variable xi ∈ Vars. In slight abuse of notation, we sometimes write x ∈ Γ if there is a type τ with x: τ ∈ Γ. We also write Γ, x: τ for Γ ∪ {x: τ} (where x /∈ Γ) and Γ, Γ for Γ ∪ Γ (where Γ and Γ have no common variables).

Fig. 3. The typing rules for expressions and functions in our language

Fig. 4. The typing rules for statements

The rules in Figs. 3 and 4 allow deriving the type of expressions, functions and statements. To state that an expression e is of type τ under a context Γ, we write Γ e: τ . Likewise, F : τ → τ indicates that F is a kernel from τ to τ . Finally, Γ P Γ states that a context Γ is transformed to Γ by a statement P. For **sampleFrom**f , we intuitively want f to map values from τ to probability density functions. To allow f to be partial, i.e., to be undefined for some values from τ , we use A ∈ Στ (and hence A ⊆ [[τ ]]) as the domain of f (see Sect. 4.3).

#### 4.3 Semantics

*Semantic Domains.* We assign to each type τ a set [[τ ]] together with an implicit σ-algebra Στ on that set. Additionally, we assign a set [[Γ]] to each context Γ = {xi : τi}i∈[n]. Concretely, we have [[**1**]] = **1** := {()} with Σ**1** = {∅,()}, [[**R**]] = **R** and Σ**R** = B. The remaining semantic domains are outlined in Fig. 5.

$$\begin{aligned} \left[\tau[\ ]\right] &= \bigcup\_{i \in \mathbb{N}} \left[\tau\ \right]^i & \Sigma\_{\tau[\ ]} \text{ is generated by } \bigcup\_{i \in \mathbb{N}} \left\{ \prod\_{j=1}^i S\_j \; \middle| \; S\_j \in \Sigma\_{\tau} \right\} \\ \left[\tau\_1 \ \times \cdots \times \tau\_n\right] &= \prod\_{i=1}^n \left[\tau\_i\right] & \Sigma\_{\tau\_1 \times \cdots \times \tau\_n} \text{ is generated by } \left\{ \prod\_{i=1}^n S\_i \; \middle| \; S\_i \in \Sigma\_{\tau\_i} \right\} \\ \left[\Gamma\right] &= \prod\_{i=1}^n \left(x\_i \colon \left\|\tau\_i\right\|\right) & \Sigma\_{\Gamma} \text{ is generated by } \left\{ \prod\_{i=1}^n \left(x\_i \colon S\_i\right) \; \middle| \; S\_i \in \Sigma\_{\tau\_i} \right\} \end{aligned}$$

### Fig. 5. Semantic domains for types

Fig. 6. The semantics of expressions. v!n stands for the n-tuple (v,..., v). t[i] stands for the i-th element (0-indexed) of the tuple t and t[i → v] is the tuple t, where the i-th element is replaced by v. |t| is the length of a tuple t. σ stands for a program state over all variables in some Γ, with σ ∈ [[Γ]].

*Expressions.* Fig. 6 assigns to each expression e typed by Γ e: τ a probability kernel [[e]]τ : [[Γ]] → [[τ ]]. When τ is irrelevant or clear from the context, we may drop it and write [[e]]. The formal interpretation of [[Γ]] → [[τ ]] is explained in Sect. 3. 3 Note that Fig. 6 is incomplete, but extending it is straightforward. When we need to evaluate multiple terms (as in (e1,...,en)), we combine the results using ×. This makes sure that in the presence of exceptions, the first exception that occurs will have priority over later exceptions. In addition, deterministic functions (like x+y) are lifted to probabilistic functions by the Dirac delta (e.g. δ(x+y)) and incomplete functions (like x/y) are lifted to complete functions via the explicit error state ⊥.

3 As a quick and intuitive reminder, κ: A → B means that for every a ∈ A, κ(a) will be a distribution over B, where B is B enriched with exception states. Hence, κ(a) may have weight on elements of B, on exception states, or on both.

Fig. 7 assigns to each function F typed by F : τ1 → τ2 a probability kernel [[F]]τ1 →τ2 : [[τ1]] → [[τ2]]. In the semantics of **flip**, δ(1): Σ**R** → [0,∞] is a measure on **R**, and p · δ(1) rescales this measure pointwise. Similarly, the sum p · δ(1) + (1 − p)· δ(0) is also meant pointwise, resulting in a measure on **R**. Finally, λp. p · δ(1)+ (1−p)·δ(0) is a kernel with source [0, 1] and target **R**. For **sampleFrom**f (e), remember that f(p)(·) is a probability density function.

$$\begin{aligned} [\mathsf{f1.ip}]\_{\mathbb{R}\rightarrow\mathbb{R}} &= \lambda p. \begin{cases} p \cdot \delta(1) + (1-p) \cdot \delta(0) & p \in [0,1] \\ \delta(\bot) & \text{otherwise} \end{cases} \\ [\mathsf{unifor}]\_{\mathbb{R}\rightarrow\mathbb{R}} &= \lambda(l,r). \begin{cases} \lambda S. \frac{1}{r-l} \lambda([l,r] \cap S) & l < r \\ \delta(\bot) & \text{otherwise} \end{cases} \\ [\mathsf{sameloFrom}\_{f}]\_{\mathbb{T}\rightarrow\mathbb{R}} &= \lambda p. \begin{cases} \lambda S. \int\_{x\in\mathbb{R}\cap S} f(p)(x) \lambda(dx) & p \in A \\ \delta(\bot) & p \notin A \end{cases} \\ [\lambda x.\{P; \mathsf{return} \ e\_{i}\}]\_{\mathbb{T}\_{1}\rightarrow\mathbb{T}\_{2}} &= \lambda v. \delta\left(\{x \mapsto v\}\right) \geqslant = \left\{ \left[P\right] \geqslant = \left[e\_{2}\right]\_{\mathbb{T}\_{2}} \right\} \end{aligned}$$

Fig. 7. The semantics of functions.

$$\begin{aligned} \left[\mathsf{skLip}\right] = \delta & \quad \left[x:=e\right] = \left[x=e\right] = \delta \overline{\times} \left[e\right] & \geqslant \lambda(\sigma, v). \delta(\sigma[x \mapsto v]) \\\\ \left[P\_1; P\_2\right] = \left[P\_1\right] & \geqslant \gg \left[P\_2\right] \quad \quad \quad \left[\left\{P\right\}\right] = \left[P\right] & \geqslant \geqslant \lambda \sigma'. \delta(\sigma'(I)) \\\\ \left[\mathsf{i}\ \mathsf{f}\ \ e\left\{P\_1\right\}\ \mathsf{else}\ \left\{P\_2\right\}\right] &= \delta \overline{\times} \left[e\right]\_{\mathbb{R}} \geqslant \geqslant \lambda(\sigma, b). \begin{cases} \left[\left[P\_1\right]\right](\sigma) & b \neq 0 \\ \left[\left[P\_2\right]\right](\sigma) & b = 0 \end{cases} \\\\ \left[\mathsf{asser}\,\mathsf{r}(e)\right] &= \delta \overline{\times} \left[e\right]\_{\mathbb{R}} \geqslant \geqslant \lambda(\sigma, b). \begin{cases} \delta(\sigma) & b \neq 0 \\ \delta(\bot) & b = 0 \end{cases} \\\\ \left[\mathsf{observe}(e)\right] &= \delta \overline{\times} \left[e\right]\_{\mathbb{R}} \geqslant \geqslant \lambda(\sigma, b). \begin{cases} \delta(\sigma) & b \neq 0 \\ \delta(\not\bot) & b = 0 \end{cases} \end{aligned}$$

Fig. 8. The semantics of programs in our probabilistic language. Here, σ[x → v] results in σ with the value stored under x updated to v. σ (Γ) selects only those variables from σ that occur in Γ, meaning {x*i* → v*i*}*i*∈I ({x*i* : τ*i*}*i*∈I- ) = {x*i* → v*i*}*i*∈I∩I-.

*Statements.* Fig. 8 assigns to each statement P with Γ P Γ a probability kernel [[P]]: [[Γ]] → [[Γ ]]. Note the use of × in δ×[[e]], which allows evaluating e while keeping the state σ in which e is being evaluated. Intuitively, if evaluating e results in an exception from X , the previous state σ is irrelevant, and the result of δ×[[e]] will be that exception from X .

*While Loop.* To define the semantics of the while loop **while** e {P}, we introduce a *kernel transformer* [[**while** e {P}]]trans : ([[Γ]] → [[Γ]]) → ([[Γ]] → [[Γ]]) that transforms the semantics for n runs of the loop to the semantics for n + 1 runs of the loop. Concretely,

$$\{\mathsf{while}\ \mathsf{e}\ \{P\}\}^{\mathsf{trans}}(\kappa) = \delta \overline{\times} [e] \succ \Rightarrow \lambda(\sigma, b). \begin{cases} \{P\}(\sigma) \ggg \kappa \, b \neq 0\\ \delta(\sigma) \qquad b = 0 \end{cases}$$

This semantics first evaluates e, while keeping the program state around using δ. If e evaluates to 0, the while loop terminates and we return the current program state σ. If e does not evaluate to 0, we run the loop body P and feed the result to the next iteration of the loop, using κ.

We can then define the semantics of **while** e {P} using a special fixed point operator fix: ((A → A) → (A → A)) → (A → A), defined by the pointwise limit fix(Δ) = limn→∞ Δn(-), where -:= λσ. δ(-) and Δn denotes the n-fold composition of Δ. Δn(-) puts all runs of the while loop that do not terminate within n steps into the state -. In the limit, only has weight on those runs of the loop that never terminate. fix(Δ) is only defined if its pointwise limit exists. Making use of fix, we can define the semantics of the while loop as follows:

$$\left[\mathsf{while}\,\mathsf{i}\,\mathsf{i}\,\mathsf{i}\,\left\{P\right\}\right] = \mathsf{fix}\left(\left[\mathsf{while}\,\mathsf{i}\,\mathsf{i}\,\mathsf{i}\,\left\{P\right\}\right]^{\mathsf{trans}}\right),$$

Lemma 7. *For* Δ *as in the semantics of the while loop, and for each* σ *and each* S*, the limit* limn→∞ Δn(-)(σ)(S) *exists.*

Lemma 7 holds because increasing n may only shift probability mass from to other states (we provide a formal proof in Appendix B). Kozen shows a different way of defining the semantics of the while loop [23], using least fixed points. Lemma 8 describes the relation of the semantics of our while loop to the semantics of the while loop of [23]. For more details on the formal interpretation of Lemma 8 and for its proof, see Appendix B.

Lemma 8. *In the absence of exception states, and using sub-probability kernels instead of distribution transformers, the definition of the semantics of the while loop from [23] is equivalent to ours.*

Theorem 2. *The semantics of each expression* [[e]] *and statement* [[P]] *is indeed a probability kernel.*

*Proof.* The proof proceeds by induction. Some lemmas that are crucial for the proof are listed in Appendix C. Conveniently, most functions that come up in our definition are continuous (like a+b) or continuous except on some countable subset (like a b ) and thus measurable.

### 4.4 Recursion

To extend our language with recursion, we apply the same ideas as for the while loop. Given the source code of a function F that uses recursion, we define its

$$\delta \overline{\times} \left[ ! \mathsf{tfu} \mathsf{i} \mathfrak{p} \left( \frac{1}{2} \right) \right] \rangle \rightleftharpoons \geqslant \lambda (\sigma, b) . \left\{ \begin{pmatrix} \kappa \overline{\times} [1] \\ \end{pmatrix} \rhd \simeq \lambda (x, y) . \delta \left( x + y \right) \right\} (\sigma) \quad b \neq 0$$

Fig. 9. Kernel transformer [[geom]]trans(κ) for geom given in Listing 11.

semantics in terms of a kernel transformer [[F]]trans. This kernel transformer takes semantics for F up to a recursion depth of n and returns semantics for F up to recursion depth n+ 1. Formally, [[F]]trans(κ) follows the usual semantics, but uses κ as the semantics for recursive calls to F (we will provide an example shortly). Finally, we define the semantics of F by [[F]] := fix [[F]]trans . Just as for the while loop, fix [[F]]trans is well-defined because stepping from recursion depth n to n+ 1 can only shift probability mass from to other states. We note that we could generalize our approach to mutual recursion.

To demonstrate how we define the kernel transformer, consider the recursive implementation of the geometric distribution in Listing 11 (to simplify presentation, Listing 11 uses early return). Given semantics κ for geom : **1** → **R** up to recursion depth n, we can define the semantics of geom up to recursion depth n + 1, as illustrated in Fig. 9.

### 4.5 Higher-Order Functions

Our language cannot express higher-order functions. When trying to give semantics to higher-order probabilistic programs, an important step is to define a σalgebra on the set of functions from real numbers to real numbers. Unfortunately, no matter which σ-algebra is picked, function evaluation (i.e. the function that takes f and x as arguments and returns f(x)) is not measurable [1]. This is a known limitation that previous work has looked into (e.g. [35] address it by restricting the set of functions to those expressible by their source code).

A promising recent approach is replacing measurable spaces by quasi-Borel spaces [16]. This allows expressing higher-order functions, at the price of replacing the well-known and well-understood measurable spaces by a new concept.

#### 4.6 Non-determinism

To extend our language with non-determinism, we may define the semantics of expressions, functions and statements in terms of sets of kernels. For an expression e typed by Γ e : τ , this means that [[e]]τ ∈ P ([[Γ]] → [[τ ]]), where P (S) denotes the power set of S. Lifting our semantics to non-determinism is mostly straightforward, except for loops. There, [[**while** e {P}]] contains all kernels of the form limn→∞(Δ1 ◦··· ◦ Δn)(-), where Δi ∈ [[**while** e {P}]]trans. Previous work has studied non-determinism in more detail, see e.g. [21,22].

### 5 Properties of Semantics

We now investigate two properties of our semantics: commutativity and associativity. These are useful in practice, e.g. because they enable rewriting programs to a form that allows for more efficient inference [5].

In this section, we write e1 e2 when expressions e1 and e2 are equivalent (i.e. when [[e1]] = [[e2]]). Analogously, we write P1 P2 for [[P1]] = [[P2]].

### 5.1 Commutativity

In the presence of exception states, our language cannot guarantee commutativity of expressions such as e1 + e2. This is not surprising, as in our semantics the first exception bypasses all later exceptions.

Lemma 9. *For function* F(){**while** 1 {**skip**}; **return** 0}*,*

$$\frac{1}{0} + F() \not\equiv F() + \frac{1}{0}$$

Formally, this is because if we evaluate 1 0 first, we only have weight on ⊥. If instead, we evaluate F() first, we only have weight on -, by an analogous calculation. A more detailed proof is included in Appendix D.

However, the only reason for non-commutativity is the presence of exceptions. Assuming that e1 and e2 cannot produce exceptions, we obtain commutativity:

Lemma 10. *If* [[e1]](σ)(X ) = [[e2]](σ)(X )=0 *for all* σ*, then* e1 ⊕ e2 e2 ⊕ e1*, for any commutative operator* ⊕*.*

The proof of Lemma 10 (provided in Appendix D) relies on the absence of exceptions and Fubini's Theorem. This commutativity result is in line with the results from [34], which proves commutativity in the absence of exceptions.

In the analogous situation for statements, we cannot assume commutativity P1; P2 P2; P1, even if there is no dataflow from P1 to P2. We already illustrated this in Listing 10, where swapping two lines changes the program semantics. However, in the absence of exceptions and dataflow from P1 to P2, we can guarantee P1; P2 P2; P1.

#### 5.2 Associativity

A careful reader might suspect that since commutativity does not always hold in the presence of exceptions, a similar situation might arise for associativity of some expressions. As an example, can we guarantee e1+(e2+e3) (e1+e2)+e3, even in the presence of exceptions? The answer is yes, intuitively because exceptions can only change the behavior of a program if the order of their occurrence is changed. This is not the case for associativity. Formally, we derive the following:

Lemma 11. e1 ⊕ (e2 ⊕ e3) (e1 ⊕ e2) ⊕ e3*, for any associative operator* ⊕*.*

We include notes on the proof of Lemma 11 in Appendix D, mainly relying on the associativity of × (Lemma 6). Likewise, sequential composition is associative: (P1; P2); P3 P1; (P2; P3). This is due to the associativity of >=> (Lemma 5).

### 5.3 Adding the **score** Primitive

Some languages include the primitive **score**, which allows to increase or decrease the probability of a certain event (or trace) [34,35].

Listing 12 shows an example program using **score**. Without normalization, it returns 0 with probability 1 2 and 1 with "probability" 1 2 · 2=1. After normalization, it returns 0 with probability 1 3 and 1 with probability 2 3 . Because **score** allows decreasing the probability of a specific event, it renders **observe** unnecessary. In general, we can replace **observe**(e) by **score**(e = 0). However, performing this replacement means losing the explicit knowledge of the weight on -

**score** can be useful to modify the shape of a given distribution. For example, Listing 13 turns the distribution of x, which is a Gaussian distribution, into the Lebesgue measure λ, by multiplying the density of x by its inverse. Hence, the density of x at any location is 1. Note that the distribution over x cannot be described by a probability measure, because e.g. the "probability" that x lies in the interval [0, 2] is 2.

Unfortunately, termination in the presence of **score** is not well-defined, as illustrated in Listing 14. In this program, the only non-terminating trace keeps changing its weight, switching between 1 and 2. In the limit, it is impossible to determine the weight of non-termination.

Hence, allowing the use of the **score** primitive only makes sense after abolishing the tracking of nontermination (-), which can be achieved by only measuring sets that do not contain non-termination. Formally, this means restricting the semantics of expressions e typed by Γ e : τ to [[e]]τ : Γ → [[τ ]] − {-} . x:=**flip**( 1 2 ); **if** x=1 { **score**(2); } **return** x; Listing 12. Using **score**

. x:=**gauss**(0,1); **score**( √2πe*x*2*/*2); **return** x;

Listing 13. Reshaping a distribution.

i:=0; **while** 1 {

```
if i=0 {
    score(2);
}else{
    score( 1
            2 );
}
i=1-i;
```
} Listing 14. **score** vs non-termination

Intuitively, abolishing non-termination means that we ignore non-terminating runs (these result in weight on non-termination). After doing this, we can give well-defined semantics to the **score** primitive.

The typing rule and semantics of **score** are:

$$\frac{\Gamma \vdash e : \mathbb{R}}{\Gamma \stackrel{\mathsf{sear}(e)}{\mathsf{succ}} \Gamma} \qquad \text{and} \quad \mathsf{[scorre(e)]} = \delta \overline{\times} [e]\_{\mathbb{R}} \iff \lambda(\sigma, c).c \* \delta(\sigma)$$

After including **score** into our language, the semantics of the language can no longer be expressed in terms of probability kernels as stated in Theorem 2, because the probability of any event can be inflated beyond 1. Instead, the semantics must be expressed in terms of s-finite kernels.

Theorem 3. *After adding the score primitive and abolishing non-termination, the semantics of each expression* [[e]] *and statement* [[P]] *is an s-finite kernel.*

Table 3. Comparison of existing semantics to ours. When adding **score** to our language (Sect. 5.3), our semantics use s-finite kernels (not probability kernels).

*Proof.* As for Theorem 2, the proof proceeds by induction. Most parts of the proof are analogous (e.g. >=> preserves s-finite kernels instead of probability kernels). For while loops, the limit still exists (Lemma 7 still holds), but it is not bounded from above anymore. The limit indeed corresponds to an s-finite kernel because the limit of strictly increasing s-finite kernels is an s-finite kernel.

**score**(2); **assert** In the presence of **score**, we can still talk about the interaction of different exceptions, assuming that we do track different types of exceptions (e.g. division by zero and out of bounds access of arrays). Then, we keep the commutativity and associativity properties studied in the previous sections, because these still hold for s-finite kernels.

Listing 15 shows an interaction of **score** with **assert**. As one would expect, our semantics will assign weight 2 to ⊥ in this case. If the two statements are switched, our semantics will ignore **score**(2) and assign weight 1 to ⊥. Hence again, commutativity does not hold.

```
assert(false);
Listing 15. Inter-
```
action of **score** and

```
while 1 {
    score(2);
    assert(flip( 1
                  2 ));
}
Listing 16. Interaction of
score, assert and loops
```
Listing 16 shows a program that keeps increasing the probability of an error. In every loop iteration, there is a "probability" of 1 of running into an error. Overall, Listing 16 results in weight ∞ on state ⊥.

### 6 Related Work

Kozen provides classic semantics to probabilistic programs [23]. We follow his main ideas, but deviate in some aspects in order to introduce additional features or to make our presentation cleaner. The semantics by Hur et al. [19] is heavily based on [23], so we do not go into more detail here. Table 3 summarizes the comparison of our approach to that of others.

*Kernels.* Like our work, most modern approaches use kernels (i.e., functions from values to distributions) to provide semantics to probabilistic programs [4,24,33, 34]. Borgström et al. [4] use sub-probability kernels on (symbolic) expressions. Staton [34] uses s-finite kernels to capture the semantics of the **score** primitive (when we discuss **score** in Sect. 5.3, we do the same).

In the classic semantics of [23], Kozen uses distribution transformers (i.e., functions from distributions to distributions). For later work [24], Kozen also switches to sub-probability kernels, which has the advantage of avoiding redundancies. A different approach uses weakest precondition to define the semantics, as in [28]. Staton et al. [35] use a different concept of measurable functions A → P(**R**≥0 ×B) (where P(S) denotes the set of all probability measures on S).

*Typing.* Some probabilistic languages are untyped [4,28], while others are limited to just a single type: **R**n [23,24] or ∞ i=1 **N**i ∪ **N**∞ [33]. Some languages provide more interesting types including sum types, distribution types and tuples [34,35]. We allow tuples and array types, and we could easily account for sum types.

*Loops.* Because the semantics of while loops is not always straightforward, some languages avoid while loops and recursion altogether [35]. Borgström et al. handle recursion instead of while loops, defining the semantics in terms of a fixed point [4]. Many languages handle while loops by least fixed points [23,24,28,33]. Staton defines while loops in terms of the counting measure [34], which is similar to defining them by a fixed point. We define the semantics of while loops in terms of a fixed point, which avoids the need to prove the least fixed point exists (still, the classic while loop semantics of [23] and our formulation are equivalent).

Most languages do not explicitly track non-termination, but lose probability weight by non-termination [4,23,24,34]. This missing weight can be used to identify the probability of non-termination, but only if other exceptions (such as **fail** in [24] or observation failure in [4]) do not also result in missing weight. The semantics of [33] are tailored to applications in networks and lose *non-terminating packet histories* instead of weight (due to a particular least fixed point construction of Scott-continuous maps on algebraic and continuous directed complete partial orders). Some works define non-termination as missing weight in the weakest precondition [28]. Specifically, the semantics in [28] can also explicitly express probability of non-termination *or* ending up in some state (using the separate construct of a weakest liberal precondition). We model nontermination by an explicit state -, which has the advantage that in the context of lost weight, we know what part of that lost weight is due to non-termination.

Kaminski et al. [21] investigate the run-time of probabilistic program with loops and **fail** (interpreted as early termination), but without observations. In [21], non-termination corresponds to an infinite run-time.

*Error States.* Many languages do not consider partial functions (like fractions a b ) and thus never run into an exception state [23,24,33]. Olmedo et al. [28] do not consider partial functions, but support the related concept of an explicit **abort**. The semantics of **abort** relies on missing weight in the final distribution. Some languages handle expressions whose evaluation may fail using sum types [34,35], forcing the programmer to deal with errors explicitly (we discuss the disadvantages of this approach at Listing 6). Formally, a sum type A + B is a disjoint union of the two sets A and B. Defining the semantics of an expression in terms of the sum type A+{⊥} allows that expression to evaluate to *either* a value a ∈ A *or* to ⊥.  Our semantics also uses sum types to handle exceptions, but the handling is implicit, by defining semantics in terms of (>=>) (which defines how exceptions propagate in a program) instead of (;).

*Constraints.* To enforce hard constraints, we use the **observe**(e) statement, which puts the program into a special failure state if it does not satisfy e. We can encode soft constraints by **observe**(e), where e is probabilistic (this is a general technique). Borgström et al. [4] allow both soft constraints that reduce the probability of some program traces and hard constraints whose failure leads to the error state **fail**. Some languages can handle generalized soft constraints: they can not only decrease the probability of certain traces using soft constraints, but also increase them, using **score**(x) [34,35]. We investigate the consequences of adding **score** to our language in Sect. 5.3. Kozen [24] handles hard (and hence soft) constraints using **fail** (which results in a sub-probability distribution). Some languages can handle neither hard nor soft constraints [23,33]. Note though that the semantics of ProbNetKAT in [33] can drop certain packages, which is a similar behavior. Olmedo et al. [28] handle hard (and hence soft) constraints by a conditional weakest precondition that tracks both the probability of not failing any observation and the probability of ending in specific states. Unfortunately, this work is restricted to discrete distributions and is specifically designed to handle observation failures and non-termination. Thus, it is not obvious how to adapt the semantics if a different kind of exception is to be added.

*Interaction of Different Exceptions.* Most existing work handles at least some exceptions by sub-probability distributions [4,23,24,33,34]. Then, any missing weight in the final distribution must be due to exceptions. However, this leads to a conflation of all exceptions handled by sub-probability distributions (for the consequences of this, see, e.g., our discussion of Listing 8). Note that semantics based on sub-probability kernels can add more exceptions, but they will simply be conflated with all other exceptions.

Some previous work does not (exclusively) rely on sub-probability distributions. Borgström et al. [4] handle errors implicitly, but still use sub-probability kernels to handle non-termination and **score**. Olmedo et al. can distinguish nontermination (which is conflated with exception failure) from failing observations by introducing two separate semantic primitives (conditional weakest precondition and conditional liberal weakest precondition) [28]. Because their solution specifically addresses non-termination, it is non-trivial to generalize this treatment to more than two exception states. By using sum types, some semantics avoid interactions of errors with non-termination or constraint failures, but still cannot distinguish the latter [34,35]. Note that semantics based on sum types can easily add more exceptions (although it is impossible to add non-termination). However, the interaction of different exceptions cannot be observed, because the programmer has to handle exceptions explicitly.

To the best of our knowledge, we are the first to give formal semantics to programs that may produce exceptions in this generality. One work investigates assertions in probabilistic programs, but explicitly disallows non-terminating loops [32]. Moreover, the semantics in [32] are operational, leaving the distribution (in terms of measure theory) of program outputs unclear. Cho et al. [8] investigate the interaction of partial programs and observe, but are restricted to discrete distributions and to only two exception states. In addition, this investigation treats these two exception states differently, making it non-trivial to extend the results to three or more exception states. Katoen et al. [22] investigate the intuitive problems when combining non-termination and observations, but restrict their discussions to discrete distributions and do not provide formal semantics. Huang [17] treats partial functions, but not different kinds of exceptions. In general, we know of no probabilistic programming language that distinguishes more than two different kinds of exceptions. Distinguishing two kinds of exceptions is simpler than three, because it is possible to handle one exception as an explicit exception state and the other one by missing weight (as e.g. in [4]).

Cousot and Monerau [9] provide a trace semantics that captures probabilistic behavior by an explicit randomness source given to the program as an argument. This allows handling non-termination by non-terminating traces. While the work does not discuss errors or observation failure, it is possible to add both. However, using an explicit randomness source has other disadvantages, already discussed by Kozen [23]. Most notably, this approach requires a distribution over the randomness source and a translation from the randomness source to random choices in the program, even though we only care about the distribution of the latter.

### 7 Conclusion

In this work we presented an expressive probabilistic programming language that supports important features such as mixing continuous and discrete distributions, arrays, observations, partial functions and while-loops. Unlike prior work, our semantics distinguishes non-termination, observation failures and error states. This allows us to investigate the subtle interaction of different exceptions, which is not possible for semantics that conflate different kinds of exceptions. Our investigation confirms the intuitive understanding of the interaction of exceptions presented in Sect. 2. However, it also shows that some desirable properties, like commutativity, only hold in the absence of exceptions. This situation is unavoidable, and largely analogous to the situation in deterministic languages.

Even though our semantics only distinguish three exception states, it can be trivially extended to handle any countable set of exception states. This allows for an even finer-grained distinction of e.g. division by zero, out of bounds array accesses or casting failures (in a language that allows type casting). Our semantics also allows enriching exceptions with the line number that the exception originated from (of course, this is not possible for non-termination). For an uncountable set of exception states, an extension is possible but not trivial.

### A Proofs for Preliminaries

In this section, we provide lemmas, proofs and some definitions that were left out or cut short in Sect. 3. For a more detailed introduction into measure theory, we recommend the book *A crash course on the Lebesgue integral and measure theory* [7].

### A.1 Measures

Definition 3. *Let* (A, ΣA) *be a measurable space and* μ: ΣA → [0,∞] *a measure on* A*.*

Note that for a σ-finite measure μ, μ(A) = ∞ is possible, even though μ(Ai) < ∞ for all i. As an example, the Lebesgue measure is σ-finite because **R** = i∈**N**[−i, i] with λ([−i, i]) = 2 ∗ i, but λ(**R**) = ∞.

Lemma 12. *The following definition of s-finite measures is equivalent to our definition of s-finite measures (the difference is that the* μi*s are only required to be finite):*

*Proof.* Since any sub-probability measure is finite, one direction is trivial. For the other direction, let μ = - i∈**N** μ i for finite measures μ i. Obviously, μ ≥ 0, μ(∅)=0 and μ( i∈**N** Ai) = - i∈**N** Ai for mutually disjoint Ai ∈ ΣA, so μ is a measure. To show that μ can be written as a sum of sub-probability measures, let ni := μ i(A). Then, μ = - i∈**N** μ i = - i∈**N** ni ni μ i = - i∈**N** - j∈[ni] 1 ni μ i. We let μi := 1 ni μ i ≤ 1.

Lemma 13. *Any* σ*-finite measure* μ: ΣA → [0,∞] *is s-finite.*

*Proof.* Since μ is σ-finite, A = i∈**N** Ai with Ai ∈ ΣA and μ(Ai) < ∞. Without loss of generality, assume that the - Ai form a partition of A. Then, μ(S) = i∈**N** μ(S ∩ Ai), with μ(· ∩ Ai) < ∞. Thus, μ is a countable sum of finite measures.

Definition 4. *The counting measure* c : B → [0,∞] *is defined by*

$$c(S) = \begin{cases} |S| & S \ finite \\ \infty & otherwise \end{cases}$$

Definition 5. *The infinity measure* μ: B → [0,∞] *is defined by*

$$\mu(S) = \begin{cases} 0 & S = \emptyset \\ \infty & otherwise \end{cases}$$

Lemma 14. *Neither the counting measure nor the infinity measure are s-finite.*

*Proof.* For the counting measure c, assume (toward a contradiction) c = - i∈**N** ci. We have **R** = {r ∈ **R** | c({r}) > 0} = i∈**N**{r ∈ **R** | ci({r}) > 0} = i∈**N** n∈**N**{r ∈ **R** | ci({r}) > 1 n }. Because **R** is uncountable, there must be i, n ∈ **N** for which S := {r ∈ **R** | ci({r}) > 1 n } is uncountable. Thus for any measurable, countably infinite S ⊆ S, ci(S ) = ∞, which means that ci is not finite. Proceed analogously for the infinity measure.

Lemma 15. *The measure* μ : B → [0,∞] *with* μ(S) = 0 λ(S)=0 ∞ λ(S) > 0 *is sfinite but not* σ*-finite.*

*Proof.* μ = - i∈**N** λ, and λ is s-finite, so μ is s-finite. Assume (toward a contradiction) that μ is σ-finite. Then **R** = i∈**N** Ai with Ai ∈ B and μ(Ai) < ∞. Thus, μ(Ai)=0 and hence μ(**R**) = μ( i∈**N** Ai) ≤ - i∈**N** μ(Ai)=0, a contradiction.

Lemma 16

$$\begin{aligned} \forall S \in \Sigma\_{A \times B} \colon (\mu \times \mu')(S) &= \int\_{a \in A} \mu'(\{b \in B \mid (a, b) \in S\}) \mu(da) \\ &= \int\_{b \in B} \mu(\{a \in A \mid (a, b) \in S\}) \mu'(db) \\\\ \forall S \in \Sigma\_{\overline{A \times B}} \colon (\mu \overline{\times} \mu')(S) &= \int\_{a \in \overline{A}} \mu'(\{b \in \overline{B} \mid \overline{(a, b)} \in S\}) \mu(da) \end{aligned}$$

$$\begin{aligned} \forall S \in \Sigma\_{\overline{A \times B}} \colon (\mu \overline{\times} \mu')(S) &= \int\_{a \in \overline{A}} \mu'(\{b \in B \mid (a, b) \in S\}) \mu(da) \\ &= \int\_{b \in \overline{B}} \mu(\{a \in \overline{A} \mid \overline{(a, b)} \in S\}) \mu'(db) \end{aligned}$$

*Proof*

$$\begin{aligned} (\mu \times \mu')(S) &= \int \int\_{a \in A} \int\_{b \in B} [(a, b) \in S] \mu'(db) \mu(da) \\ &= \int \int\_{a \in A} \int \left[ b \in \{ b' \in B \mid (a, b') \in S \} \right] \mu'(db) \mu(da) \\ &= \int\_{a \in A} \mu'(\{ b' \in B \mid (a, b') \in S \}) \mu(da) \end{aligned}$$

172 B. Bichsel et al.

$$\begin{aligned} (\mu \times \mu')(S) &= \int\limits\_{a \in A} \int\limits\_{b \in B} [(a, b) \in S] \mu'(db) \mu(da) \\ &= \int\limits\_{b \in B} \int\limits\_{a \in A} [(a, b) \in S] \mu(da) \mu'(db) \\ &= \dots \\ &= \int\limits\_{b \in B} \mu(\{a' \in A \mid (a', b) \in S\}) \mu'(db) \end{aligned} \text{Fubini}$$

In the second line, we have used that (a, b) ∈ S ⇐⇒ b ∈ {b ∈ B | (a, b ) ∈ S}. The proof works analogously for ×.

Lemma 17. *Let* δ : A → A*,* κ: A → B*. Then,*

$$(\delta \overline{\times} \kappa)(a)(S) = \kappa(a)(\{b \in \overline{B} \mid \overline{(a,b)} \in S\})$$

*Proof*

$$\begin{aligned} (\delta \overline{\times} \kappa)(a)(S) &= \int\_{b \in \overline{B}} \delta(a)(\{a' \in A \mid \overline{(a',b)} \in S\}) \kappa(a)(db) & \text{Lemma 16} \\ &= \int\_{b \in \overline{B}} \overline{[(a,b) \in S]} \kappa(a)(db) \\ &= \kappa(a)(\{b \in \overline{B} \mid \overline{(a,b)} \in S\}) \end{aligned}$$

Lemma 1. *For measures* μ: ΣA → [0,∞]*,* μ : ΣB → [0,∞]*, let* S ∈ ΣA *and* T ∈ ΣB*. Then,* (μ × μ )(S × T) = μ(S) · μ (T)*.*

*Proof*

$$\begin{aligned} (\mu \times \mu')(S \times T) &= \int\_{a \in A} \mu'(\{b \in B \mid (a, b) \in S \times T\}) \mu(da) & \text{Lemma 16} \\ &= \int\_{a \in A} \mu' \left( \left\{ \begin{aligned} T \ a \in S \\ \emptyset \text{ otherwise} \end{aligned} \right\} \right) \mu(da) \\ &= \int\_{a \in S} \mu'(T) \mu(da) \\ &= \mu(S) \* \mu'(T) \end{aligned}$$

Lemma 2. × *and* × *for s-finite measures are associative, left- and right-distributive and preserve (sub-)probability and s-finite measures.*

*Proof.* Remember that (μ × μ )(S) = a∈A b∈B[(a, b) ∈ S]μ (db)μ(da) and that (μ×μ )(S) = a∈A b∈B[(a, b) ∈ S]μ (db)μ(da). Preservation of (sub-)probability measures is trivial. Distributivity and preservation of s-finite measures are easily established by properties of the Lebesgue integral in Lemma 19.

For associativity, let μ: ΣA → [0,∞], μ: ΣB → [0,∞] and μ: ΣC → [0,∞].

$$\begin{split} & (\mu \vee \mu') \times \mu')(S) \\ &= \int\_{c \in C} (\mu \times \mu')(\{t \in A \times B \mid (t, c) \in S\}) \mu''(dc) \qquad \text{Lemma 16} \\ &= \int\_{c \in C} \int\_{a \in A} \int\_{b \in B} [(a, b) \in \{t \in A \times B \mid (t, c) \in S\}] \mu'(db) \mu(da) \mu''(dc) \\ &= \int\_{c \in C} \int\_{a \in A} \int\_{b \in B} [(a, b, c) \in S] \mu'(db) \mu(da) \mu''(dc) \\ &= \int\_{a \in A} \int\_{b \in B} \int\_{c \in C} [(a, b, c) \in S] \mu''(dc) \mu'(db) \mu(da) \qquad \text{Fubini} \\ &= \int\_{a \in A} \int\_{b \in B} \int\_{c \in C} [(b, c) \in \{t \in B \times C \mid (a, t) \in S\}] \mu''(dc) \mu'(db) \mu(da) \\ &= \int\_{a \in A} (\mu' \times \mu'')(\{t \in B \times C \mid (a, t) \in S\}) \mu(da) \\ &= (\mu \times (\mu' \times \mu''))(S) \mu(da) \qquad \text{Lemma 16} \end{split}$$

The proof proceeds analogously for ×.

Lemma 18. *Let* (A, ΣA) *and* (B,ΣB) *be measurable spaces. Consider measures* μ, μ1, μ2 : ΣA → [0,∞] *and* ν, ν1, ν2 : ΣB → [0,∞]*. We assume that* ν1 ≤ ν2 *and* μ1 ≤ μ2 *hold pointwise. Then,*

$$
\begin{aligned}
\mu \overline{\times} \nu\_1 &\leq \mu \overline{\times} \nu\_2 \\
\mu\_1 \overline{\times} \nu &\leq \mu\_2 \overline{\times} \nu
\end{aligned}
$$

*Proof.* Let S ∈ ΣA×B and ν1 ≤ ν2. Then, we have

$$\begin{aligned} \nu\_1 &\le \nu\_2\\ \Longrightarrow \underbrace{\int\_{b\in\overline{B}} [\overline{(a,b)} \in S] \nu\_1(db)}\_{=:f(a)} &\le \underbrace{\int\_{b\in\overline{B}} [\overline{(a,b)} \in S] \nu\_2(db)}\_{=:g(a)} &\le \text{Lemma 19.}\\ \Longrightarrow \int\_{a\in\overline{A}} f(a)\mu(da) &\le \int\_{a\in\overline{A}} g(a)\mu(da) &\text{Lemma 19.}\\ \Longrightarrow (\mu\overline{\times}\nu\_1)(S) &\le (\mu\overline{\times}\nu\_2)(S) \end{aligned}$$

The proof for μ1×ν ≤ μ2×ν is similar.

#### A.2 Lebesgue Integral

Lemma 19.  The following properties can be proven for simple functions and limits of simple functions (this suffices):

$$\int\_{a \in E} f(a) \left( \sum\_{n=1}^{\infty} \mu\_i \right) (da) = \sum\_{n=1}^{\infty} \int\_{a \in E} f(a) \mu\_i(da)$$

$$\mu \le \nu \implies \int\_{a \in E} f(a) \mu(da) \le \int\_{a \in E} f(a) \nu(da)$$

$$\mu \le \mu(\nu) \text{ a constant } \mu \text{ is not equal}$$

 a∈E f(a)δ(x)(da) = f(x) is straightforward. For the other properties, see [31].

Theorem 1 (Fubini's theorem). *For s-finite measures* μ: ΣA → [0,∞] *and* μ : ΣB → [0,∞] *and any measurable function* f : A × B → [0,∞]*,*

$$\int\_{a \in A} \int\_{b \in B} f(a, b) \mu'(db) \mu(da) = \int\_{b \in B} \int\_{a \in A} f(a, b) \mu(da) \mu'(db)$$

*For s-finite measures* μ: ΣA → [0,∞] *and* μ : ΣB → [0,∞] *and any measurable function* f : A × B → [0,∞]*,*

$$\int\_{a \in \overline{A}} \int\_{b \in \overline{B}} \overline{f}(a,b)\mu'(db)\mu(da) = \int\_{b \in \overline{B}} \int\_{a \in \overline{A}} \overline{f}(a,b)\mu(da)\mu'(db)\mu'(db)$$

*Proof.* Let μ = - i∈**N** μi and μ = - i∈**N** μ i for bounded measures μi and μ i.

$$\begin{aligned} &\int\_{a\in A} \int\_{b\in B} f(a,b)\mu'(db)\mu(da) \\ &= \sum\_{i,j\in\mathbb{N}} \int\_{a\in A} \int\_{b\in B} f(a,b)\mu'\_j(db)\mu\_i(da) \quad \text{Lemma 19} \\ &= \sum\_{i,j\in\mathbb{N}} \int\_{b\in B} \int\_{a\in A} f(a,b)\mu\_i(da)\mu'\_j(db) \quad \text{Fubini for} \sigma\text{-finite measures } \mu\_i, \mu'\_j \\ &= \int\_{b\in B} \int\_{a\in A} f(a,b)\mu(da)\mu'(db) \end{aligned}$$

The proof in the presence of exception state is analogous.

Lemma 20. *Fubini does not hold for the counting measure* c : B → [0,∞] *and the Lebesgue measure* λ: B → [0,∞] *(because* c *is not s-finite).*

*Proof*

$$\begin{aligned} \int\_{x \in [0,1]} \int\_{y \in [0,1]} [x = y] c(dy) \lambda(dx) &= \int\_{x \in [0,1]} 1 \lambda(dx) = 1 \\ \int\_{y \in [0,1]} \int\_{x \in [0,1]} [x = y] \lambda(dx) c(dy) &= \int\_{y \in [0,1]} 0 c(dy) = 0 \end{aligned}$$

#### A.3 Kernels

Lemma 21. *Let* κ1, κ 1 : A → B *and* κ2, κ 2 : B → C *be s-finite kernels. If* κ1 ≤ κ 1 *holds pointwise, then*

$$
\kappa\_1 \gg \gg \kappa\_2 \le \kappa'\_1 \gg \gg \kappa\_2
$$

*If* κ2 ≤ κ 2 *holds pointwise, then*

$$
\kappa\_1 \gg \gg \kappa\_2 \le \kappa\_1 \gg \dots \gg \kappa\_2'
$$

*Proof.* Assume κ2 ≤ κ 2. Thus, κ2 ≤ κ 2. Now, let a ∈ A, S ∈ ΣC .

$$\begin{aligned} (\kappa\_1 \gg \dots \gg \kappa\_2)(a)(S) &= \int\_{b \in \overline{B}} \overline{\kappa\_2}(b)(S) \, \kappa\_1(a)(db) \\ &\le \int\_{b \in \overline{B}} \overline{\kappa\_2'}(b)(S) \, \kappa\_1(a)(db) \qquad \overline{\kappa\_2} \le \overline{\kappa\_2'}, \text{Lemma 19} \\ &= (\kappa\_1 \gg \dots \gg \kappa\_2')(a)(S) \end{aligned}$$

The proof for κ1 >=> κ2 ≤ κ 1 >=> κ2 works analogously.

Lemma 3. (;) *is associative, left- and right-distributive, has neutral element*4 δ *and preserves (sub-)probability and s-finite kernels.*

4 δ is a neutral element of (;) if (δ;κ)=(κ;δ) = κ for all kernels κ.

*Proof.* Remember that (f;g)(a)(S) = b∈B g(b)(S) f(a)(db). Left- and right-distributivity and the neutral element δ follow from properties of the Lebesgue integral in Lemma 19.

Associativity and preservation of (sub-)probability kernels is well known (see for example [12]). For s-finite kernels f = - i∈**N** fi and g = - i∈**N** - gi and h = i∈**N** hi, we have (for sub-probability kernels fi, gi, hi)

$$\begin{aligned} (f;g);h &= \left(\left(\sum\_{i\in\mathbb{N}} f\_i\right); \left(\sum\_{j\in\mathbb{N}} g\_j\right)\right); \sum\_{k\in\mathbb{N}} h\_k = \sum\_{i,j,k\in\mathbb{N}} (f\_i;g\_j);h\_k\\ &= \sum\_{i,j,k\in\mathbb{N}} f\_i; (g\_j;h\_k) = f;(g;h) \end{aligned}$$

(;) preserves s-finite kernels because for s-finite kernels f and g, we have (for sub-probability kernels fi, gi) f;g = - i,j∈**N** fi;gi, a sum of kernels.

Lemma 4. *For* f : A → B *and* g : B → C*,* a ∈ A *and* S ∈ ΣC *,*

$$(f \succ \succ g)(a)(S) = (f;g)(a)(S) + \sum\_{x \in \mathcal{X}} \delta(x)(S)f(a)(\{x\})$$

*Proof*

$$\begin{aligned} (f \underset{\succeq \mathcal{B}}{\geqslant} g)(a)(S) &= \int\_{b \in \mathcal{B}} \overline{g}(b)(S) \, f(a)(db) \\ &= \int\_{b \in B} \overline{g}(b)(S) \, f(a)(db) + \int\_{b \in \mathcal{X}} \overline{g}(b)(S) \, f(a)(db) \\ &= \int\_{b \in B} g(b)(S) \, f(a)(db) + \sum\_{b \in \mathcal{X}} \overline{g}(b)(S) \, f(a)(\{x\}) \\ &= (f;g)(a)(S) + \sum\_{x \in \mathcal{X}} \delta(x)(S) f(a)(\{x\}) \end{aligned}$$

Lemma 5. >=> *is associative, left-distributive (but not right-distributive), has neutral element* δ *and preserves (sub-)probability and s-finite kernels.*

*Proof.* Remember that (f >=> g)(a)(S) = b∈B g(b)(S) f(a)(db). Left-distributivity follows from the properties of the Lebesgue integral in Lemma 19. Rightdistributivity does not necessarily hold because g1 + g2(⊥) = g1(⊥) + g2(⊥). Associativity for f : A → B, g : B → C and h: C → D can be derived by

((f >=> g) >=> h)(a)(S) = f >=> g ;h (a)(S) +  x∈X δ(x)(S)(f >=> g)(a)({x}) = f;g + λa .λS . x∈X δ(x)(S )f(a )({x}) ;h (a)(S) +  x∈X δ(x)(S)(f >=> g)(a)({x}) = (f;g;h)(a)(S) + λa .λS . x∈X δ(x)(S )f(a )({x}) ;h (a)(S) =0((;)integrates over non-exception states) +  x∈X δ(x)(S)(f >=> g)(a)({x}) = (f;g;h)(a)(S) +  x∈X δ(x)(S) (f;g)(a)({x}) +  x-∈X δ(x )({x})f(a)({x }) = (f;g;h)(a)(S) +  x∈X δ(x)(S) (f;g)(a)({x}) + f(a)({x}) = (f;g;h)(a)(S) +  x∈X δ(x)(S)(f;λa .λS .g(a )(S ))(a)({x}) +  x∈X δ(x)(S)f(a)({x}) = (f;g;h)(a)(S) + f; λa .λS . x∈X δ(x)(S )g(a )({x}) (a)(S) +  δ(x)(S)f(a)({x})

$$\begin{aligned} &= \left( f; \left( g; h + \lambda a'. \lambda S'. \sum\_{x \in \mathcal{X}} \delta(x) (S') g(a') (\{x\}) \right) \right) (a) (S) + \sum\_{x \in \mathcal{X}} \delta(x) (S) f(a) (\{x\}) \\ &= \left( f; \left( g \gg > h \right) \right) (a) (S) + \sum\_{x \in \mathcal{X}} \delta(x) (S) f(a) (\{x\}) \\ &= (f \gg \gg (g \gg \gg h)) (a) (S) \end{aligned}$$

Here, we have used Lemma 4, left- and right-distributivity of (;).

To show that f >=> g preserves s-finite kernels, let f : A → B and g : B → C be s-finite kernels. Then, for sub-probability kernels fi,

$$\begin{aligned} (f \gg \gg g)(a)(S) &= (f;g)(a)(S) + \sum\_{x \in \mathcal{X}} \delta(x)(S)f(a)(\{x\}) \\ &= (f;g)(a)(S) + \sum\_{x \in \mathcal{X}} \sum\_{i \in \mathbb{N}} \delta(x)(S)f\_i(a)(\{x\}) \end{aligned}$$

Note that for each x ∈ X and i ∈ **N**, λa.λS.δ(x)(S)fi(a)({x}) is a sub-probability kernel. Thus, f >=> g is a sum of s-finite kernels and hence s-finite.

Proving that for sub-probability kernels f and g, f >=> g is also a (sub-) probability kernel is trivial, since we only need to show that (f >=> g)(a)(C)=1 (or ≤ 1).

Lemma 22. *Let* (A, ΣA) *and* (B,ΣB) *be measurable spaces. Let* f : A × B → [0,∞] *be measurable and* κ: A → B *be a sub-probability kernel. Then,* f : A → [0,∞] *defined by*

$$f'(a) := \int\_{b \in B} f(a, b) \kappa(a) (db)$$

*is measurable.*

*Proof.* See Theorem 20 of [30].

Lemma 23. × *and* × *preserve (sub-)probability kernels.*

*Proof.* Let κ: A → B and κ : A → C be (sub-)probability kernels. The fact that (κ × κ )(a)(·) for all a ∈ A is a (sub-)probability measure is inherited from Lemma 2. It remains to show that (κ×κ )(·)(S) is measurable for all S ∈ ΣB×C , with

$$(\kappa \times \kappa')(a)(S) = \int\_{b \in B} \int\_{c \in C} [(b, c) \in S] \kappa'(a)(dc) \kappa(a)(db)$$

By Lemma 22, f : A×B → [0,∞] defined by f (a, b) = c∈C [(b, c) ∈ S]κ (a)(dc) is measurable, using the measurable function f : (A × B) × C → [0,∞] defined by f((a, b), c) = [(b, c) ∈ S]. Again by Lemma 22, b∈B c∈C [(b, c) ∈ S]κ (a)(dc)κ(a)(db) is measurable.

Proving that for (sub-)probability kernels κ: A → B and κ : A → C, κ×κ is a (sub-)probability kernel proceeds analogously.

Lemma 6. × *and* × *for kernels preserve (sub-)probability and s-finite kernels, are associative, left- and right-distributive.*

*Proof.* Associativity, left- and right-distributivity are inherited from respective properties of the product of measures established by Lemma 2. Sub-probability kernels are preserved by Lemma 23.

S-finite kernels are preserved because κ × κ = (- i∈**N** κi) × ( - i∈**N** κ i - ) = i,j∈**N** κi × κ j (analogously for ×).

### B Proofs for Semantics

Lemma 7. *For* Δ *as in the semantics of the while loop, and for each* σ *and each* S*, the limit* limn→∞ Δn(-)(σ)(S) *exists.*

*Proof.* In general, 0 ≤ Δn(-)(σ)(S) ≤ 1. First, we restrict the allowed arguments for limn→∞ Δn(-)(σ)(S) to only those S with -∈ S. We prove by induction that Δn+1(-) ≤ Δn(-), meaning ∀σ : ∀S : -∈ S =⇒ Δn+1(-)(σ)(S) ≤ Δn(-)(σ)(S). Hence, Δn(-) is monotone decreasing in n and lower bounded by 0, which means that the limit must exist.

As a base case, we have Δ1(-)(σ)(S) ≤ 1 = δ-(S) = Δ0(-)(σ)(S), because -∈ S. We proceed by induction with

$$\begin{split} \Delta^{n+1}(\mathsf{O})(\sigma)(S) &= \left(\delta \overline{\times} [e] \right) \underset{\delta}{\iff} \lambda(\sigma, b). \left\{ \begin{bmatrix} P \\ \delta(\sigma) \end{bmatrix} \begin{matrix} b \not\equiv \Delta^{n}(\mathsf{O}) \ b \not\equiv 0 \\ b = 0 \end{matrix} \right\}(\sigma)(S) \\ &\leq \left(\delta \overline{\times} [e] \underset{\Rightarrow}{\iff} \lambda(\sigma, b). \left\{ \begin{matrix} [P](\sigma) \ \gg = \Delta^{n-1}(\mathsf{O}) \ b \not\equiv 0 \\ \delta(\sigma) \end{matrix} \right\}(\sigma)(S) \\ &= \Delta^{n}(\mathsf{O})(\sigma)(S) \end{split} \right\}(\sigma)(S) \end{split}$$

In the second line, we have used the induction hypothesis. This application is valid because κ2 ≤ κ 2 implies κ1 >=> κ2 ≤ κ1 >=> κ 2 (Lemma 21).

We proceed analogously when we restrict the allowed arguments for the kernel limn→∞ Δn(-)(σ)(S) to only those S with -∈/ S, proving Δn+1(-) ≥ Δn(-) for that case.

Lemma 8. *In the absence of exception states, and using sub-probability kernels instead of distribution transformers, the definition of the semantics of the while loop from [23] is equivalent to ours.*

Definition 6. *In [23], Kozen shows a different way of defining the semantics of the while loop. In our notation, and in terms of probability kernels instead of distribution transformers, that definition becomes*

$$\left[\mathsf{while1}\mathsf{i}\mathsf{e}\ \{P\}\right] = \sup\_{n\in\mathbb{N}}\sum\_{k=0}^{n}\left(\left[\mathsf{if1}\mathsf{i}\mathsf{i}\mathsf{i}\mathsf{i}\mathsf{e}\right]\right) \mathrel{\geqslant} \Longrightarrow \left[P\right]\right)^{k} \mathrel{\geqslant} \mathrel{\geqslant} \left[\mathsf{if1}\mathsf{i}\mathsf{i}\mathsf{i}\mathsf{e}\left(\neg e\right)\right]$$

*Here, exponentiation is in terms of Kleisli composition, i.e.* κ0 = δ *and* κn+1 = κ >=> κn*. The sum and limit are meant pointwise. Furthermore, we define filter by the following expression (note that* [[*filter*(e)]] *and* [[*filter*(¬e)]] *are only sub-probability kernels, not probability kernels).*

$$\begin{aligned} \[\mathsf{filt}\,\mathsf{ter}(e)\} &= \delta \overline{\times} [e] \; \rhd \Rightarrow \lambda(\sigma, b) . \begin{cases} \delta(\sigma) \; b \neq 0 \\ \mathbf{0} & b = 0 \end{cases} \\\ [\mathsf{filt}\,\mathsf{ter}(\neg e)] &= \delta \overline{\times} [e] \; \rhd \Rightarrow \lambda(\sigma, b) . \begin{cases} \delta(\sigma) \; b = 0 \\ \mathbf{0} & b \neq 0 \end{cases} \end{aligned}$$

To justify Lemma 8, we prove the more formal Lemma 24. Note that in the presence of exceptions (e.g. P is just **assert**(0)), Definition 6 does not make sense, because if

Lemma 24. *For all* S *with* S ∩ X = ∅

$$\left(\sum\_{k=0}^{n} \left(\left[\mathsf{f11\mathsf{ter}}(e)\right] >> \left[\mathsf{I}P\right]\right)^{k} >> \left[\mathsf{f11\mathsf{ter}}(\neg e)\right]\right)(\sigma)(S) = \Delta^{n+1}(\mathsf{O})(\sigma)(S)$$

*Proof.* For n = 0, we have

$$\begin{split} & \left( \sum\_{k=0}^{0} \left( [\mathtt{filter}(e)] \asymp \triangleright [P] \right)^{k} \asymp \triangleright [\mathtt{filter}(\neg e)] \right) (\sigma)(S) \\ &= \left( \left( [\mathtt{filter}(e)] \asymp \triangleright [P] \right)^{0} \asymp \triangleright [\mathtt{filter}(\neg e)] \right) (\sigma)(S) \\ &= \left( \delta \asymp \triangleright [\mathtt{filter}(\neg e)] \right) (\sigma)(S) \\ &= [\mathtt{filter}(\neg e)] (\sigma)(S) \\ &= \left( \delta \overline{\times} [e] \asymp \triangleright \lambda(\sigma', b) . \left\{ \begin{array}{l} \delta(\sigma') \ b = 0 \\ \mathbf{0} \end{array} \right\} (\sigma)(S) \\ &= \left( \delta \overline{\times} [e] \asymp \triangleright \lambda(\sigma', b) . \left\{ \begin{array}{l} \delta(\sigma') \ b = 0 \\ \mathbf{0} \left( \sigma' \right) \ b \neq 0 \end{array} \right\} (\sigma)(S) \\ &= \left( \delta \overline{\times} [e] \asymp \triangleright \lambda(\sigma', b) . \left\{ \begin{array}{l} \delta(\sigma') \qquad b = 0 \\ ([P] \ \asymp \odot \mathbf{O})(\sigma') \ b \neq 0 \end{array} \right\} (\sigma)(S) \\ &= \Delta^{1}(\mathbf{O}) \\ &= \Delta^{1}(\mathbf{O}) \end{split}$$

For n ≥ 0, we have

n +1 k=0 [[**filter**(e)]] >=> [[P ]]k >=> [[**filter**(¬e)]] (σ)(S) = n k=0 [[**filter**(e)]] >=> [[P ]]k+1 + ([[**filter**(e)]] >=> P ) 0 >=> [[**filter**(¬e)]] (σ)(S) = n k=0 [[**filter**(e)]] >=> [[P ]]k+1 + δ >=> [[**filter**(¬e)]] (σ)(S) = n k=0 [[**filter**(e)]] >=> [[P ]]k+1 >=> [[**filter**(¬e)]] (σ)(S) since S ∩ X = ∅ + (δ >=> [[**filter**(¬e)]]) (σ)(S) = n k=0 [[**filter**(e)]] >=> [[P ]]k+1 >=> [[**filter**(¬e)]] (σ)(S) + [[**filter**(¬e)]](σ)(S) = [[**filter**(e)]] >=> [[P ]] >=> n k=0 ([[**filter**(e)]] >=> [[P ]])k >=> [[**filter**(¬e)]] (σ)(S) + [[**filter**(¬e)]](σ)(S) = [[**filter**(e)]] >=> [[P ]] >=> n k=0 ([[**filter**(e)]] >=> [[P ]])k >=> [[**filter**(¬e)]] (σ)(S) + [[**filter**(¬e)]](σ)(S) = [[**filter**(e)]] >=> [[P ]] >=> Δn+1(-) (σ)(S) + [[**filter**(¬e)]](σ)(S) = δ×[[e]] >=> λ(σ- , b). [[P ]](σ- ) = Δn+1(-) b = 0 δ(σ- ) b = 0 (σ)(S) = Δn+2(-)(σ)(S)

In particular, have have used that left-distributivity does hold in this case since S ∩ X = ∅.

### C Probability Kernel

In the following, we list lemmas that are crucial to prove Theorem 2 (restated for convenience).

Theorem 2. *The semantics of each expression* [[e]] *and statement* [[P]] *is indeed a probability kernel.*

Lemma 25. *Any measurable function* f : A → [0,∞] *can be viewed as an sfinite kernel* f : A → **1***, defined by* f(x)(∅)=0 *and* f(x)(**1**) = f(x)*.*

*Proof.* We prove that f is an s-finite kernel. Let A∞ := {x ∈ A | f(x) = ∞}. Since f is measurable, the set A∞ must be measurable. f(x)(S) = - i∈**N**[x ∈ A∞][() ∈ S] + - i∈**N** f(x)[i ≤ f(x) < i + 1][() ∈ S], which is a sum of finite kernels because the sets A∞ and {x | i ≤ f(x) < i + 1} = f −1([i, i + 1)) are measurable. Note that any sum of finite kernels can be rewritten as a sum of sub-probability kernels.

Lemma 26. *Let* κ : X → Y *and* κ : X → Y *be kernels, and* f : X → **R** *measurable. Then,*

$$\kappa(x)(S) = \begin{cases} \kappa'(x)(S) & \text{if } f(x) = 0\\ \kappa''(x)(S) & \text{otherwise} \end{cases}$$

*is a kernel.*

*Proof.* Let f=0(x) := [f(x) = 0], f=0(x) := [f(x) = 0]. Then, κ = f=0 × κ + f=0 × κ. Viewing f=0 and f=0 as kernels X → **1** immediately gives the desired result.

Lemma 27. *Let* (A, ΣA) *and* (B,ΣB) *be measurable spaces. Let* {Ai}i∈I *be a partition of* A *into measurable sets, for a countable set of indices* I*. Consider a function* f : A → B*. If* f|Ai : Ai → B *is measurable for each* i ∈ I*, then* f *is measurable.*

Lemma 28. *Let* f : A → B *be measurable. Then* κ: A → B *with* κ(a) = δ(f(a)) *is a kernel.*

The following lemma is important to show that the semantics of the while loop is a probability kernel.

Lemma 29. *Suppose* {κn}n∈**N** *is a sequence of (sub-)probability kernels* A → B*. Then, if the limit* κ = limn→∞ κn *exists, it is also a (sub-)probability kernel. Here, the limit is pointwise in the sense* ∀a ∈ A: ∀S ∈ ΣB : κ(a, S) = limn→∞ κn(a)(S)*.*

*Proof.* For every a ∈ A, κ(a, ·) is a measure, because the pointwise limit of finite measures is a measure. For every S ∈ ΣB, κ(·, S) is measurable, because the pointwise limit of measurable functions fn : A → **R** (with B as the σ-algebra on **R**) is measurable.

### D Proofs for Consequences

In this section, we provide some proofs of consequences of our semantics, explained in Sect. 5.

Lemma 9. *For function* F(){**while** 1 {**skip**}; **return** 0}*,*

$$\frac{1}{0} + F() \not\equiv F() + \frac{1}{0}$$

*Proof.* If we evaluate 1 0 first, we will only have weight on ⊥.

$$\begin{aligned} &\left[\frac{1}{0} + F()\right] \\ &= \left[\frac{1}{0}\right] \overline{\times} [F()] \gg \ge \lambda(x, y). \delta(x + y) \\ &= \delta(\perp) \overline{\times} [F()] \gg \gg \lambda(x, y). \delta(x + y) \\ &= \delta(\perp) \gg \gg \lambda(x, y). \delta(x + y) \\ &= \delta(\perp) \end{aligned}$$

If instead, we first evaluate F(), we only have weight on -, by an analogous calculation.

Lemma 10. *If* [[e1]](σ)(X ) = [[e2]](σ)(X )=0 *for all* σ*, then* e1 ⊕ e2 e2 ⊕ e1*, for any commutative operator* ⊕*.*

*Proof*

$$\begin{aligned} [e\_1 \oplus e\_2](\sigma)(S) &= [e\_1] \overline{\times} [e\_2] \ \rangle \ \rangle \ \geqslant \ \lambda(x, y) \ \delta(x \oplus y) \\ &= \int\_{z \in \overline{\mathbb{R} \times \mathbb{R}}} \overline{\lambda(x, y) \cdot \delta(x \oplus y)}(z) (S) ([e\_1] \overline{\times} [e\_2])(\sigma)(dz) \\ &= \int\_{(x, y) \in \mathbb{R} \times \mathbb{R}} \delta(x \oplus y)(S) ([e\_1] \times [e\_2])(\sigma)(d(x, y)) \\ &= \int\_{(y, x) \in \mathbb{R} \times \mathbb{R}} \delta(y \oplus x)(S) ([e\_2] \times [e\_1])(\sigma)(d(y, x)) \\ &= [e\_2 \oplus e\_1](\sigma)(S) \end{aligned}$$

Here, we crucially rely on the absence of exceptions (for the third equality) and Fubini's Theorem (for the fourth equality).

Lemma 11. e1 ⊕ (e2 ⊕ e3) (e1 ⊕ e2) ⊕ e3*, for any associative operator* ⊕*.* *Proof.* The important steps of the proof are the following.

$$\begin{aligned} \left[e\_1 \oplus (e\_2 \oplus e\_3)\right] &= \left[e\_1\right] \overline{\times} \left[e\_2 \oplus e\_3\right] \gg \implies \lambda(x,s).\delta(x \oplus s) \\ &= \left[e\_1\right] \overline{\times} \left(\left[e\_2\right] \overline{\times} \left[e\_3\right] \gg \implies \lambda(y,z).\delta(y \oplus z)\right) \gg \implies \lambda(x,s).\delta(x \oplus s) \\ &= \left[e\_1\right] \overline{\times} \left(\left[e\_2\right] \overline{\times} \left[e\_3\right]\right) \gg \implies \lambda(x,\left(y,z\right)).\delta(x \oplus y \oplus z) \\ &= \left(\left[e\_1\right] \overline{\times} \left[e\_2\right]\right) \overline{\times} \left[e\_3\right] \gg \implies \lambda((x,y),z).\delta(x \oplus y \oplus z) \\ &= \left[\left(e\_1 \oplus e\_2\right) \oplus e\_3\right] \end{aligned}$$

Here, we make crucial use of associativity for the lifted product of measures in Lemma 6.

### References

Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **How long, O Bayesian network, will I sample thee? A program analysis perspective on expected sampling times**

Kevin Batz(B) , Benjamin Lucien Kaminski(B) , Joost-Pieter Katoen(B) , and Christoph Matheja(B)

> RWTH Aachen University, Aachen, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>
**Abstract.** Bayesian networks (BNs) are probabilistic graphical models for describing complex joint probability distributions. The main problem for BNs is inference: Determine the probability of an event given observed evidence. Since exact inference is often infeasible for large BNs, popular approximate inference methods rely on sampling.

We study the problem of determining the expected time to obtain a single valid sample from a BN. To this end, we translate the BN together with observations into a probabilistic program. We provide proof rules that yield the exact expected runtime of this program in a fully automated fashion. We implemented our approach and successfully analyzed various real–world BNs taken from the Bayesian network repository.

**Keywords:** Probabilistic programs · Expected runtimes Weakest preconditions · Program verification

### **1 Introduction**

*Bayesian networks* (BNs) are *probabilistic graphical models* representing joint probability distributions of sets of random variables with conditional dependencies. Graphical models are a popular and appealing modeling formalism, as they allow to succinctly represent complex distributions in a human–readable way. BNs have been intensively studied at least since 1985 [43] and have a wide range of applications including machine learning [24], speech recognition [50], sports betting [11], gene regulatory networks [18], diagnosis of diseases [27], and finance [39].

*Probabilistic programs* are programs with the key ability to draw values at random. Seminal papers by Kozen from the 1980s consider formal semantics [32] as well as initial work on verification [33,47]. McIver and Morgan [35] build on this work to further weakest–precondition style verification for imperative probabilistic programs.

c The Author(s) 2018 A. Ahmed (Ed.): ESOP 2018, LNCS 10801, pp. 186–213, 2018. https://doi.org/10.1007/978-3-319-89884-1\_7

The interest in probabilistic programs has been rapidly growing in recent years [20,23]. Part of the reason for this d´ej`a vu is their use for representing probabilistic graphical models [31] such as BNs. The full potential of modern probabilistic programming languages like Anglican [48], Church [21], Figaro [44], R2 [40], or Tabular [22] is that they enable rapid prototyping and obviate the need to manually provide inference methods tailored to an individual model.

*Probabilistic inference* is the problem of determining the probability of an event given observed evidence. It is a major problem for both BNs and probabilistic programs, and has been subject to intense investigations by both theoreticians and practitioners for more than three decades; see [31] for a survey. In particular, it has been shown that for probabilistic programs exact inference is highly undecidable [28], while for BNs both *exact inference* as well as *approximate inference* to an arbitrary precision are NP–hard [12,13]. In light of these complexity– theoretical hurdles, a popular way to analyze probabilistic graphical models as well as probabilistic programs is to gather a large number of independent and identically distributed (i.i.d. for short) samples and then do statistical reasoning on these samples. In fact, all of the aforementioned probabilistic programming languages support sampling based inference methods.

*Rejection sampling* is a fundamental approach to obtain valid samples from BNs with observed evidence. In a nutshell, this method first samples from the joint (unconditional) distribution of the BN. If the sample complies with all evidence, it is valid and accepted; otherwise it is rejected and one has to resample.

Apart from rejection sampling, there are more sophisticated sampling techniques, which mainly fall in two categories: Markov Chain Monte Carlo (MCMC) and importance sampling. But while MCMC requires heavy hand–tuning and suffers from slow convergence rates on real–world instances [31, Chapter 12.3], virtually all variants of importance sampling rely again on rejection sampling [31,49].

A major problem with rejection sampling is that for poorly conditioned data, this approach might have to reject and resample very often in order to obtain just a single accepting sample. Even worse, being poorly conditioned need not be immediately evident for a given BN, let alone a probabilistic program. In fact, Gordon et al. [23, p. 177] point out that

"the main challenge in this setting [i.e. sampling based approaches] is that many samples that are generated during execution are ultimately rejected for not satisfying the observations."

If too many samples are rejected, the expected sampling time grows so large that sampling becomes infeasible. The expected sampling time of a BN is therefore a key figure for deciding whether sampling based inference is the method of choice.

*How Long, O Bayesian Network, will I Sample Thee?* More precisely, we use techniques from program verification to give an answer to the following question:

Given a Bayesian network with observed evidence, how long does it take in expectation to obtain a *single* sample that satisfies the observations?

**Fig. 1.** A simple Bayesian network.

As an example, consider the BN in Fig. 1 which consists of just three nodes (random variables) that can each assume values 0 or 1. Each node X comes with a conditional probability table determining the probability of X assuming some value given the values of all nodes Y that X depends on (i.e. X has an incoming edge from Y ), see [3, Appendix A.1] for detailed calculations. For instance, the probability that G assumes value 0, given that S and R are both assume 1, is 0.2. Note that this BN is paramterized by a ∈ [0, 1].

Now, assume that our observed evidence is the event G = 0 and we apply rejection sampling to obtain *one* accepting sample from this BN. Then our approach will yield that a rejection sampling algorithm will, on average, require

$$\frac{200a^2 - 40a - 460}{89a^2 - 69a - 21}$$

guard evaluations, random assignments, etc. until it obtains a single sample that complies with the observation G = 0 (the underlying runtime model is discussed in detail in Sect. 3.3). By examination of this function, we see that for large ranges of values of a the BN is rather well–behaved: For a ∈ [0.08, 0.78] the expected sampling time stays below 18. Above a = 0.95 the expected sampling time starts to grow rapidly up to 300.

While 300 is still moderate, we will see later that expected sampling times of real–world BNs can be much larger. For some BNs, the expected sampling time even exceeded 1018, rendering sampling based methods infeasible. In this case, exact inference (despite NP–hardness) was a viable alternative (see Sect. 6).

*Our Approach.* We apply weakest precondition style reasoning a l´a McIver and Morgan [35] and Kaminski et al. [30] to analyze both expected outcomes and *expected runtimes* (ERT) of a *syntactic fragment of* pGCL, which we call the *Bayesian Network Language* (BNL). Note that since BNL is a syntactic fragment of pGCL, every BNL program is a pGCL program but *not vice versa*. The main restriction of BNL is that (in contrast to pGCL) loops are of a special form that prohibits undesired data flow across multiple loop iterations. While this restriction renders BNL incapable of, for instance, counting the number of loop iterations1, BNL is expressive enough to encode Bayesian networks with observed evidence.

For BNL, we develop dedicated proof rules to determine *exact* expected values and the *exact* ERT of any BNL program, including loops, without any user– supplied data, such as invariants [30,35], ranking or metering functions [19], (super)martingales [8–10], etc.

As a central notion behind these rules, we introduce f–*i.i.d.–ness* of probabilistic loops, a concept closely related to stochastic independence, that allows us to *rule out undesired parts of the data flow across loop iterations*. Furthermore, we show how every BN with observations is translated into a BNLprogram, such that

As a consequence, exact expected sampling times of BNs can be inferred by means of weakest precondition reasoning in a fully automated fashion. This can be seen as a first step towards formally evaluating the quality of a plethora of different sampling methods (cf. [31,49]) on source code level.

*Contributions.* To summarize, our main contributions are as follows:

*Outline.* We discuss related work in Sect. 2. Syntax and semantics of the probabilistic programming language pGCL are presented in Sect. 3. Our proof rules are introduced in Sect. 4 and applied to BNs in Sect. 5. Section 6 reports on experimental results and Sect. 7 concludes.

1 An example of a program that is *not* expressible in BNL is given in Example 1.

### **2 Related Work**

While various techniques for formal reasoning about runtimes and expected outcomes of probabilistic programs have been developed, e.g. [6,7,17,25,38], none of them explicitly apply formal methods to reason about Bayesian networks on source code level. In the following, we focus on approaches close to our work.

*Weakest Preexpectation Calculus.* Our approach builds upon the expected runtime calculus [30], which is itself based on work by Kozen [32,33] and McIver and Morgan [35]. In contrast to [30], we develop specialized proof rules for a clearly specified program fragment *without* requiring user–supplied invariants. Since finding invariants often requires heavy calculations, our proof rules contribute towards simplifying and automating verification of probabilistic programs.

*Ranking Supermartingales.* Reasoning about almost–sure termination is often based on ranking (super)martingales (cf. [8,10]). In particular, Chatterjee et al. [9] consider the class of affine probabilistic programs for which linear ranking supermartingales exist (Lrapp); thus proving (positive2) almost–sure termination for all programs within this class. They also present a doubly–exponential algorithm to approximate ERTs of Lrapp programs. While all BNL programs lie within Lrapp, our proof rules yield *exact* ERTs as *expectations* (thus allowing for compositional proofs), in contrast to a single number for a fixed initial state.

*Bayesian Networks and Probabilistic Programs.* Bayesian networks are a—if not the most—popular probabilistic graphical model (cf. [4,31] for details) for reasoning about conditional probabilities. They are closely tied to (a fragment of) probabilistic programs. For example, Infer.NET [36] performs inference by compiling a probabilistic program into a Bayesian network. While correspondences between probabilistic graphical models, such as BNs, have been considered in the literature [21,23,37], we are not aware of a formal soudness proof for a translation from classical BNs into probabilistic programs including conditioning.

Conversely, some probabilistic programming languages such as Church [21], Stan [26], and R2 [40] directly perform inference on the program level using sampling techniques similar to those developed for Bayesian networks. Our approach is a step towards understanding sampling based approaches formally: We obtain the exact expected runtime required to generate a sample that satisfies all observations. This may ultimately be used to evaluate the quality of a plethora of proposed sampling methods for Bayesian inference (cf. [31,49]).

### **3 Probabilistic Programs**

We briefly present the probabilistic programming language that is used throughout this paper. Since our approach is embedded into weakest-precondition style approaches, we also recap calculi for reasoning about both expected outcomes and expected runtimes of probabilistic programs.

2 Positive almost–sure termination means termination in finite expected time [5].

#### **3.1 The Probabilistic Guarded Command Language**

We enhance Dijkstra's Guarded Command Language [14,15] by a probabilistic construct, namely a random assignment. We thereby obtain a *probabilistic Guarded Command Language* (for a closely related language, see [35]).

Let Vars be a finite set of *program variables*. Moreover, let Q be the set of rational numbers, and let D (Q) be the set of discrete probability distributions over Q. The set of *program states* is given by Σ = { σ | σ : Vars → Q }.

A *distribution expression* μ is a function of type μ: Σ → D (Q) that takes a program state and maps it to a probability distribution on values from Q. We denote by μσ the distribution obtained from applying σ to μ.

The probabilistic guarded command language (pGCL) is given by the grammar

where x ∈ Vars is a program variable, μ is a distribution expression, and ϕ is a Boolean expression guarding a choice or a loop. A pGCL program that contains neither diverge, nor while, nor repeat − until loops is called loop–free.

For σ ∈ Σ and an arithmetical expression E over Vars, we denote by σ(E) the evaluation of E in σ, i.e. the value that is obtained by evaluating E after replacing any occurrence of any program variable x in E by the value σ(x). Analogously, we denote by σ(ϕ) the evaluation of a guard ϕ in state σ to either true or false. Furthermore, for a value v ∈ Q we write σ [x → v] to indicate that we set program variable x to value v in program state σ, i.e.3

$$
\sigma \left[ x \mapsto v \right] \;= \; \lambda \, y \\
\bullet \; \begin{cases} v, & \text{if } y = x \\ \sigma(y), & \text{if } y \neq x \end{cases}
$$

We use the Iverson bracket notation to associate with each guard its according indicator function. Formally, the Iverson bracket [ϕ] of ϕ is thus defined as the function [ϕ] = λ σ**.** σ(ϕ).

Let us briefly go over the pGCL constructs and their effects: skip does not alter the current program state. The program diverge is an infinite busy loop, thus takes infinite time to execute. It returns no final state whatsoever.

The random assignment x :≈ μ is (a) the only construct that can actually alter the program state and (b) the only construct that may introduce random

3 We use λ–expressions to construct functions: Function λX**.** applied to an argument α evaluates to in which every occurrence of X is replaced by α.

behavior into the computation. It takes the current program state σ, then *samples* a value v from probability distribution μσ, and then assigns v to program variable x. An example of a random assignment is

$$\{x \mathrel{\approx} 1/2 \cdot \langle 5 \rangle + 1/6 \cdot \langle y+1 \rangle + 1/3 \cdot \langle y-1 \rangle \ .\}$$

If the current program state is σ, then the program state is altered to either σ [x → 5] with probability 1/2, or to σ [x → σ(y) + 1] with probability 1/6, or to σ [x → σ(y) − 1] with probability 1/3. The remainder of the pGCL constructs are standard programming language constructs.

In general, a pGCL program C is executed on an input state and yields a *probability distribution* over final states due to possibly occurring random assignments inside of C. We denote that resulting distribution by -Cσ. Strictly speaking, programs can yield *subdistributions*, i.e. probability distributions whose total mass may be below 1. The "missing"probability mass represents the probability of nontermination. Let us conclude our presentation of pGCLwith an example:

*Example 1 (Geometric Loop).* Consider the program C*geo* given by

$$\begin{aligned} x &\approx 0; & c &\approx 1/2 \cdot \langle 0 \rangle + 1/2 \cdot \langle 1 \rangle; \\ \mathtt{while} \, \mathtt{e} \, (c &= 1) \left\{ x : \approx x + 1; \, c &\approx 1/2 \cdot \langle 0 \rangle + 1/2 \cdot \langle 1 \rangle \right\} \end{aligned}$$

This program basically keeps flipping coins until it flips, say, heads (c = 0). In x it counts the number of unsuccessful trials.4 In effect, it almost surely sets c to 0 and moreover it establishes a geometric distribution on x. The resulting distribution is given by

$$\left[\left\|C\_{geo}\right\|\_{\sigma}(\tau)\right] = \sum\_{n=0}^{\omega} \left[\tau = \sigma\left[c, x \mapsto 0, n\right]\right] \cdot \frac{1}{2^{n+1}}\,. \tag{7}$$

#### **3.2 The Weakest Preexpectation Transformer**

We now present the weakest preexpectation transformer wp for reasoning about expected outcomes of executing probabilistic programs in the style of McIver and Morgan [35]. Given a random variable f mapping program states to reals, it allows us to reason about the expected value of f after executing a probabilistic program on a given state.

**Expectations.** The random variables the wp transformer acts upon are taken from a set of so-called expectations, a term coined by McIver and Morgan [35]:

4 This counting is also the reason that C*geo* is an example of a program that is not expressible in our BNL language that we present later.

**Definition 1 (Expectations).** *The set of expectations* E *is defined as*

E = f f : Σ → R∞ ≥0 .

*We will use the notation* f[x/E] *to indicate the* replacement *of every occurrence of* x *in* f *by* E*. Since* x*, however, does not actually occur in* f*, we more formally define* f[x/E] = λσ*.* f(σ [x → σ(E)])*.*

*A complete partial order* ≤ *on* E *is obtained by point–wise lifting the canonical total order on* R∞ ≥0*, i.e.*

$$f\_1 \preceq f\_2 \quad \text{iff} \quad \forall \sigma \in \Sigma \colon \quad f\_1(\sigma) \le \ f\_2(\sigma) \dots$$

*Its least element is given by* λσ*.* 0 *which we (by slight abuse of notation) also denote by* 0*. Suprema are constructed pointwise, i.e. for* S ⊆ E *the supremum* sup S *is given by* sup S = λσ*.* supf∈S f(σ)*.*

We allow expectations to map only to positive reals, so that we have a complete partial order readily available, which would not be the case for expectations of type Σ → R ∪ {−∞, +∞}. A wp calculus that *can* handle expectations of such type needs more technical machinery and cannot make use of this underlying natural partial order [29]. Since we want to reason about ERTs which are by nature non–negative, we will not need such complicated calculi.

Notice that we use a slightly different definition of expectations than McIver and Morgan [35], as we allow for *unbounded* expectations, whereas [35] requires that expectations are *bounded*. This however would prevent us from capturing ERTs, which are potentially unbounded.

**Expectation Transformers.** For reasoning about the expected value of f ∈ E after execution of C, we employ a backward–moving weakest preexpectation transformer wp-C : E → E, that maps a *postexpectation* f ∈ E to a *preexpectation* wp -C (f) ∈ E, such that wp -C (f) (σ) is the expected value of f after executing C on initial state σ. Formally, if C executed on input σ yields final distribution -Cσ, then the *weakest preexpectation* wp -C (f) *of* C *with respect to postexpectation* f *is given by*

$$\text{wp}\left[\![C]\!\right](f)\left(\sigma\right) = \int\_{\Sigma} f \, d\left[\![C]\!\right]\_{\sigma} \,, \tag{1}$$

where we denote by A h dν the expected value of a random variable h: A → R∞ ≥0 with respect to a probability distribution ν : A → [0, 1]. Weakest preexpectations can be defined in a very systematic way:

**Definition 2 (The wp Transformer** [35]**).** *The weakest preexpectation transformer* wp: pGCL → E → E *is defined by induction on all* pGCL *programs according to the rules in Table 1. We call* Ff (X)=[¬ϕ] · f + [ϕ] · wp -C (X) *the* wp*–* characteristic functional *of the loop* while (ϕ) {C} *with respect to postexpectation* f*. For a given* wp*–characteristic function* Ff *, we call the sequence* {F n f (0)}n∈N *the* orbit of Ff *.*

**Table 1.** Rules for the wp–transformer.

Let us briefly go over the definitions in Table 1: For skip the program state is not altered and thus the expected value of f is just f. The program diverge will never yield any final state. The distribution over the final states yielded by diverge is thus the null distribution ν0(τ ) = 0, that assigns probability 0 to *every* state. Consequently, the expected value of f after execution of diverge is given by Σ f dν0 = τ∈Σ 0 · f(τ ) = 0.

The rule for the random assignment x :≈ μ is a bit more technical: Let the current program state be σ. Then for every value v ∈ Q, the random assignment assigns v to x with probability μσ(v), where σ is the current program state. The value of f after assigning v to x is f(σ [x → v]) = f[x/v](σ) and therefore the expected value of f after executing the random assignment is given by

$$\sum\_{v \in \mathbb{Q}} \mu\_{\sigma}(v) \cdot f[x/v](\sigma) \, = \int\_{\mathbb{Q}} \left(\lambda v \bullet f[x/v](\sigma)\right) \, d\mu\_{\sigma} \, .$$

Expressed as a function of σ, the latter yields precisely the definition in Table 1.

The definition for the conditional choice if (ϕ) {C1} else {C2} is not surprising: if the current state satisfies ϕ, we have to opt for the weakest preexpectation of C1, whereas if it does not satisfy ϕ, we have to choose the weakest preexpectation of C2. This yields precisely the definition in Table 1.

The definition for the sequential composition C1; C2 is also straightforward: We first determine wp -C2 (f) to obtain the expected value of f after executing C2. Then we mentally prepend the program C2 by C1 and therefore determine the expected value of wp -C2 (f) after executing C1. This gives the weakest preexpectation of C1; C2 with respect to postexpectation f.

The definition for the while loop makes use of a least fixed point, which is a standard construction in program semantics. Intuitively, the fixed point iteration of the wp–characteristic functional, given by 0, Ff (0), F2 f (0), F3 f (0), ..., corresponds to the portion the expected value of f after termination of the loop, that can be collected within at most 0, 1, 2, 3, ... loop guard evaluations. The Kleene Fixed Point Theorem [34] ensures that this iteration converges to the least fixed point, i.e.

$$\sup\_{n \in \mathbb{N}} F\_f^n(0) = \text{lfp } F\_f = \text{wp} \left[ \text{whi1e} \left( \varphi \right) \{ C \} \right](f) \dots$$

By inspection of the above equality, we see that the least fixed point is exactly the construct that we want for while loops, since supn∈N F n f (0) in principle allows the loop to run for any number of iterations, which captures precisely the semantics of a while loop, where the number of loop iterations is—in contrast to e.g. for loops—not determined upfront.

Finally, since repeat {C} until (ϕ) is syntactic sugar for C; while (ϕ) {C}, we simply define the weakest preexpectation of the former as the weakest preexpectation of the latter. Let us conclude our study of the effects of the wp transformer by means of an example:

*Example 2.* Consider the following program C:

$$\begin{aligned} c &\approx \mathbf{1} \langle 3 \cdot \langle 0 \rangle + 2 \langle 3 \cdot \langle 1 \rangle \rangle; \\ \mathbf{if} \,(c=0) &\{ x \approx 1/2 \cdot \langle 5 \rangle + 1/6 \cdot \langle y+1 \rangle + 1/3 \cdot \langle y-1 \rangle \} \mathbf{e1sse \{skip\}}. \end{aligned}$$

Say we wish to reason about the expected value of x + c after execution of the above program. We can do so by calculating wp -C (x + c) using the rules in Table 1. This calculation in the end yields wp -C (x + c) = 3y+26/18 The expected valuation of the expression x + c after executing C is thus 3y+26/18. Note that x + c can be thought of as an expression that is evaluated in the final states after execution, whereas 3y+26/18 must be evaluated in the initial state before execution of C.

**Healthiness Conditions of wp.** The wp transformer enjoys some useful properties, sometimes called *healthiness conditions* [35]. Two of these healthiness conditions that we will heavily make use of are given below:

**Theorem 1 (Healthiness Conditions for the wp Transformer** [35]**).** *For all* C ∈ pGCL*,* f1, f2 ∈ E*, and* a ∈ R≥0*, the following holds:*

#### **3.3 The Expected Runtime Transformer**

While for deterministic programs we can speak of *the* runtime of a program on a given input, the situation is different for probabilistic programs: For those we instead have to speak of the *expected runtime* (ERT). Notice that the ERT can be finite (even constant) while the program may still admit infinite executions. An example of this is the geometric loop in Example 1.

A wp–like transformer designed specifically for reasoning about ERTs is the ert transformer [30]. Like wp, it is of type ert-C : E → E and it can be shown that

**Table 2.** Rules for the ert–transformer.

ert -C (0) (σ) is precisely the *expected runtime of executing* C *on input* σ. More generally, if f : Σ → R∞ ≥0 measures the time that is needed after executing C (thus f is evaluated in the final states after termination of C), then ert -C (f) (σ) is the expected time that is needed to run C on input σ and then let time f pass. For a more in–depth treatment of the ert transformer, see [30, Sect. 3]. The transformer is defined as follows:

**Definition 3 (The ert Transformer** [30]**).** *The expected runtime transformer* ert: pGCL → E → E *is defined by induction on all* pGCL *programs according to the rules given in Table 2. We call* Ff (X) = 1+[¬ϕ]·f +[ϕ]·wp -C (X) *the* ert*–* characteristic functional *of the loop* while (ϕ) {C} *with respect to postexpectation* f*. As with* wp*, for a given* ert*–characteristic function* Ff *, we call the sequence* {F n f (0)}n∈N *the* orbit of Ff *. Notice that*

$$\text{vert}\,\text{[wh1e}\,\text{e}\,(\varphi)\,\{C\}\text{]}\,(f) \,=\,\text{|f\!p}\,F\_f \,=\,\text{sup}\,\,\{F\_f^n(0)\}\_{n\in\mathbb{N}}.$$

The rules for ert are very similar to the rules for wp. The runtime model we assume is that skip statements, random assignments, and guard evaluations for both conditional choice and while loops cost one unit of time. This runtime model can easily be adopted to count only the number of loop iterations or only the number of random assignments, etc. We conclude with a strong connection between the wp and the ert transformer, that is crucial in our proofs:

**Theorem 2 (Decomposition of ert** [41]**).** *For any* C ∈ pGCL *and* f ∈ E*,*

$$\text{ert}\begin{bmatrix} \left[ C \right] \left( f \right) \end{bmatrix} = \text{ert}\begin{bmatrix} \left[ C \right] \left( 0 \right) + \mathbf{w} \mathbf{p} \begin{bmatrix} \left[ C \right] \left( f \right) \end{bmatrix}.$$

### **4 Expected Runtimes of i.i.d Loops**

We derive a proof rule that allows to determine *exact ERTs of independent and identically distributed loops* (or *i.i.d. loops* for short). Intuitively, a loop

**Fig. 2.** An i.i.d. loop sampling a point within a circle uniformly at random using rejection sampling. The picture on the right–hand side visualizes the procedure: In each iteration a point (×) is sampled. If we obtain a point within the white area inside the square, we terminate. Otherwise, i.e. if we obtain a point within the gray area outside the circle, we resample.

is i.i.d. if the distributions of states that are reached at the end of different loop iterations are equal. This is the case whenever there is no data flow across different iterations. In the non–probabilistic case, such loops either terminate after exactly one iteration or never. This is different for probabilistic programs.

As a running example, consider the program C*circle* in Fig. 2. C*circle* samples a point within a circle with center (5, 5) and radius r = 5 uniformly at random using rejection sampling. In each iteration, it samples a point (x, y) ∈ [0,..., 10]2 within the square (with some fixed precision). The loop ensures that we resample if a sample is not located within the circle. Our proof rule will allow us to systematically determine the ERT of this loop, i.e. the average amount of time required until a single point within the circle is sampled.

Towards obtaining such a proof rule, we first present a syntactical notion of the i.i.d. property. It relies on expectations that are not affected by a pGCL program:

**Definition 4.** *Let* C ∈ pGCL *and* f ∈ E*. Moreover, let Mod* (C) *denote the set of all variables that occur on the left–hand side of an assignment in* C*, and let Vars*(f) *be the set of all variables that "occur in* f*", i.e. formally*

$$x \in \mathsf{Vars}(f) \qquad \text{iff} \qquad \exists \sigma \, \exists \, v, v' \colon \quad f(\sigma \, [x \mapsto v]) \neq \, f(\sigma \, [x \mapsto v']).$$

*Then* f *is* unaffected *by* C*, denoted* f -C*, iff Vars*(f) ∩ *Mod* (C) = ∅*.*

We are interested in expectations that are unaffected by pGCL programs because of a simple, yet useful observation: If g - C, then g *can be treated like a constant* w.r.t. the transformer wp (i.e. like the a in Theorem 1 (1)). For our running example C*circle* (see Fig. 2), the expectation f = wp -C*body* ([x + y ≤ 10]) is unaffected by the loop body C*body* of C*circle* . Consequently, we have wp -C*body* (f) = f · wp -C*body* (1) = f. In general, we obtain the following property: