suppReferences

# Expectation Programming

Tim Reichelt1, 2 Adam Goliński1 Luke Ong2 Tom Rainforth1

1 Department of Statistics, University of Oxford
2 Department of Computer Science, University of Oxford

###### Abstract

Building on ideas from probabilistic programming, we introduce the concept of
an _expectation programming framework_ (EPF) that automates the calculation of
expectations. Analogous to a probabilistic program, an expectation program is
comprised of a mix of probabilistic constructs and deterministic calculations
that define a conditional distribution over its variables. However, the focus
of the inference engine in an EPF is to directly estimate the resulting
expectation of the program return values, rather than approximate the
conditional distribution itself. This distinction allows us to achieve
substantial performance improvements over the standard probabilistic
programming pipeline by tailoring the inference to the precise expectation we
care about. We realize a particular instantiation of our EPF concept by
extending the probabilistic programming language _Turing_ to allow so-called
_target-aware_ inference to be run automatically, and show that this leads to
significant empirical gains compared to conventional posterior-based
inference.

##  1 Introduction

Calculating expectations is at the center of many scientific workflows. For
example, the decision theoretic foundations of most statistical paradigms,
e.g. Bayesian modeling, are rooted in calculating the expectation of a loss
function [1]. Moreover, these expectations are often taken with respect to
distributions that require inference to draw samples from or evaluate the
density of.

One such context is probabilistic programming [2, 3]. Here our program is
typically specified (often indirectly) through an unnormalized density
$\gamma(x)$. The role of system’s inference engine is now to approximate the
distribution specified by the normalized density $\pi(x)=\gamma(x)/Z$, where
$Z$ is an unknown normalizing constant. Here the standard pipeline to
calculate some expectation $\mathbb{E}_{\pi(x)}[f(x)]$ is to construct such an
approximation $\hat{\pi}(x)$ independently of $f(x)$, before using this
approximation to, in turn, estimate the expectation.

Sometimes $f(x)$ is not yet known at the point of inference, such that this
pipeline is inevitable. However, if $f(x)$ is known upfront, then we can
improve on this pipeline [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. In particular,
substantial gains can be achieved by using information in $f(x)$ to construct
target-aware Bayesian inference (TABI) estimators directly geared towards the
target expectation [14, 15].

At the moment, none of the contemporary probabilistic programming languages
(PPL) (e.g. Bingham et al. [16], Carpenter et al. [17], Cusumano-Towner et al.
[18], Ge et al. [19], Salvatier et al. [20], Tran et al. [21], Wood et al.
[22], Mansinghka et al. [23], Goodman et al. [24], Goodman and Stuhlmüller
[25], Tavares et al. [26], Fierens et al. [27], Murray and Schön [28],
Narayanan et al. [29], Minka et al. [30]) allow one to explicitly incorporate
the function $f(x)$ into their inference calculations for a given program.
Namely, they all implicitly rely on the standard pipeline of approximating
$\pi(x)$ and then separately using this approximation in turn to estimate
expectations. As such, no existing PPLs currently automate the estimation of
expectations in a _target-aware_ fashion.

To address this, we introduce, and formalize, the concept of an _expectation
programming framework_ (EPF): a framework that allows users to easily define
expectations directly, and then automate their estimation using target-aware
inference algorithms. In a manner analogous to how PPLs can be viewed as tools
for assisting and automating the process of approximate Bayesian inference,
EPFs assist and automate the estimation of conditional expectations. They
further provide a convenient and expressive mechanism for expressing such
problems in a programmatic manner.

To demonstrate this, we further introduce a specific implementation of an EPF,
called EPT (Expectation Programming in Turing), built upon the Turing PPL
[19]. EPT takes as input a Turing-style program and uses program
transformations to create a new set of three valid Turing programs that can be
used to construct target-aware estimators via the TABI approach of [15].
Estimates are generated for these new programs separately, before being
recombined into an overall estimate. As the transformations result in valid
Turing programs, we can repurpose any native Turing inference algorithm that
provides a marginal likelihood estimate into a target-aware inference
strategy.

To fully realize the potential of EPT, we further implement an annealed
importance sampling (AnIS) [31] inference engine for Turing. This allows for
effective marginal likelihood estimation in a much wider array of problems
than Turing’s previous support inference strategies. Using this, we show that
EPT can be used to express and run effective inference for a number of
problems, finding that it produces empirical estimates that are significantly
more accurate than conventional usage of Turing. We further confirm the
statistical soundness of EPT theoretically, showing that it produces
consistent estimates under nominal assumptions.

##  2 Background and Preliminaries

We are concerned with computing expectation of the form
$\mathbb{E}_{\pi(x)}[f(x)]$ where $f(x)$ is known, but $\pi(x)$ cannot be
directly evaluated or sampled from. Namely, $\pi(x)=\gamma(x)/Z$ where
$\gamma(x)$ is a known unnormalized distribution, but $Z$ is an unknown
normalization constant (also sometimes referred to as the marginal
likelihood). The standard approach to estimating $\mathbb{E}_{\pi(x)}[f(x)]$
is to first approximate $\pi(x)$ (e.g. with Monte Carlo samples) and then, in
turn, use this to approximate the expectation.

###  2.1 Target-Aware Bayesian Inference

Unfortunately, this standard pipeline ignores information about $f$ and is
therefore suboptimal if this information is available at inference time [14].
The recently proposed Target-Aware Bayesian Inference (TABI) framework of
Rainforth et al. [15] provides a means of getting around this by directly
creating an estimator for $\mathbb{E}_{\pi(x)}[f(x)]$ that incorporates
relevant information in $f$. It does this by breaking down the expectation
into three parts:

$\mathbb{E}_{\pi(x)}[f(x)]=(Z^{+}_{1}-Z^{-}_{1})/Z_{2}$ (1)

where

$\displaystyle Z^{+}_{1}=\int\gamma(x)f^{+}(x)dx,\quad
Z^{-}_{1}=\int\gamma(x)f^{-}(x)dx,\quad Z_{2}=\int\gamma(x)dx,$ $\displaystyle
f^{+}(x)=\max(f(x),0),\quad f^{-}(x)=-\min(f(x),0).$

To estimate the original expectation one separately estimates each term and
then recombines these estimates to an overall estimator as
$\mathbb{E}_{\pi(x)}[f(x)]\approx(\hat{Z}_{1}^{+}-\hat{Z}_{1}^{-})/\hat{Z}_{2}.$

Critically, each individual term can often be estimated more accurately in
isolation than the original expectation. To see this, first note that the
three subcomponents can be seen as the normalization constants of the three
densities $\gamma_{1}^{+}(x)\propto\gamma(x)f^{+}(x)$,
$\gamma_{1}^{-}(x)\propto\gamma(x)f^{-}(x)$, and $\gamma_{2}(x)=\gamma(x)$
respectively. The TABI framework now allows us to define a separate estimator
_tailored_ to each of these problems. For example, in an importance sampling
context we can use distinct proposals designed to be as accurate as possible
for each. By contrast, targeting the original expectation directly with a
single sampler is implicitly using the same proposal for each component and
this proposal cannot simultaneously be optimal for all of them. More
generally, the TABI framework allows us to repurpose any algorithm which
provides estimates of the normalization constant into a target-aware inference
algorithm separately applying it to each of $\gamma_{1}^{+}(x)$,
$\gamma_{1}^{-}(x)$, and $\gamma_{2}(x)$.

TABI’s most significant theoretical property is that it can achieve an
arbitrarily low error for any fixed sample budget ($\geq 3$). This is in stark
contrast to standard approaches such as self-normalized importance sampling
(SNIS) or Markov Chain Monte Carlo (MCMC); their expected error is lower
bounded, even when using an optimal proposal/sampler. Rainforth et al. [15]
further show that TABI methods can provide substantial empirical gains in
practice, often even outperforming theoretically optimal traditional
estimators, with gains increasing with the degree of mismatch between $\pi(x)$
and $\pi(x)f(x)$.

###  2.2 Turing Programs as Densities

⬇

@model function model(y)

x |$\sim$| Normal(0, 1)

@addlogprob!(0.1)

y |$\sim$| Normal(x, 1)

end

Figure 1: Turing example.

We now present a brief overview of the Turing PPL [19] (see
https://turing.ml/dev/docs/using-turing/ for full documentation). We will use
the simple Turing program in Figure 1 as a running example. We further provide
our own new precise formalism for the densities of Turing programs by
extending the approach of [32, §4.3], necessitated by some technical
intricacies of the expectation programming approach.

A Turing program is defined similarly to a normal Julia function [33]: the `
`macro indicates the definition of a Turing model which means we can use tilde
statements inside the function body, e.g. \eptx|$\sim$|Normal(0, 1), to denote
probabilistic model components. Observed data can be passed in as a formal
argument to the function. If the variable name on the left-hand side of the
tilde statement is not part of the arguments of the functions then it is
interpreted as a random variable. Let $x_{1:n}$ denote the set of direct
outputs from sampling statements and $y_{1:m}$ the observed data. We can view
Turing programs as defining an unnormalized density $\gamma(x_{1:n})$ (with an
implicit appropriate reference measure). For a given $x_{1:n}$, the density
gets computed by executing the program just like a normal Julia program, while
keeping track of the density of the current execution. More specifically, when
Turing reaches a tilde statement corresponding to a random variable it samples
a value for $x_{i}$, evaluates the density of this draw, and factors this into
the overall execution density. We denote the density of the draw as
$g_{i}(x_{i}|\eta_{i})$, where $g_{i}$ denotes the form of the sampling
statement and $\eta_{i}$ its parameters (which are a subset of the variables
currently in scope). For the tilde statements corresponding to the observed
data, it simply evaluates the density function $h_{j}(y_{j}|\phi_{j})$—where
$h_{j}$ and $\phi_{j}$ are analogous to $g_{i}$ and $\eta_{i}$
respectively—and factors the overall density accordingly. Sometimes a user
might want to add additional factors to the density without using a tilde
statement. For this, Turing provides the ` `primitive which multiplies the
density of the current execution by an arbitrary value \eptexp(log_p). We use
$\psi_{1},\dots,\psi_{K}$ to denote all the terms that are added to the
density using `.` Putting these together, the density defined by any valid
program trace can be written as \beginequation γ(x_1:n) = ∏_i=1^ng_i(x_i|η_i)
∏_j=1^mh_j(y_j|ϕ_j) ∏_k=1^K exp(ψ_k). Figure 1 thus defines the density
$\gamma(x)\\!=\\!\exp(0.1)\mathcal{N}\\!(x;0,1)\mathcal{N}\\!(y;x,1),$ with a
fixed input $y$. Note here that everything (i.e.
$n,x_{1:n},\eta_{1:n},g_{1:n},m,y_{1:m},\phi_{1:m},h_{1:m},K,\psi_{1:K}$) can
be a random variable because of potential stochasticity in the program path.
However, using the program itself, everything is deterministically calculable
from $x_{1:n}$, which can thus be thought of as the ‘raw’ random draws that
dictate all the randomness of the program; everything else is a pushforward of
these.

##  3 Expectation Programming

We now introduce the concept of an _expectation programming framework_ (EPF).
At a high level, this is a system that helps automate the estimation of
expectations in a target-aware manner. More precisely, we introduce the
concept of an expectation program as being analogous to a probabilistic
program, but where the quantity of interest is the expectation of its return
values under the program’s conditional distribution, rather than the
conditional distribution itself. We then define an EPF as a system which
provides automated inference engines that estimate this expectation in a
target-aware manner. In mathematical notation, if the unnormalized conditional
distribution of a program is $\gamma(x)$ and its return value is $f(x)$, then
the expectation program defines the expectation $\mathbb{E}_{\pi(x)}[f(x)]$
where $\pi(x)=\gamma(x)/Z$ is the normalized version of $\gamma(x)$. We will
assume that $f(x)$ is a scalar for now.

⬇

@expectation function expt_prog(y)

x |$\sim$| Normal(0, 1) # $x \sim\mathcal{N}(x; 0, 1)$

y |$\sim$| Normal(x, 1) # $y \sim\mathcal{N}(y; x, 1)$

return x^3 |$\phantom{\sim}$|# $f(x) = x^3$

end

expct_estimate, diagnostics = estimate_expectation(expt_prog(2), TABI(

marginal_likelihood_estimator=TuringAlgorithm(AnIS(), num_samples=1000)))

Figure 2:  An example of estimating an expectation with EPT. Here `` is our
“do inference” call which takes in expectation program `` (which itself is
given input $y=2$) and an inference method to apply (here a TABI estimator
using annealed importance sampling with $1000$ samples), and returns an
estimate for expected return value of ``. Appendix C shows how to estimate
this expectation using plain Turing in a non-target-aware manner.

###  3.1 Expectation Programming in Turing

Our particular realization of the EPF concept, _Expectation Programming in
Turing_ (EPT), allows users to specify $\gamma(x)$ as they normally would
using Turing’s ` `macro, and it interprets the \eptreturn semantics as
(implicitly) defining $f(x)$. We then use information about $f(x)$ contained
in the program itself to construct TABI–based estimators. Using the return
values provides a very convenient mechanism of defining $f(x)$ as it
automatically ensures all required input variables are in scope and allows
arbitrary functions to be defined by the deterministic components of the
program. The key component of the EPT implementation is to split up the
estimation of the desired expectation as per the TABI framework introduced in
Section 2.1. To do so we use source-code transformations to generate three
different Turing programs, one for each of the densities $\gamma_{1}^{+}(x)$,
$\gamma_{1}^{-}(x)$, and $\gamma_{2}(x)$. We then estimate the expectation by
individually estimating the normalization constant of each of these densities
and then combining them as per Equation (1). Generating valid Turing programs
has the advantage that we are able to leverage any inference algorithm in
Turing that provides marginal likelihood estimates to estimate the quantities
$Z^{+}_{1}$, $Z^{-}_{1}$, and $Z_{2}$. This modularity means that we do not
have to implement custom inference algorithms that would only work with EPT.
Estimating expectations with EPT is done in two stages. First, users define an
expectation program, say `,` with the \ept@expectation macro, which is a drop-
in replacement for `.` Using code transformations, \ept@expectation
automatically generates the three Turing programs representing the densities
$\gamma_{1}^{+}(x)$, $\gamma_{1}^{-}(x)$, and $\gamma_{2}(x)$. This happens
behind the scenes and the user does not need to deal with the transformed
programs directly. To estimate the expectation, the user then calls `,` where
\eptmethod specifies the estimation approach to be used. At present, the only
supported class of methods is `,` which implements the previously explained
TABI estimators, but the syntax is designed to allow for easy addition of
hypothetical other approaches (e.g.~control variates). EPT then estimates the
different normalization constants
($Z^+_1,Z^-_1$,and$Z_2$)byrunningaTuringinferencealgorithmoneachTuringprogramgeneratedby{\small\verb{}{@expectation}a{}{@expectation}ndcombiningthenormalizationconstantestimatestoformanestimateoftheexpectation.IntheexampleinFigure~\reffig:ep_{p}rogram}$,
we use ` `with annealed importance sampling \eptAnIS

, which is a new plain Turing inference algorithm that we have added to the
system for the purposes of this paper. ` `is a thin-wrapper object storing the
necessary information that allows \eptTABI to use a Turing inference method. `
`can be substituted with any other Turing inference algorithm that returns a
marginal likelihood estimate. Here \eptAnIS() implies the use of some
arbitrary default AnIS parameters regarding the Markov chain transition
kernel, and the number and spacing of intermediate potentials used.

###  3.2 Program Transformations

We now consider how to generate the Turing programs corresponding to each of
the TABI densities. Note that expectation programs in EPT are also valid
Turing models, i.e., replacing ` `with \ept@model yields a valid Turing
program. Such a program corresponds to the unnormalized density
$\gamma_{2}(x)\\!=\\!\gamma(x)$ without requiring any transformation. To
create a Turing program corresponding to $\gamma_{1}^{+}(x)$, we need to
multiply the unnormalized density of the unaltered Turing program $\gamma(x)$
by $\max(f(x),0)$. This is achieved using Turing’s aforementioned `
`primitive, such that we can think of it as adding a new
$\exp(\psi_K+1)=max(f(x_1:n),0)$factortotheprogramdensitydefinitionin~{}\eqref{eq:ppl_density}.Ourtransformationsarepatternmatchingproceduresthatfindallthe{\small\verb{}{return
expr} {}{return
expr}statementsinthefunctionbodyandthena)createanewlocalvariable\epttmp=expr}$
(where ` `is a unique identifier generated using Julia’s \eptgensym()), b)
insert a statement `b`efore the \eptreturn, and c) change the return statement
itself to `.` A concrete example of the transformation is presented in
Figure~\reflst:program_transformation. The transformation for
$\gamma_{1}^{-}(x)$ is analogous but inserts a statement ` `instead.
\beginfigure*[t]  ⬇ @expectation function expt_prog(y)  x |$\sim$| Normal(0,
1)  y |$\sim$| Normal(x, 1)  return x^3  end ⬇ @model function expt_prog(y)  x
|$\sim$| Normal(0, 1)  y |$\sim$| Normal(x, 1)  tmp = x^3
@addlogprob!(log(max(tmp, 0)))  return tmp end The results of one of the three
program transformations applied to the EPT `` program from Figure 3 [left].
Presented is the transformation into a valid Turing `` program [right]
corresponding to the density $\gamma_{1}^{+}(x)\propto\gamma(x)f^{+}(x)$. The
transformed code fragment is highlighted. The full transformation is slightly
more complex due to Turing’s internals. Appendix D shows the full source code
transformation for this model.  Users can define multiple expectations by
specifying multiple return values, while each individual return value needs to
almost surely be a numerical scalar. This ensures that each target expectation
is well defined and individually identified. For each return expression, we
apply our program transformation separately and derive a corresponding TABI
estimator for each. For example, if we have `,` the program transformation for
$\γ_1^+(x)}_2$wouldaddthestatement{\small\verb{}{@addlogprob!(log(max(expr2,
0)))}.{}{@addlogprob!(log(max(expr2,
0)))}Appendix~\refapd:multiple_{e}xpectations}$ shows a full example of how to
estimate multiple expectations with EPT.

###  3.3 Demonstration of Validity

We now demonstrate the statistical correctness of the EPT approach. For
simplicity, we will assume throughout that programs almost surely return a
single scalar value.111That is, the probability that the return value fails to
be a well-defined scalar (e.g. that it returns a char) is 0. Generalization to
programs with multiple return values is straightforward (provided the number
of return values is fixed) by considering each return value separately in
isolation (as EPT does itself). To start, we formally define probabilistic
programs and how they define expectations in an EPF.

######  Definition 1.

A probabilistic program $\mathcal{P}$ in a probabilistic programming language
defines an unnormalized density $\gamma(x_{1:n})$ over the _raw random draws_
$x_{1:n}\in\mathcal{X}$ of the program, which collectively we refer to as the
program trace,222Here it suffices to assume that all program traces have the
same length, but the semantic framework extends to programs with stochastic
support that generate traces of varying (even unbounded) lengths. along with
an implicitly defined reference measure $\mu$.
$\pi(x_{1:n})=\gamma(x_{1:n})/Z$ denotes the normalized density with
$Z=\int_{\mathcal{X}}\gamma(x_{1:n})d\mu(x_{1:n})$. $\pi(x_{1:n})$ and $\mu$
combined implicitly define the probability distribution specified by
$\mathcal{P}$, which we denote
$\mathbb{P}(A)=\int_{A}\pi(x_{1:n})d\mu(x_{1:n})$. Additionally, $\mathcal{P}$
defines a mapping from the raw random draws to the return values which we
denote $f:\mathcal{X}\rightarrow\mathbb{R}$.

Most PPLs, including Turing, focus solely on approximating $\pi(x_{1:n})$ and
ignore the definition of $f(x_{1:n})$. However, we can view probabilistic
programs as defining a random variable over return values, which we denote as
$F:=f(x_{1:n})$, and the expectation of this random variables is

$\mathbb{E}[F]=\int_{\mathcal{X}}f(x_{1:n})d\mathbb{P}(x_{1:n})=\int_{\mathcal{X}}f(x_{1:n})\pi(x_{1:n})d\mu(x_{1:n}),$
(2)

which is how they are interpreted under EPF. Note that the formal definition
of the function we are taking the expectation of is that it is the full
mapping from the raw random draws to the returned values rather than what is
lexically written in the ` `statement(s). This is why, for instance, it is
still valid to have multiple different \eptreturn statements in a program,
provided each ` `statement defines the same number of return values. In
practice, this is not something we need to worry about when writing either
models or inference engines as the law of the unconscious statistician
relieves us from explicitly delineating the random variable defined by our
function (the expectation of this random variable does not vary if we change
the parameterization of our model). However, the distinction is important for
ensuring validity and understanding the breakdown TABI uses
(i.e.~$\gamma^\pm$). To ensure that the induced probability measure of a
program is well-defined we generally assume that we have a valid unnormalized
density. In the case of Turing, this guarantees that there is a valid
probability distribution the inference algorithm can converge to. We formalize
this as follows. \begindefinition[Valid Turing Program] A Turing program
$\mathcal{P}$ is _valid_ (and defines a valid unnormalized probabilistic
program density $\gamma(x_{1:n})$) if all of the following hold:
$\gamma(x_{1:n})\geq 0,\forall x_{1:n}\in\mathcal{X}$;
$0<\int_{\mathcal{X}}\gamma(x_{1:n})d\mu(x_{1:n})<\infty$; and
$\gamma(x_{1:n})$ factorizes according to Equation (1). Programs in EPT are
largely equivalent to Turing programs but because they define expectations
they require a slightly different set of assumptions to ensure that a given
EPT program is valid.

######  Definition 2 (Valid EPT Program).

An EPT program $\mathcal{E}$ with density $\gamma(x_{1:n})$ on program traces
$x_{1:n}\in\mathcal{X}$ and return value $F=f(x_{1:n})$ is _valid_ if both of
the following hold: $\gamma(x_{1:n})$ is a valid unnormalized probabilistic
program density; and $F$ is integrable, i.e. $\mathbb{E}[|F|]<\infty$.

The assumptions in Definition 2 are actually strictly weaker than those made
implicitly by most existing PPLs (cf. Appendix B). Note that the integrability
of $F$ implies both $\mathbb{E}[F^{+}]$ and $\mathbb{E}[F^{-}]$ are finite.
This is important to establish the statistical correctness of the EPT approach
as follows.

######  Theorem 1.

Let $\mathcal{E}$ be a valid expectation program in EPT with unnormalized
density $\gamma(x_{1:n})$, defined on possible traces $x_{1:n}\in\mathcal{X}$,
and return value $F=f(x_{1:n})$. Then
$\gamma_{1}^{+}(x_{1:n}):=\gamma(x_{1:n})\max(0,f(x_{1:n}))$,
$\gamma_{1}^{-}(x_{1:n}):=-\gamma(x_{1:n})\min(0,f(x_{1:n}))$, and
$\gamma_{2}(x_{1:n}):=\gamma(x_{1:n})$ are all valid unnormalized
probabilistic program densities. Further, if $\\{\hat{Z}_{1}^{+}\\}_{n}$,
$\\{\hat{Z}_{1}^{-}\\}_{n}$, $\\{\hat{Z}_{2}\\}_{n}$ are sequences of
estimators for $n=\mathbb{N}^{+}$ such that

$\\{\hat{Z}_{1}^{\pm}\\}_{n}\overset{p}{\to}\int_{\mathcal{X}}\gamma^{\pm}_{1}(x_{1:n})d\mu(x_{1:n}),\quad\\{\hat{Z}_{2}\\}_{n}\overset{p}{\to}\int_{\mathcal{X}}\gamma_{2}(x_{1:n})d\mu(x_{1:n})$

where $\overset{p}{\to}$ means convergence in probability as $n\to\infty$,
then
$(\\{\hat{Z}_{1}^{+}\\}_{n}-\\{\hat{Z}_{1}^{-}\\}_{n})/\\{\hat{Z}_{2}\\}_{n}\overset{p}{\to}\mathbb{E}[F].$

Theorem 1, which is proved in Appendix B, shows that if we have programs with
the desired densities and we use consistent marginal likelihood estimators for
each, then our resulting expectation estimates will themselves be consistent.
The latter requirement is covered by the consistency of Turing’s own inference
engines. The former requires that our transformed programs are valid Turing
programs with the intended densities. This is indeed the case as we now show.
Given an input EPT program $\mathcal{E}$ we apply transformations to get the
three Turing programs $\mathcal{P}_{1}^{+}$, $\mathcal{P}_{1}^{-}$, and
$\mathcal{P}_{2}$ with $\gamma_{1}^{+}(x_{1:n})$, $\gamma_{1}^{-}(x_{1:n})$,
and $\gamma_{2}(x_{1:n})$ as their respective densities. To ensure that the
transformations for $\gamma_{1}^{+}(x_{1:n})$ and $\gamma_{1}^{-}(x_{1:n})$
are correct, we need to ensure that a) the inserted code in our
transformations is itself valid, b) the transformation does not have any
unintended side effects, and c) the new density terms add valid factors to the
program density. The first is true as the operation of the transformed
sections of code are identical to the originals except for the new ` `terms,
which themselves produce no outputs and, by construction, use only the
variables that are in scope. The second is guaranteed by ensuring that the
\epttmp variables are given unique identifiers that cannot clash with each
other or any other variables in the program. The third follows from the
restriction that each return value must almost surely be a numerical scalar,
coupled with the fact that the added density factors (namely ` `and \ept-
min(tmp, 0)) are non-negative by construction. Thus, we see that EPT will
produce a consistent estimation of program expectations, under the assumptions
of Definition 2 and the consistency of the Turing base inference algorithms.

##  4 Related Work

Though a few papers [2, 34] have formalized the expectation defined by a
probabilistic program, none does this from the perspective of directly
targeting this expectation as the quantity to estimate. As a pedagogical
example, Gordon et al. [2] introduce a small PPL which exemplifies how return
statements can be interpreted as defining expectations. Their main goal,
however, is still posterior inference. Hakaru [29] provides the primitive
`expect(m, f)} \citep`zinkov2017Composinga which takes as input a program
defining a measure and a function and returns a program that denotes the
integral of the function w.r.t. the measure. This primitive is then used to
construct algorithms for posterior inference. The end goal of Hakaru is to
provide a modular interface for inference algorithms by computing conditional
distributions analytically. In contrast, our framework has the goal of
computing _expectations_ as accurately as possible via automated target-aware
inference. Some PPLs provide syntactic sugars for forming expectation
estimates from the samples produced by inference, but these do not adjust the
inference itself to exploit target function information. For example, in Stan
[17] users can apply target functions to posterior samples using the
generated_quantities block. Similarly, in Pyro [16] the return values are
stored along with MCMC posterior samples, thus allowing expectations to be
estimated by taking empirical averages. PyMC3 [20] allows users to track
deterministic transformations of the latent variables. Turing itself also
provides a generated_quantities function, similar to Stan (see Appendix C for
an example).

Figure 3: Relative squared error and effective sample size for the Gaussian
posterior predictive experiment. The solid lines show the median of the
estimator while the shaded region show the 25 % and 75 % quantiles. Medians
and quantiles are computed over 10 separate runs with different random seed
for the posterior predictive problem. For the ESS plot we are plotting
$\min(\text{ESS}_{Z_{1}},\text{ESS}_{Z_{2}})$ for TAAnIS and
$\min(\text{ESS}_{\text{AnIS retargeted}},\text{ESS}_{\text{AnIS}})$ for AnIS.

##  5 Experiments

We demonstrate the effectiveness of the EPT target-aware inference methods on
three problems: a synthetic numerical example, an SIR epidemiology model, and
a Bayesian hierarchical model of Radon concentration. The code for all
experiments can be found at https://git.io/JZOqN. The performance of EPT
depends on the performance of the chosen marginal likelihood estimator. At the
time of writing, Turing provides implementations (only with the prior as the
proposal) of Sequential Monte Carlo [35] and Importance Sampling (IS) as
inference algorithms that provide marginal likelihood estimates. We
implemented a new Turing inference engine that uses Annealed Importance
Sampling (AnIS) [31] (see Appendix A). We chose AnIS because it has proved to
be effective at estimating the normalization constant even for high-
dimensional models [36, 37, 38]. AnIS requires setting two hyperparameters: an
annealing schedule and a transition kernel. Currently, users can choose
between two transition kernels: Metropolis-Hastings (MH) implemented in
AdvancedMH.jl [39] and Hamiltonian Monte Carlo (HMC) [40, 41, 42] in
AdvancedHMC.jl [43]. We will use AnIS and TAAnIS to refer to the non-target-
aware and the target-aware AnIS, respectively. To ensure a fair comparison we
use the same setup and hyperparameters for each estimator. We also compare
directly to MCMC targeting the posterior and using the same type of transition
kernel as AnIS and TAAnIS. This transition kernel is MH in Section 5.1 and HMC
elsewhere. Detailed configurations for each experiment are given in Appendix
F. To compare the performance of the estimators we look at the effective
sample size (ESS) and the relative squared error (RSE)
$\hat{\delta}:=(\hat{\mu}-\mu)^{2}/\mu^{2}$ , where $\mu$ denotes the ground-
truth value and $\hat{\mu}$ is the estimate. All our experiments correspond to
target functions which are always positive, so henceforth we use $Z_{1}$ to
refer to $Z^{+}_{1}$ as we always have $Z^{-}_{1}=0$. Appendix E shows how to
configure EPT to avoid the computation for $Z^{-}_{1}$ when this is the case.
Both TAAnIS and AnIS produce unnormalized weighted samples $w_{i}$ so the ESS
is calculated as $(\sum_{i}w_{i})^{2}/\sum_{i}w_{i}^{2}$. For MCMC we omit the
ESS comparison and focus solely on the RSE, as MCMC ESS estimates, which are
based on autocorrelation [44], are not directly comparable with those based on
importance weights, and unreliable if the RSE is large. For TAAnIS we can look
at the ESS of the two AnIS runs which we denote $\text{ESS}_{Z_{1}}$ and
$\text{ESS}_{Z_{2}}$. Similarly, $\text{ESS}_{\text{AnIS}}$ is the ESS for
AnIS. We use $\text{ESS}_{\text{AnIS retargeted}}$ to denote what the ESS
would be if we use the AnIS samples to estimate $Z_{1}$ by multiplying the
AnIS importance weights by $f(x)$. When comparing our estimators based on RSE
and ESS we ensure that we give each estimator the same computational budget
i.e. same number of likelihood evaluations. Thus if we have an AnIS estimator
with $n$ intermediate distributions and $K$ samples, then we compare it to the
TAAnIS estimator which uses $n$ intermediate distributions and $K/2$ samples
for each of the two separate terms. In turn, we compare AnIS and TAAnIS to an
MCMC estimator with $n\cdot K$ samples.

###  5.1 Gaussian Posterior Predictive

The first problem considered is calculating the posterior predictive
distribution of a Gaussian model with an unknown mean, where
$\gamma(\mathbf{x})=\mathcal{N}(\mathbf{x};0,I)\,\mathcal{N}(\mathbf{y};\mathbf{x},I)$
and $f(\mathbf{x})=\mathcal{N}(-\mathbf{y};\mathbf{x},\frac{1}{2}I)$ are the
unnormalised density and target function, respectively. We assume our observed
data is $\mathbf{y}=(3.5/\sqrt{10})\mathbf{1}$ where $\mathbf{1}$ is a
10-dimensional vector of ones. Using EPT we can express this expectation in
just 5 lines of code—the full model is given in Appendix J. This problem is
amenable to an analytic solution so allows us to compute the error of the
estimates. Figure 3 compares the performance of TAAnIS, AnIS, and MCMC (here
MH). We see a clear benefit to using the target-aware inference algorithm to
estimate the expectation. TAAnIS achieves a lower RSE, and the ESS highlights
the advantage of using separate marginal likelihood estimators for $Z_{1}$ and
$Z_{2}$.

###  5.2 SIR Epidemiological Model

Figure 4: Relative squared error and effective sample size for the SIR
experiment. Conventions as in Figure 3. Results are computed over 5 runs for
different random seeds.

Our second problem setting is a more applied example based on the Susceptible-
Infected-Recovered (SIR) model of [45] from the field of epidemiology. Assume
we face an outbreak of a contagious disease. The government has provided us
with a function yielding the expected cost of the disease which depends on the
basic reproduction rate $R_{0}$, which indicates the expected number of people
one infected person will infect in a population where everyone is susceptible.
We seek to infer $R_{0}$ and the expected cost of the outbreak. The SIR model
divides the population into three compartments: people who are susceptible to
the disease, those who are currently infected, and those who have already
recovered. The dynamics of the outbreak are modelled by a set of differential
equations with parameters $\beta$ and $\gamma$:

$\displaystyle\frac{dS}{dt}=-\beta S\frac{I}{N},\quad\frac{dI}{dt}=\beta
S\frac{I}{N}-\gamma I,\quad\frac{dR}{dt}=\gamma I,$ (3)

where $S$, $I$ and $R$ correspond to the number of people susceptible,
infected and recovered, respectively. The size of the total population is
$N=S+I+R$. Roughly, $\beta$ models the constant rate of infectious contact
between people, while $\gamma$ is the constant recovery rate of infected
individuals. From these parameters we can calculate the basic reproduction
rate $R_{0}=\beta/\gamma$. We assume $\gamma$ to be known, and we want to
infer $\beta$ and the initial number of infected people $I_{0}$. The full
statistical model along with the cost function which is defined with respect
to $R_{0}$ is given in Appendix G. This scenario is a good use case for EPT
because we are interested in estimating a specific expectation with high
accuracy. Furthermore, our cost function has some outcomes which might have
low probability under the posterior but which incur a very high cost. These
are liable to be missed by non-target-aware schemes, leading to skew
estimators that tend to underestimate the expectation. Figure 4 compares the
performance of the estimators. Since this problem is not amenable to an
analytic solution, we estimate the ground-truth using a customized IS
estimator with orders of magnitude more samples than estimates presented in
the plot (see Appendix F). We see that our approach substantially improves on
the baselines, with MCMC (here HMC) failing to provide any meaningful
estimate; it produces no samples where $f(x)$ is significant. EPT with TAAnIS
overcomes this through its use of a separate estimator for $\gamma(x)f(x)$.
The fact that MCMC does far worse than AnIS, despite neither being target-
aware, stems from the latter producing a greater diversity of (weighted)
samples, a small number of which land in regions of high $f(x)$ by chance (cf.
Appendix L).

###  5.3 Hierarchical Model of Radon Concentration

Our third problem setting is a Bayesian hierarchical model for the radon
concentration in households in different counties, adapted from Gelman and
Hill [46]. In this, we observe data from $N$ counties which are indexed by $i$
and there are $M_{i}$ houses in the $i^{\text{th}}$ county. For the
$j^{\text{th}}$ house in county $i$, we would like to predict the log radon
concentration $y_{ij}$ inside the house. For each house we have a covariate
$x_{ij}$ which is $0$ if the house has a basement, and $1$ if it does not.
With this setup, the model is defined as

$\displaystyle\mu_{\alpha}$ $\displaystyle\sim\mathcal{N}(0,10),$
$\displaystyle\alpha_{i}$
$\displaystyle\sim\mathcal{N}(\mu_{\alpha},\sigma_{\alpha}),$
$\displaystyle\mu_{\beta}$ $\displaystyle\sim\mathcal{N}(0,10),$ (4a)
$\displaystyle\beta_{i}$
$\displaystyle\sim\mathcal{N}(\mu_{\beta},\sigma_{\beta}),$
$\displaystyle\epsilon$ $\displaystyle\sim\text{HalfCauchy}(0,5),$
$\displaystyle y_{ij}$
$\displaystyle\sim\mathcal{N}(\alpha_{i}+\beta_{i}x_{ij},\epsilon).$ (4b)

We set $\sigma_{\alpha}=0.12$ and $\sigma_{\beta}=0.22$, with these values
informed by preliminary analysis of the data.

Table 1: Final estimates for the Radon experiments. The mean and standard deviation are estimated over 10 runs. Method | Final Estimate
---|---
TAAnIS | $3.74\mathrm{e}{-8}\,\pm\,2.39\mathrm{e}{-9}$
AnIS | $1.15\mathrm{e}{-9}\,\pm\,3.02\mathrm{e}{-9}$
MCMC | $7.79\mathrm{e}{-18}\,\pm\,2.46\mathrm{e}{-17}$

We are interested in a decision problem similar to the one presented in [47,
Chapter 9]. Concretely, we want to find out whether the radon level in all
households is below an acceptable level. Informed by Gelman et al. [47] we
assume that a radon concentration below 4pCi/L does not require an
intervention. The probability of this event is equal to the expectation under
the posterior of a step function $f(x)$. However, to allow the use of HMC
transition kernels we use a logistic function as a continuous relaxation of
this step function. See Appendix H for details.

Figure 5: ESS plots for the Radon experiment. Conventions as in Figure 3.
Results are computed over 10 runs/seeds.

This problem cannot be solved analytically and, unlike before, estimating the
ground-truth with sufficient accuracy is computationally infeasible. We,
therefore, resort to comparing TAAnIS (i.e. EPT) and AnIS based on their ESSs,
noting that a low ESS almost exclusively means a poor inference estimate,
while a high ESS is a strong (but not absolute) indicator of good performance.
As we can see in Figure 5, EPT outperforms standard AnIS by several orders of
magnitude. Additionally, Table 1 presents the final expectation estimates for
each method. All methods differ in their estimates. However, TAAnIS is the
only one where the standard deviation of the estimate is small relative to its
mean estimate, which, coupled with our ESS results, provides strong evidence
that it is significantly outperforming the baselines. In particular, it seems
clear that the MCMC (here HMC) estimate is very poor, as its estimate is many
orders of magnitude smaller than the others, indicating it has failed to
produce any samples in regions where $f(x)$ is non-negligible.

##  6 Discussion

Limitations. The efficiency and ease of use of EPT depend largely on the
underlying marginal likelihood estimator. AnIS is a powerful algorithm but it
has several hyperparameters that need to be adjusted for a given model. This
is not a fundamental limitation of EPT, because any new marginal likelihood
estimators that emerge, which are more efficient and have fewer free
parameters, can easily be incorporated into EPT. However, marginal likelihood
estimation can itself sometimes be a more challenging task than generating
posterior samples, as MCMC methods are suitable for the latter but not the
former. As mentioned in Section 2.1, EPT inherits the property of TABI that
its performance gain depends on the degree of mismatch between $\pi(x)$ and
$\pi(x)f(x)$. Consequently, its use of TABI-based estimators may not always be
beneficial for problems where this mismatch is low. Conclusions. We have
introduced the concept of an expectation programming framework (EPF). An EPF
allows users to encode expectations through the return values of a
probabilistic program and then employs target-aware inference strategies to
automatically estimate this effectively. We have presented a specific
realization of an EPF, called EPT, based on Turing and shown that it
accurately estimates expectations in practice, through a combination of
program transformations and TABI-based estimators. The modularity of EPT means
that it can easily be built on by others through alternative base inference
algorithms or target-aware strategies. All in all, we believe the general
concept of EPF is an exciting addition to the literature that paves the way
for further advancements.

## Acknowledgments and Disclosure of Funding

We would like to thank Sheheryar Zaidi for helpful discussions on configuring
Annealed Importance Sampling. Tim Reichelt and Adam Goliński are supported by
the UK EPSRC CDT in Autonomous Intelligent Machines and Systems.

## References

  * Robert and Casella [2004] Christian Robert and George Casella.  _Monte Carlo Statistical Methods_.  Springer Texts in Statistics. Springer-Verlag, New York, second edition, 2004.  ISBN 978-0-387-21239-5.  doi: 10.1007/978-1-4757-4145-2.
  * Gordon et al. [2014] Andrew D. Gordon, Thomas A. Henzinger, Aditya V. Nori, and Sriram K. Rajamani.  Probabilistic Programming.  In _Future of Software Engineering Proceedings_ , FOSE 2014, pages 167–181, New York, NY, USA, May 2014. Association for Computing Machinery.  ISBN 978-1-4503-2865-4.  doi: 10.1145/2593882.2593900.
  * van de Meent et al. [2018] Jan-Willem van de Meent, Brooks Paige, Hongseok Yang, and Frank Wood.  An Introduction to Probabilistic Programming.  _arXiv:1809.10756 [cs, stat]_ , September 2018.
  * Torrie and Valleau [1977] G.M. Torrie and J.P. Valleau.  Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling.  _Journal of Computational Physics_ , 23(2):187 – 199, 1977.  ISSN 0021-9991.  doi: https://doi.org/10.1016/0021-9991(77)90121-8.  URL http://www.sciencedirect.com/science/article/pii/0021999177901218.
  * Hesterberg [1988] Timothy Classen Hesterberg.  _Advances in Importance Sampling_.  PhD thesis, Stanford University, 1988.
  * Wolpert [1991] Robert L Wolpert.  Monte Carlo Integration in Bayesian Statistical Analysis.  _Contemporary Mathematics_ , 115:101–116, 1991.
  * Oh and Berger [1992] Man-Suk Oh and James O Berger.  Adaptive Importance Sampling in Monte Carlo Integration.  _Journal of Statistical Computation and Simulation_ , 41(3-4):143–168, 1992.
  * Evans et al. [1995] Michael Evans, Tim Swartz, et al.  Methods for Approximating Integrals in Statistics with Special Emphasis on Bayesian Integration Problems.  _Statistical science_ , 10(3):254–272, 1995.
  * Meng and Wong [1996] Xiao-Li Meng and Wing Hung Wong.  Simulating Ratios of Normalizing Constants via a Simple Identity: a Theoretical Exploration.  _Statistica Sinica_ , pages 831–860, 1996.
  * Chen et al. [1997] Ming-Hui Chen, Qi-Man Shao, et al.  On Monte Carlo Methods for Estimating Ratios of Normalizing Constants.  _The Annals of Statistics_ , 25(4):1563–1594, 1997.
  * Gelman and Meng [1998] Andrew Gelman and Xiao-Li Meng.  Simulating Normalizing Constants: From Importance Sampling to Bridge Sampling to Path Sampling.  _Statistical science_ , pages 163–185, 1998.
  * Lacoste-Julien et al. [2011] Simon Lacoste-Julien, Ferenc Huszár, and Zoubin Ghahramani.  Approximate Inference for the Loss-Calibrated Bayesian.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 416–424, 2011.
  * Owen [2013] Art B. Owen.  _Monte Carlo Theory, Methods and Examples_.  2013\.
  * Golinski et al. [2019] Adam Golinski, Frank Wood, and Tom Rainforth.  Amortized Monte Carlo Integration.  In _International Conference on Machine Learning (ICML)_ , pages 2309–2318. PMLR, May 2019.
  * Rainforth et al. [2020] Tom Rainforth, Adam Golinski, Frank Wood, and Sheheryar Zaidi.  Target–Aware Bayesian Inference: How to Beat Optimal Conventional Estimators.  _Journal of Machine Learning Research_ , 21(88):1–54, 2020.
  * Bingham et al. [2019] Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman.  Pyro: Deep Universal Probabilistic Programming.  _Journal of Machine Learning Research_ , 20(28):1–6, 2019.  ISSN 1533-7928.
  * Carpenter et al. [2017] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: A probabilistic programming language.  _Journal of statistical software_ , 76(1), 2017.
  * Cusumano-Towner et al. [2019] Marco F. Cusumano-Towner, Feras A. Saad, Alexander K. Lew, and Vikash K. Mansinghka.  Gen: A General-Purpose Probabilistic Programming System with Programmable Inference.  In _Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation_ , PLDI 2019, pages 221–236, New York, NY, USA, June 2019. Association for Computing Machinery.  ISBN 978-1-4503-6712-7.  doi: 10.1145/3314221.3314642.
  * Ge et al. [2018] Hong Ge, Kai Xu, and Zoubin Ghahramani.  Turing: A Language for Flexible Probabilistic Inference.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 1682–1690. PMLR, March 2018.
  * Salvatier et al. [2016] John Salvatier, Thomas V Wiecki, and Christopher Fonnesbeck.  Probabilistic Programming in Python Using PyMC3.  _PeerJ Computer Science_ , 2:e55, 2016.
  * Tran et al. [2016] Dustin Tran, Alp Kucukelbir, Adji B Dieng, Maja Rudolph, Dawen Liang, and David M Blei.  Edward: A Library for Probabilistic Modeling, Inference, and Criticism.  _arXiv preprint arXiv:1610.09787_ , 2016.
  * Wood et al. [2014] Frank Wood, Jan Willem Meent, and Vikash Mansinghka.  A New Approach to Probabilistic Programming Inference.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 1024–1032. PMLR, April 2014.
  * Mansinghka et al. [2014] Vikash K. Mansinghka, Daniel Selsam, and Yura N. Perov.  Venture: a higher-order probabilistic programming platform with programmable inference.  _CoRR_ , abs/1404.0099, 2014.  URL http://arxiv.org/abs/1404.0099.
  * Goodman et al. [2008] Noah D. Goodman, Vikash K. Mansinghka, Daniel Roy, Keith Bonawitz, and Joshua B. Tenenbaum.  Church: A language for generative models.  In _Proceedings of the Twenty-Fourth Conference on Uncertainty in Artificial Intelligence_ , UAI’08, page 220–229, Arlington, Virginia, USA, 2008. AUAI Press.  ISBN 0974903949.
  * Goodman and Stuhlmüller [2014] Noah D Goodman and Andreas Stuhlmüller.  The Design and Implementation of Probabilistic Programming Languages.  http://dippl.org, 2014.  Accessed: 2021-5-21.
  * Tavares et al. [2019] Zenna Tavares, Xin Zhang, James Koppel, Ria Das, and Armando Solar Lezama.  A language for counterfactual generative models.  2019\.
  * Fierens et al. [2015] Daan Fierens, Guy Van den Broeck, Joris Renkens, Dimitar Shterionov, Bernd Gutmann, Ingo Thon, Gerda Janssens, and Luc De Raedt.  Inference and learning in probabilistic logic programs using weighted boolean formulas.  _Theory and Practice of Logic Programming_ , 15(3):358–401, 2015.  doi: 10.1017/S1471068414000076.
  * Murray and Schön [2018] Lawrence M. Murray and Thomas B. Schön.  Automated learning with a probabilistic programming language: Birch.  _Annual Reviews in Control_ , 46:29–43, 2018.  ISSN 1367-5788.  doi: https://doi.org/10.1016/j.arcontrol.2018.10.013.  URL https://www.sciencedirect.com/science/article/pii/S1367578818301202.
  * Narayanan et al. [2016] Praveen Narayanan, Jacques Carette, Wren Romano, Chung-Chieh Shan, and Robert Zinkov.  Probabilistic Inference by Program Transformation in Hakaru (System Description).  In Oleg Kiselyov and Andy King, editors, _Functional and Logic Programming_ , volume 9613, pages 62–79. Springer International Publishing, Cham, 2016.  ISBN 978-3-319-29603-6 978-3-319-29604-3.  doi: 10.1007/978-3-319-29604-3_5.
  * Minka et al. [2018] Tom Minka, John M. Winn, John P. Guiver, Yordan Zaykov, Dany Fabian, and John Bronskill.  /Infer.NET 0.3, 2018.  Microsoft Research Cambridge. http://dotnet.github.io/infer.
  * Neal [2001] Radford M. Neal.  Annealed Importance Sampling.  _Statistics and Computing_ , 11(2):125–139, April 2001.  ISSN 0960-3174.  doi: 10.1023/A:1008923215028.  URL https://doi.org/10.1023/A:1008923215028.
  * Rainforth [2017] Tom Rainforth.  _Automating Inference, Learning, and Design Using Probabilistic Programming_.  http://purl.org/dc/dcmitype/Text, University of Oxford, 2017.
  * Bezanson et al. [2017] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah.  Julia: A Fresh Approach to Numerical Computing.  _SIAM Review_ , 59(1):65–98, January 2017.  ISSN 0036-1445.  doi: 10.1137/141000671.
  * Zinkov and Shan [2017] Robert Zinkov and Chung-Chieh Shan.  Composing Inference Algorithms as Program Transformations.  _Proceedings of Uncertainty in Artificial Intelligence (UAI)_ , page 10, 2017.
  * Del Moral et al. [2006] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra.  Sequential Monte Carlo Samplers.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(3):411–436, 2006.
  * Wallach et al. [2009] Hanna M. Wallach, Iain Murray, Ruslan Salakhutdinov, and David Mimno.  _Evaluation Methods for Topic Models_ , page 1105–1112.  Association for Computing Machinery, New York, NY, USA, 2009.  ISBN 9781605585161.  URL https://doi.org/10.1145/1553374.1553515.
  * Salakhutdinov and Larochelle [2010] Ruslan Salakhutdinov and Hugo Larochelle.  Efficient Learning of Deep Boltzmann Machines.  In Yee Whye Teh and Mike Titterington, editors, _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , volume 9 of _Proceedings of Machine Learning Research_ , pages 693–700, Chia Laguna Resort, Sardinia, Italy, 13–15 May 2010. JMLR Workshop and Conference Proceedings.  URL http://proceedings.mlr.press/v9/salakhutdinov10a.html.
  * Wu et al. [2017] Yuhuai Wu, Yuri Burda, R. Salakhutdinov, and Roger B. Grosse.  On the Quantitative Analysis of Decoder-Based Generative Models.  _International Conference on Learning Representations (ICLR)_ , 2017\.
  * Turing Development Team [2020] The Turing Development Team.  TuringLang/AdvancedMH.jl.  The Turing Language, October 2020.  URL https://github.com/TuringLang/AdvancedMH.jl.
  * Neal [2011] Radford M. Neal.  MCMC Using Hamiltonian Dynamics.  In _Handbook of Markov Chain Monte Carlo_ , pages 113–162. Chapman & Hall / CRC Press, 2011.
  * Hoffman and Gelman [2014] M. Hoffman and A. Gelman.  The No-U-turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15:1593–1623, 2014\.
  * Betancourt [2018] Michael Betancourt.  A Conceptual Introduction to Hamiltonian Monte Carlo, 2018\.
  * Xu et al. [2020] Kai Xu, Hong Ge, Will Tebbutt, Mohamed Tarek, Martin Trapp, and Zoubin Ghahramani.  AdvancedHMC.jl: A Robust, Modular and Efficient Implementation of Advanced HMC Algorithms.  In _Symposium on Advances in Approximate Bayesian Inference (AABI)_ , pages 1–10. PMLR, February 2020.
  * Vehtari et al. [2020] Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner.  Rank-Normalization, Folding, and Localization: An Improved $\widehat{R}$ for Assessing Convergence of MCMC.  _Bayesian Analysis_ , 2020.  doi: 10.1214/20-BA1221.  URL https://doi.org/10.1214/20-BA1221.  Advance publication.
  * Kermack et al. [1927] William Ogilvy Kermack, A. G. McKendrick, and Gilbert Thomas Walker.  A Contribution to the Mathematical Theory of Epidemics.  _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , 115(772):700–721, August 1927.  doi: 10.1098/rspa.1927.0118.
  * Gelman and Hill [2006] Andrew Gelman and Jennifer Hill.  _Data Analysis Using Regression and Multilevel/Hierarchical Models_.  Analytical Methods for Social Research. Cambridge University Press, 2006\.  doi: 10.1017/CBO9780511790942.
  * Gelman et al. [2013] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin.  _Bayesian Data Analysis_.  CRC press, 2013.
  * Grinsztajn et al. [2020] Leo Grinsztajn, Elizaveta Semenova, Charles C. Margossian, and Julien Riou.  Bayesian Workflow for Disease Transmission Modeling in Stan, 2020.

##  Appendix A Annealed Importance Sampling

Annealed importance sampling (AnIS) [31] is an inference algorithm which was
developed with the goal of efficiently estimating the normalization constant
$Z$ of an unnormalized density $\gamma(x)$. It works by defining a sequence of
annealing distributions $\pi_{0}(x),\dots,\pi_{n}(x)$ which interpolate
between a simple base distribution $\pi_{0}(x)$ (typically the prior for a
Bayesian model) and the complex target density $\pi_{n}(x)=\gamma(x)$. The
most common scheme is to take

$\pi_{i}(x)\propto\lambda_{i}(x)=\pi_{0}(x)^{1-\beta_{n}}\gamma(x)^{\beta_{n}},$
(5)

with $0\\!=\\!\beta_{0}\\!<\\!\dots\\!<\\!\beta_{n}\\!=\\!1$. The algorithm
further requires the definition of Markov chain transition kernels
$\tau_{1}(x,x^{\prime}),\dots,\tau_{n-1}(x,x^{\prime})$ and proceed to
generate the $j^{\text{th}}$ weighted sample as follows First, sample initial
particle $x^{(1)}_{j}\sim\pi_{0}(x)$, then for $i=1,\ldots,(n-1)$, generate
$x^{(i+1)}_{j}\sim\tau_{i}(x^{(i)}_{j},\cdot)$ and, finally, return sample
$x^{(n)}_{j}$ with weight

$w_{j}=\frac{\lambda_{1}(x_{j}^{(1)})\lambda_{2}(x_{j}^{(2)})\dots\lambda_{n}(x_{j}^{(n)})}{\pi_{0}(x_{j}^{(1)})\lambda_{1}(x_{j}^{(2)})\dots\lambda_{n-1}(x_{j}^{(n)})}$
(6)

We can estimate expectations with the weights and samples just as in
importance sampling. Thus we can estimate the expectation and the
normalization constant as

$\displaystyle\begin{split}\mathbb{E}_{\pi(x)}[f(x)]\approx\frac{\sum_{j=1}^{N}w_{j}f(x^{(n)}_{j})}{\sum_{j=1}^{N}w_{j}}\quad\text{and}\quad
Z\approx\frac{1}{N}\sum_{j=1}^{N}w_{j}.\end{split}$

##  Appendix B Theoretical Details

###  B.1 Assumptions in Definition 2

To ensure correctness most PPLs assume that a particular inference algorithm
will converge to the distribution of $F$ (i.e. the distribution over return
values). A standard PPL Monte Carlo inference engine will now produce a
sequence of samples $F_{n},~{}n=1,2,\dots$ and consistency requires that
$F_{n}$ converges in distribution to $F$ as $n\to\infty$. This is equivalent
to requiring that for _any_ integrable function $h$,
$\mathbb{E}[h(F_{n})]\to\mathbb{E}[h(F)]$; and it presupposes that the
distribution of $F$ is a finite measure, i.e., $\mathbb{E}[F]$ is finite. We
thus see our assumption is strictly weaker than that of standard PPLs that
allow return values from programs: we only need convergence in the case where
$h$ is the identity mapping, not all integrable functions.

###  B.2 Proof for Theorem 1

See 1

###### Proof.

We start by noting that as $\gamma_{2}(x_{1:n})$ is identical to
$\gamma(x_{1:n})$, it is by assumption a valid unnormalized program density.
Meanwhile, by construction,
$\gamma(x_{1:n})^{+}_{1},\gamma(x_{1:n})^{-}_{1}\geq 0,\forall
x_{1:n}\in\mathcal{X}$. Further, each can be written in the form of (1) by
taking the correspond definition of $\gamma(x_{1:n})$ and adding in factors
$\exp(\psi_{K+1})=\max(0,f(x_{1:n}))$ and
$\exp(\psi_{K+1})=-\min(0,f(x_{1:n}))$ for $\gamma(x_{1:n})^{+}_{1}$ and
$\gamma(x_{1:n})^{-}_{1}$ respectively. To finish the proof that
$\gamma^{\pm}(x_{1:n})$ are valid densities, we show that
$0<Z_{1}^{\pm}<\infty$. Starting with the standard definition of an
expectation for arbitrary random variables, we can express $\mathbb{E}[F]$ as

$\displaystyle\int_{\mathcal{X}}f(x_{1:n})d\mathbb{P}(x_{1:n})$
$\displaystyle=\int_{\mathcal{X}}f^{+}(x_{1:n})d\mathbb{P}(x_{1:n})-\int_{\mathcal{X}}f^{-}(x_{1:n})d\mathbb{P}(x_{1:n}).$
(7) Noting that if $F$ is integrable then by the definition of the Lebesgue
integral $\int_{\mathcal{X}}f^{+}(x_{1:n})d\mathbb{P}(x_{1:n})<\infty$ and
$\int_{\mathcal{X}}f^{-}(x_{1:n})d\mathbb{P}(x_{1:n})<\infty$. Now inserting
the distribution the program defines over $x_{1:n}$,
$\displaystyle=\int_{\mathcal{X}}f^{+}(x_{1:n})\pi(x_{1:n})d\mu(x_{1:n})-\int_{\mathcal{X}}f^{-}(x_{1:n})\pi(x_{1:n})d\mu(x_{1:n})$
(8) and noting that $\gamma(x_{1:n})\geq 0$ for all $x_{1:n}\in\mathcal{X}$
and $0<\int_{\mathcal{X}}\gamma(x_{1:n})d\mu(x_{1:n})<\infty$,
$\displaystyle=\frac{\int_{\mathcal{X}}f^{+}(x_{1:n})\gamma(x_{1:n})d\mu(x_{1:n})-\int_{\mathcal{X}}f^{-}(x_{1:n})\gamma(x_{1:n})d\mu(x_{1:n})}{\int_{\mathcal{X}}\gamma(x_{1:n})d\mu(x_{1:n})}$
(9)
$\displaystyle=\frac{\int_{\mathcal{X}}\gamma^{+}_{1}(x_{1:n})d\mu(x_{1:n})-\int_{\mathcal{X}}\gamma_{1}^{-}(x_{1:n})d\mu(x_{1:n})}{\int_{\mathcal{X}}\gamma_{2}(x_{1:n})d\mu(x_{1:n})}=:\frac{Z^{+}_{1}-Z^{-}_{1}}{Z_{2}}.$
(10)

In our theorem statement we have assumed that
$\\{\hat{Z}_{1}^{+}\\}_{n}\overset{p}{\to}Z_{1}^{+}$,
$\\{\hat{Z}_{1}^{-}\\}_{n}\overset{p}{\to}Z_{1}^{-}$, and
$\\{\hat{Z}_{2}\\}_{n}\overset{p}{\to}Z_{2}$, from which it now follows by
Slutsky’s Theorem that

$\displaystyle\frac{\\{\hat{Z}_{1}^{+}\\}_{n}-\\{\hat{Z}_{1}^{-}\\}_{n}}{\\{\hat{Z}_{2}\\}_{n}}\overset{p}{\to}\frac{Z^{+}_{1}-Z^{-}_{1}}{Z_{2}}=\mathbb{E}[F]$
(11)

as required. ∎

##  Appendix C Estimating Expectations in Turing

###  C.1 Standard approach

⬇

@model function model(y=2)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

end

\parnum_samples = 1000

posterior_samples = sample(model(), NUTS(0.65), num_samples)

\parf(x) = x^3

posterior_x = Array(posterior_samples[:x])

expectation_estimate = mean(map(f, posterior_x))

Full example of the estimation of an expectation with the Turing language. The
user first defines the model, then conditions it on some observed data,
computes posterior samples and then uses these samples to compute a Monte
Carlo estimate of the expectation.

###  C.2 Using ``

When we designed the API Turing largely ignored the ` `statements in the model
definition. In the meantime Turing introduced a convenience function
\eptgenerated_quantities. Given a model and $N$ samples it returns a list of
the $N$ return values generated by running the program on each sample. Note
that ` `reruns the entire \eptmodel function for each posterior sample to
compute the return value. This means that for models which have an expensive
likelihood computation the use of ` `might incur a significant overhead. It is
important to note that \eptgenerated_quantities is merely a convenience
function and does not change how Turing interprets model definitions. In fact,
the ` `function provides complimentary functionality and Turing models
generated with EPT can use this function without problems. The example from
Section~\refapd:turing_monte_carlo can be rewritten to use `:`
\beginminted[breaklines,escapeinside=||,mathescape=true,numbersep=3pt,gobble=2]eptlexer.py:EPTLexer
-x @model function model(y=2) x |$\sim$| Normal(0, 1) y |$\sim$| Normal(x, 1)
return x^3 end num_samples = 1000 posterior_samples = sample(model(),
NUTS(0.65), num_samples) expectation_estimate =
mean(generated_quantities(model(), posterior_samples))

##  Appendix D Full Example of `` Macro Transformation

The expectation

⬇

@expectation function expt_prog(y)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

return x^3

end

gets transformed into

⬇

@model function gamma1_plus(y)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

tmp = x^3

if _context isa Turing.DefaultContext

@addlogprob!(log(max(tmp, 0)))

end

return tmp

end

\par@model function gamma1_minus(y)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

tmp = x^3

if _context isa Turing.DefaultContext

@addlogprob!(log(-min(tmp, 0)))

end

return tmp

end

\par@model function gamma2(y)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

return x^3

end

\parexpt_prog = Expectation(

gamma1_plus,

gamma1_minus,

gamma2

)

The type ` `is simply used to have one common object which stores the three
different Turing models. Notice that for \eptgamma2 the function body is
identical to the original function. For ` `and \eptgamma1_minus we also have
to check in what ` `the model is executed in. Turing allows to execute the
model with different contexts which change the model behaviour. For example,
there is a \eptPriorContext which essentially ignores the tilde statements
which have observed data on the LHS. This is useful for evaluating the prior
probability of some parameters. However, by default the ` `macro ignores the
model context. As a consequence if a Turing model includes an \ept@addlogprob
macro and is executed with a ` `then it no longer calculates the log prior
probability but instead the log prior probability plus whatever value was
added with the \ept@addlogprob statement. Since we want to use the Turing
model with Annealed Importance Sampling we need to be able to extract the
prior from our model and hence we need to ensure that we do not call ` `when
executed in a \eptPriorContext. This is what the added `` clause ensures.

##  Appendix E Different Estimators for $Z_{1}^{+}$, $Z_{1}^{-}$ and $Z_{2}$

The target function $f(x)=x^{2}$ in the following expectation is always
positive:

⬇

@expectation function expt_prog(y)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

return x^2

end

Therefore, we already know that $Z_{1}^{-}=0$, so it would be wasteful to
spend computational resources on estimating $Z_{1}^{-}$. EPT allows users to
specify the marginal likelihood estimator for each of the terms in TABI
separately which means if the user knows that the target function is always
positive they can specify that 0 samples should be used to estimate
$Z_{1}^{-}$:

⬇

expct_estimate, diagnostics = estimate_expectation(

expt_prog(2), TABI(

TuringAlgorithm(AnIS(), num_samples=1000), # $Z_1^+$

TuringAlgorithm(AnIS(), num_samples=0), # $Z_1^-$

TuringAlgorithm(AnIS(), num_samples=1000) # $Z_2$

))

It is easy to see how this can be adapted to the case in which we have
$Z_{1}^{+}=0$. This interface is not just useful for avoiding unnecessary
computation, in some cases the user might also want to have different marginal
likelihood estimators for each term. This allows user to further tailor the
inference algorithm for the given target function $f(x)$.

##  Appendix F Hyperparameters for Experiments

Target-aware annealed importance sampling (TAAnIS) runs standard annealed
importance sampling twice: one time to estimate $Z^{+}_{1}$ and the other time
to estimate $Z_{2}$. For each of the problems we always use the same
hyperparameters for the annealed importance sampling algorithm both to run
AnIS and for the two estimates in TAAnIS. All of our experiments were run on a
server with 56 cores and Intel(R) Xeon(R) Gold 5120 CPUs. Running times for
the individual scripts are given in the code supplementary.

###  F.1 Posterior Predictive

For the annealed importance sampling, we use a MH transition kernel with an
isotropic Gaussian with covariance $0.5I$ as a proposal and $5$ MH steps on
each annealing distribution. We use $100$ uniformly spaced annealing
distributions. For the MCMC, we collect $5\cdot 10^{6}$ samples in total. To
parallelise sampling we run $500$ chains with $10^{4}$ samples each in
parallel, discarding the first $10^{3}$ samples as burn-in. We use a MH
transition kernel with standard Normal proposal.

###  F.2 SIR Model

For the annealed importance sampling estimators we use HMC transition kernels
with a step size of 0.05, 10 leapfrog steps and 10 MCMC steps on each
annealing distribution. We use 100 geometrically spaced annealing
distributions. For the MCMC model we collect $10^{6}$ samples in total with
Turing’s implementation of NUTS and a target acceptance rate of
65%.333https://turing.ml/dev/docs/library/#Turing.Inference.NUTS We
parallelise sampling over $10^{2}$ chains with $10^{4}$ samples and discard
the first $10^{3}$ samples as burn-in. The ground truth is computed using
importance sampling with $10^{8}$ samples and the prior as a proposal
distribution. See Equation (12) for the full SIR model including the priors.
The observed data was generated from the model described in (12) with
$\beta=0.25$, $I_{0}=100$, $N=10^{4}$ and $\phi=10$ as the overdispersion
parameter of the SIR model. We generate data for 15 time steps.

###  F.3 Radon model

We run TAAnIS and AnIS with 200 intermediate distributions and one step of the
dynamic HMC transition kernel [42, 41] on each intermediate distribution with
a step size of $0.044$. The step size was informed by running adaptive MCMC on
the target distribution.

##  Appendix G SIR Experiment

We assume we are given data in the form of observations $y_{i}$, the number of
observed newly infected people on day $i$. Fixing $\gamma=0.25$, this gives us
the statistical model

$\displaystyle\beta$
$\displaystyle\sim\text{TruncatedNormal}(2,1.5^{2},[0,\infty]),$
$\displaystyle I_{0}$
$\displaystyle\sim\text{TruncatedNormal}(100,100^{2},[0,10000]),$ (12a)
$\displaystyle S_{0}$ $\displaystyle=10000-I_{0},$ $\displaystyle R_{0}$
$\displaystyle=0,$ (12b) $\displaystyle\mathbf{x}$
$\displaystyle=\texttt{ODESolve}(\beta,\gamma,S_{0},I_{0},R_{0}),$
$\displaystyle y_{i}$
$\displaystyle\sim\text{NegativeBinomial}(\mu=x_{i},\phi=0.5).$ (12c)

Here ODESolve indicates a call to a numerical ODE solver which solves the set
of equations (3). It outputs $x_{i}$, the predicted number of newly infected
people on day $i$. We assume the observation process is noisy and model it
using a negative binomial distribution, which is parametrised by a mean $\mu$
and an overdispersion coefficient $\phi$. For an in-depth discussion about
doing Bayesian parameter inference in the SIR model we refer the reader to the
case study of Grinsztajn et al. [48]. We are further given a cost function in
terms of $R_{0}$, $\text{cost}(R_{0})=10^{12}*\text{logistic}(10R_{0}-30)$.
Intuitively, the cost initially increases exponentially with $R_{0}$. However,
the total cost also saturates for very large $R_{0}$ (as the entire population
becomes infected).

##  Appendix H Hierarchical Radon Model

The data for this problem was taken from: https://github.com/pymc-devs/pymc-
examples/blob/main/examples/data/radon.csv (the repository uses an MIT
license; the data contains no personally identifiable information). The
original data contains information about houses in 85 counties. In order to
make estimating normalization constants more tractable we reduce the number of
counties to 20. Our target function is a function of predicted radon levels
$y_{i}$ for a typical house with a basement (i.e. $x_{i}=0$) in county $i$;
$y_{i}$ is calculated using the predictive equation given in (4b). We apply
the function

$\displaystyle f(y_{i})=\frac{1}{1+\exp(5(y_{i}-4))}$

to all the predicted radon levels and then take the product of all the
$f_{i}$. Finally, to avoid floating point underflow we set a minimum value of
$1\mathrm{e}{-200}$.

##  Appendix I Multiple Expectations and Restrictions on $f(\cdot)$

The user is not restricted to defining only one expectation per model. By
specifying multiple return values the user can specify multiple expectations.
The ` `macro can recognise multiple return values and generates an expectation
for each of them. The user can then estimate each expectation independently
using \eptestimate_expectation:

⬇

@expectation function expt_prog(y)

x |$\sim$| Normal(0, 1)

y |$\sim$| Normal(x, 1)

return x, x^2, x^3

end

y_observed = 3

expt_prog1, expr_prog2, expt_prog3 = expt_prog

expct1 = expt_prog1(y_observed)

expct1_estimate, diagnostics = estimate_expectation(

expct1, method=TABI(marginal_likelihood_estimator=TuringAlgorithm(

AnIS(), num_samples=1000)))

##  Appendix J Posterior Predictive Model in EPT

The expectation from Section 5.1 can be defined in just 5 lines of code with
EPT:

⬇

@expectation function expt_prog(y)

x |$\sim$| MvNormal(zeros(length(y)), I) # $\mathbf{x}
\sim\mathcal{N}(\mathbf{x}; 0, I)$

y |$\sim$| MvNormal(x, I) # $\mathbf{y} \sim\mathcal{N}(\mathbf{y};
\mathbf{x}, I)$

return pdf(MvNormal(x, 0.5*I), -y) |$\phantom{\sim}$|# $f(\mathbf{x}) =
\mathcal{N}(-\mathbf{y}; \mathbf{x}, \frac{1}{2}I)$

end

##  Appendix K Syntax Design

Prior works have considered two families of syntax design corresponding to the
semantics required by EPF. Gordon et al. [2] define the semantics for
expectation computation via the syntax of probabilistic program’s return
expression, which is the approach we adopted in the design of EPT. Zinkov and
Shan [34] take a different route and define the expectation semantics via the
use of syntax

    expect(m, f)}
    where \mintinline

eptlexer.py:EPTLexer -xm is the program defining a measure and

    f} is the target function.
    While designing the interface of EPT we considered two different design for
    defining the target function: either letting users specify the target function
    implicitly through the return values of the function or allowing users to
    specify a target function \ept

f externally. The external function could then be passed to the ` `function
explicitly. For EPT, we decided to adopt the former of the two designs mainly
due to the simplicity of the resulting user interface and implementation. In
particular, it allows for simple to execute program transformations of the
\mintinlineeptlexer.py:EPTLexer -x@expectation macro into valid Turing
programs to represent the individual densities, and thus the ability to use
native Turing inference algorithms. Adopting the other approach would
additionally require designing and specifying the interface between the
function signature

    f(.)} and
    the values of the named random draws performed by the model
    \mintinline

eptlexer.py:EPTLexer -xm. This would result in a more complex user-facing
interface, at the slight advantage of improved compositionality of models and
functions.

##  Appendix L SIR Discussion

(a) MCMC samples.

(b) MCMC samples including burn-in samples (in black).

(c) AnIS samples.

(d) TAAnIS samples for $Z_{1}$.

Figure 6: Samples from the different algorithms for the SIR model. Note that
for Figure 5(b) some burn-in samples lie outside the boundaries of the plot
but we adjusted the axis limits so that they are the same for all plots to
allow for easier comparison.

In the SIR experiment AnIS achieved a significantly lower RSE than MCMC even
though both are non-target-aware. Figure 6 shows samples from the different
algorithms. The TAAnIS samples for $Z_{1}$ visualise well in which regions of
parameter space both the posterior and the target function have sufficient
mass ($\beta\in[0.5,2.0]$). The samples from AnIS and MCMC suggest that most
of the posterior mass is located in the interval $\beta\in[0.3,0.7]$. However,
AnIS also generates a significant amount of samples in the parameter region
$\beta\in[1.0,1.5]$. The samples in this second “mode” are directly in the
region of the target-aware samples. Further, the plots suggest that AnIS
generates more samples in this regions than MCMC which is what allows AnIS to
achieve a lower RSE. However, it seems that the AnIS represents the second
“mode” disproportionally. Specifically looking at the burn-in samples from
MCMC in Figure 5(b) shows that MCMC will converge to the parameter space in
$\beta\in[0.3,0.7]$ even if the initial parameter samples are around
$\beta\in[1.0,1.5]$. This indicates that this is not a failure of MCMC to
detect another mode but rather that there is negligible posterior mass in that
parameter region. Therefore the better performance of AnIS compared to MCMC
seems to occur mostly because AnIS got lucky by accidentally generating
samples in the right parameter region.

###  L.1 A Note on MCMC ESS

The SIR experiment provides a good example of how the MCMC ESS [44] is
unreliable for our use case. As detailed in Section F.2 for MCMC we run $100$
chains with $10,000$ samples each. This is replicated $5$ times to get
estimates on the variability in behaviour. After discarding the burn-in
samples for each chain the $5$ replications give us the following final ESS
estimates: $[631,360;\,805,868;\,873,269;\,665,683;\,5,114]$. We observe that
all but one replication give disproportionally high ESS estimates. We found
that the replication which gives a more conservative ESS estimate of $5,114$
is the replication which generated samples in the parameter region
$\beta\in[1.0,1.5]$ (see Figure 5(a)). More importantly, the MCMC ESS
estimates do not seem to show any correlation with the RSE values (see Figure
4) which is the more important metric because it directly measures the error
in our estimate. Therefore, we decided against using the MCMC ESS in our
evaluation because it can give the impression that MCMC is performing well
when it is actually failing dramatically (in terms of RSE).

##  Appendix M Effective Sample Size

(a) Gaussian Posterior Predictive.

(b) SIR.

(c) Radon.

Figure 7: Individual ESS values as defined in Section 5 for the three
different experiments. Instead of taking
$\min(\text{ESS}_{Z_{1}},\text{ESS}_{Z_{2}})$ for TAAnIS and
$\min(\text{ESS}_{\text{AnIS retargeted}},\text{ESS}_{\text{AnIS}})$ for AnIS
we plot each value individually. Plotting each ESS value separately shows that
the performance of AnIS is severely limited by its ability to generate samples
in regions in which the target function $f(x)$ is large. This is indicated by
the low values for $\text{ESS}_{\text{AnIS retargeted}}$.

##  Appendix N Positive and Negative Target Functions

To demonstrate that EPT is also beneficial for target functions which are
positive and negative we provide a brief description of a synthetic
experiment. We assume the following model which gives us a banana shaped
density (see Figure 8):

⬇

@expectation function banana()

x1 |$\sim$| Normal(0, 4)

x2 |$\sim$| Normal(0, 4)

@addlogprob!(banana_density(x1, x2))

return banana_f(x1, x2)

end

\parbanana_density(x1, x2) = -0.5*(0.03*x1^2+(x2/2+0.03*(x1^2-100))^2)

Note that there is no observed data in this experiment which is why we chose
to express the banana distribution as an unnormalized density (i.e. use the `
`primitive). Our target function is given by \vspace-8pt

⬇

function banana_f(x1, x2)

cond = 1 / (1 + exp(50 * (x2 + 5)))

return cond * (x1 - 2)^3

end

Note that the target function can be positive and negative. Figure 8 shows the
RSE for TAAnIS and AnIS. We used an MH transition kernel and 200 intermediate
potentials for the Annealed Importance Sampling estimators. The RSE of AnIS
does not improve because it fails to generate samples in the regions in which
the target $f(x)$ is large. [15] provide a comparison to MCMC on a similar
problem so we omit it here.

Figure 8: Banana experiment. [Left] Heatmap of the density of the model.
[Right] Relative Squared Error for TAAnIS and AnIS.