# Functional probabilistic programming for scalable Bayesian modelling

Jonathan Law  School of Mathematics, Statistics & Physics
Newcastle University, U.K.  National Innovation Centre for Data  Darren J.
Wilkinson  School of Mathematics, Statistics & Physics
Newcastle University, U.K.  The Alan Turing Institute

###### Abstract

Bayesian inference involves the specification of a statistical model by a
statistician or practitioner, with careful thought about what each parameter
represents. This results in particularly interpretable models which can be
used to explain relationships present in the observed data. Bayesian models
are useful when an experiment has only a small number of observations and in
applications where transparency of data driven decisions is important.
Traditionally, parameter inference in Bayesian statistics has involved
constructing bespoke MCMC (Markov chain Monte Carlo) schemes for each newly
proposed statistical model. This results in plausible models not being
considered since efficient inference schemes are challenging to develop or
implement. Probabilistic programming aims to reduce the barrier to performing
Bayesian inference by developing a domain specific language (DSL) for model
specification which is decoupled from the parameter inference algorithms. This
paper introduces functional programming principles which can be used to
develop an embedded probabilistic programming language. Model inference can be
carried out using any generic inference algorithm. In this paper Hamiltonian
Monte Carlo (HMC) is used, an efficient MCMC method requiring the gradient of
the un-normalised log-posterior, calculated using automatic differentiation.
The concepts are illustrated using the Scala programming language.

A git repository containing runnable example code for the probabilistic models
and an illustrative implementation of forward and reverse mode automatic
differentiation is available at https://git.io/probprog.

##  1 Introduction

Bayesian inference using probabilistic programming is becoming of interest to
large industry-backed artificial intelligence (AI) labs such as Google, Uber
and Stripe (TensorFlow, 2018; Uber, 2018; Bryant, 2018), despite deep learning
continuing to solve practical problems with large amounts of labelled data.
Deep neural networks are typically used as a “black box”, and state of the art
networks can involve many parameters (Huang et al., 2017). Some dense
convolutional neural networks can have millions of parameters. The values of
these parameters are impossible to reason about independently and require a
large amount of data to learn.

In a Bayesian model, the likelihood and prior distributions for each parameter
must is specified using the judgement of the statistician and/or subject
matter experts. This results in an interpretable model which can be used to
explain relationships present in the data and perform accurate predictions
with appropriately quantified uncertainty. Bayesian methods are especially
useful when an experiment has only a small number of observations, for
instance in an A/B test using a small number of highly valuable customers or
in a clinical trial. Bayesian methods provide a natural encoding for
uncertainty; if a small number of observations are available to inform the
likelihood then the posterior distribution will reflect the prior. Specifying
uncertainty honestly along with predictions can allow practitioners to make
informed decisions. Well specified Bayesian models ensure industries such as
insurance and finance, which price products and shares based on mathematical
and statistical models, can be transparent in the decisions that they make.

Traditionally, parameter inference in Bayesian statistics has involved
constructing bespoke MCMC schemes for each newly proposed statistical model.
This results in new models being explored more slowly than practitioners would
like, sometimes leading to plausible models not being considered since
efficient inference schemes are challenging to develop or implement.
Probabilistic programming aims to reduce the barrier to performing Bayesian
inference by developing a domain specific language for model specification
which does not rely on a detailed understanding of the parameter inference
algorithms. This means that new models can be explored easily, without
significant programming effort.

There exist many domain specific languages (DSLs) for performing Bayesian
inference, the most well known within the statistical community are BUGS (Lunn
et al., 2000) and JAGS (Plummer et al., 2003) in which a simple syntax can be
used to specify complex hierarchical models. BUGS and JAGS both use variants
of Gibbs sampling to perform inference. Stan is another popular DSL which uses
similar syntax to BUGS and JAGS. Stan uses Hamiltonian Monte Carlo (HMC) with
automatic differentiation for efficient sampling (Carpenter et al., 2016) but
is restricted to models with continuous parameters. Stan also implements
approximate variational inference schemes.

Many existing frameworks for probabilistic programming, including BUGS, JAGS
and Stan use DSLs which are in turn compiled to lower-level languages. For
instance in a Stan script, arbitrary C++ can be written to define a model or
inference scheme, but in a C++ program the Stan DSL for defining a model
cannot be easily used. This is in contrast to an embedded DSL, where the DSL
is written directly in the target language and can be easily deployed inside a
larger program.

TensorFlow is a program created at Google for deep learning and neural
networks (Abadi et al., 2015). Google has developed TensorFlow Probability
(formally known as Edward), a Python library for Bayesian inference utilising
the TensorFlow architecture (TensorFlow, 2018). TensorFlow probability uses
automatic differentiation and GPU acceleration from TensorFlow for variational
inference and HMC. Ride-sharing company Uber have developed a Python library
called Pyro (Uber, 2018) for probabilistic programming marrying deep learning
and Bayesian inference using PyTorch (Paszke et al., 2017) for automatic
differentiation and GPU acceleration. Pyro is implemented in Python, a popular
general purpose programming language. These probabilistic programming
languages can be easily deployed in existing Python code bases, however the
syntax for implementing a model is more verbose than that of BUGS, JAGS and
Stan. Hakaru (Narayanan et al., 2016) is an example of a probabilistic
programming language with both an embedded DSL written in Haskell and a more
accessible DSL written in the style of Python. This design choice was to
minimise the learning curve of those coming from languages similar to Python.

Within the theoretical computer science community, an effort has been made to
formalise the semantics of Bayesian inference in category theory by defining a
probability monad grounded in measure theory (Giry, 1982; Lawvere, 1962). A
practical implementation of Bayesian updating using the sampling based
probability monad, built upon the State monad with applications to robot
localisation and mapping is presented in Park et al. (2005). A stochastic
lambda calculus is introduced in Ramsey and Pfeffer (2002) along with a review
of the probability monad as presented in earlier literature, using a Haskell-
like notation. More recently, the use of the Giry monad as a foundation for
probabilistic programming has been questioned. The issue is that the Giry
monad is defined on the category of measurable spaces, but this category is
not Cartesian closed. Since closure is desirable, alternative probability
monads over Cartesian closed categories are being explored, such as quasi-
Borel spaces (Heunen et al., 2017) and the Kantorivich monad (Fritz and
Perrone, 2017). The development of a composable probabilistic programming
language using modern functional programming approaches is described in
Ścibior et al. (2018a).

Rainier (Bryant, 2018) is an embedded probabilistic programming language
written in Scala (Odersky et al., 2004) developed by the online payment
processing company, Stripe. Rainier provides a compositional, monadic syntax
for the specification of models. Each model has an associated compute graph
which can be used to calculate the un-normalised log-posterior and its
gradient using automatic differentiation. The log-posterior and its gradient
is then used for HMC (Duane et al., 1987). The monadic interface is familiar
to functional programmers and means that well researched and developed
features of existing functional programming languages can be used when
developing probabilistic programs. Crucially, models are are regular values
within the host programming language, which can be manipulated and composed
using regular language features. This is demonstrated in an example presented
in this paper, a random effects model (see section 8.3). Each class of the
random effects model is modelled using an independent linear model, and the
linear model code from example 8.1 is reused. This code reuse is
straightforward when using a well developed programming language such as Scala
and this is one of the main advantages of embedded probabilistic programming
languages.

Scala is a programming language which runs on the Java Virtual Machine (JVM)
and provides support for functional programming. Scala makes heavy use of
function composition, higher-order functions and type classes to build large
applications from simple functions. These features allow programmers to build
complex programs whilst maximising code-reuse and minimising lines-of-code.
This style of programming is a departure from object-oriented or imperative
programming and as such can have a steep learning curve for those coming from
languages which emphasise those approaches. This paper introduces the concepts
required to understand monads in functional programming and how they can be
used to build an embedded DSL for Bayesian statistical modelling.

Embedded DSLs are straightforward to write in functional languages using, for
example, free monads and generalised algebraic data types (Swierstra, 2008).
Ścibior et al. (2015) implements an algebra in Haskell for operating on
distributions using a Generalised Algebraic Datatype and the Free Monad
(Ścibior et al., 2015). This allows the implementation details of complex
programs (in this case statistical inference algorithms) to be decoupled from
the model specification. Since the DSL is implemented in the same programming
language, it is straightforward to use within existing programs. The
importance of building DSLs for scientific computing is expanded upon in an
article advocating the use of Lisp in bioinformatics (Khomtchouk et al.,
2016). Lisps are a collection of dialects of dynamically typed functional
programming languages, which share some similarities with functional
programming in strongly typed functional languages such as Scala and Haskell.

##  2 Category Theory

This section introduces the aspects of category theory (Awodey, 2010; Barr and
Wells, 1990) required to build a probabilistic programming language. The
concepts and their uses are reified in the Scala language in section 3.
Category theory is concerned with describing abstract structures in a very
general way. For an introduction to category theory aimed at programmers, see
Milewski (2018).

######  Definition.

A _category_ $\mathcal{C}$ is a collection of objects and morphisms which go
between them such that for the objects $X,Y,Z\in\mathcal{C}$ and morphisms
$g:X\rightarrow Y$ and $h:Y\rightarrow Z$,

  1. 1.

There must exist a morphism $f:X\rightarrow Z$ which is the composition of $g$
and $h$, $f=h\circ g$

  2. 2.

Each object in a category must have an identity morphism written as
$\textrm{id}_{X}:X\rightarrow X$

  3. 3.

Composition of morphisms must be associative. For all $f:X\rightarrow Y$,
$g:Y\rightarrow Z$ and $h:Z\rightarrow W$ then $h\circ(f\circ g)=(h\circ
f)\circ g=h\circ f\circ g$

  4. 4.

For every morphism $g:X\rightarrow Y$ then $\textrm{id}_{Y}\circ
g=g=g\circ\textrm{id}_{X}$.

The collection of objects in a category $\mathcal{C}$ is often denoted
$\textrm{obj}(\mathcal{C})$ and morphisms $\textrm{hom}(\mathcal{C})$. The set
of morphisms from $X$ to $Y$ is written as $\textrm{hom}_{\mathcal{C}}(X,Y)$.

An example of a category is Set, the category of sets: each set is an object
and the morphisms are functions between the sets, the identity morphism is the
identity function. In functional programming, the category most often under
consideration is similar to the category of sets: each object corresponds to a
specific datatype and the morphisms represent functions between the types.

###  2.1 Functors

A functor is a structure preserving mapping from one category to another.

######  Definition.

If $\mathcal{C}$ and $\mathcal{D}$ are categories, then a _functor_
$F:\mathcal{C}\rightarrow\mathcal{D}$ takes each object $X$ in $\mathcal{C}$
to an object $F(X)$ in the category $\mathcal{D}$ and each morphism,
$f:X\rightarrow Y$ in $\mathcal{C}$ to a morphism in $\mathcal{D}$,
$F(f):F(X)\rightarrow F(Y)$ such that

  1. 1.

$F(\textrm{id}_{X})=\textrm{id}_{F(X)}$ for each $X$ in $\mathcal{C}$

  2. 2.

$F(g\circ f)=F(g)\circ F(f)$ for all morphisms, $f$ and $g$ in $\mathcal{C}$.

Since functional programming is mainly concerned with the category Set,
functors are typically endofunctors, $F:\textbf{Set}\rightarrow\textbf{Set}$.
An endofunctor is a functor from a category to itself,
$F:\mathcal{C}\rightarrow\mathcal{C}$. Note that every category $\mathcal{C}$
has an associated identity endofunctor,
$\textrm{Id}_{\mathcal{C}}:\mathcal{C}\rightarrow\mathcal{C}$ mapping all
objects and morphisms to themselves.

###  2.2 Natural Transformation

A functor is a morphism between two categories (which itself has morphisms
between its objects). A _natural transformation_ is a morphism between
functors which preserves the structure of the categories.

######  Definition.

If $F:\mathcal{C}\rightarrow\mathcal{D}$ and
$G:\mathcal{C}\rightarrow\mathcal{D}$ are both functors between the categories
$\mathcal{C}$ and $\mathcal{D}$ then a _natural transformation_
$\alpha:F\Rightarrow G$ is a family of morphisms such that

  1. 1.

For all $X\in\mathcal{C}$ then $\alpha_{X}:F(X)\rightarrow G(X)$ is a morphism
in $\mathcal{D}$ called the component of $\alpha$ at $X$

  2. 2.

For each morphism $f\in\mathcal{C}$, $f:X\rightarrow Y$ then $\alpha_{Y}\circ
F(f)=G(f)\circ\alpha_{X}$.

Figure 1 shows a commutative diagram of the second natural transformation law.
Following the arrows of the commutative diagram for any path from any given
start to end point must be equivalent.

Figure 1: Commutative diagram of the natural transformation

###  2.3 Monads

######  Definition.

A monad on the category $\mathcal{C}$ is an endofunctor,
$T:\mathcal{C}\rightarrow\mathcal{C}$ with two natural transformations,
$\eta:\textrm{Id}_{c}\rightarrow T$ called unit and $\mu:T\circ T\rightarrow
T$ called multiplication (or counit) such that

  1. 1.

$\mu\circ T\mu=\mu\circ\mu T$

  2. 2.

$\mu\circ T\eta=\mu\circ\eta T=\textrm{Id}_{T}$

$T\mu:T(T\circ T)\rightarrow T$ is a natural transformation defined by taking
the natural transformation $\mu$ and the endofunctor,
$T:\mathcal{C}\rightarrow\mathcal{C}$ then $(T\mu)_{C}=T\mu_{C}$. The natural
transformation $\mu T$ is a similar transformation with $T$ composed from the
right, $\mu T:(T\circ T)T\rightarrow T$. Figure 2 shows commutative diagrams
of the monad laws.

(a)  (b)

Figure 2: (a) Commutative diagram expressing the first condition of the monad
laws (b) Commutative diagram expressing the second condition of the monad laws

##  3 Principles of Functional Programming

In functional programming, a program is comprised of many small functions each
of which are easy to test and reason about in isolation. These small functions
are composed together to form a larger program, which in practice is more
difficult to reason about and test. The programmer gains confidence in the
correctness of the program by ensuring each individual function in the
composition is correct. This section provides an introduction to the
functional programming principles required to build a probabilistic
programming language guided by category theory, with examples in Scala. For a
more comprehensive introduction to functional programming in Scala, see
Chiusano and Bjarnason (2014).

###  3.1 Referential Transparency and Pure Functions

_Referential transparency_ is a fundamental principle of functional
programming. A function is referentially transparent if it produces the same
output value for the same function parameters upon repeated evaluation, and
produces no side effects. This enables the programmer to directly replace any
referentially transparent function with the value of that function without
changing the output of the program. This is called reasoning by substitution.
A function which is referentially transparent is often called deterministic.

A _pure function_ is a function which is referentially transparent, inculpable
and total. A function is total if it returns a value (not an exception) for
every possible input. Pure functions are more straightforward to test, reason
about and reuse by composing with other pure functions. A function is called
_inculpable_ if it is free from side effects.

Functions which have _side effects_ are not referentially transparent. An
example of a side effect is printing to the screen, reading from a file or
database, accepting input from the user or mutating a non-local variable. Some
or all of these actions are needed for a program to be useful. Functional
programming makes side effects explicit and good functional design ensures
that side effects are only performed when needed, typically at the edges of a
program.

An example of a common side effect in statistical programming is generating a
random number. In Scala (which is not a _pure_ functional language) a random
double between 0 and 1 can be generated using scala.util.Random.nextDouble().
When calling this function it is impossible to know which random number will
appear next without knowing the current state of the random number generator.
Additionally, when this function is called multiple times in a program it will
return different values. This is due to the fact that the function has the
side effect of mutating the global random number state. The state monad can be
used to implement a purely functional random number generator and this is
expanded upon in section 3.10.

###  3.2 Static Types

Functions can accept a range of parameters and each of these parameters has a
type. For instance the function $f(x)=x+1$ is a mathematical function which
adds one to its parameter, $x$. In this case, the parameter type is implicitly
a numeric type as the method $+$ is associated with numbers. The return type
of the function $f$ is also a numeric type. Types can be either inferred (as
in the function $f$) or they can be made explicit. The example Scala code
below shows a the function $f$ written with an explicit type declaration of
Int => Int.

⬇

val f: Int => Int = x => x + 1

If this function was applied to a String data-type, the program would not
compile. In a language with dynamic types, the program would run and possibly
try to convert the String into an Int this can lead to runtime errors which
are typically difficult to debug as they can potentially go unnoticed.

Static types are not exclusive to functional programming, nor are they
required. However, programs using static types are typically safer and easier
to reason about; Scala and Haskell are examples of functional programming
languages with static types. However, both Scala and Haskell have _type
inference_ , which means that explicit type annotations are not required in
situations where the compiler can unambiguously resolve all types; this
provides most of the advantages of static typing without the verbosity
associated with statically typed languages lacking type inference.

###  3.3 Immutable Data

_Immutable data_ is an important aspect of functional programming. An
immutable value or data structure cannot be modified in place. This is a
strict requirement which ultimately makes code easier to reason about. However
this means that imperative programming constructs, such as for-loops and
while-loops can not be used since, for instance, the index of the loop can not
be advanced (by mutating the variable representing the index) as would be
required in a for-loop.

###  3.4 Higher Kinded Types

Collection types such as a singly linked lists, represented by List, are
sometimes referred to as a type constructors (or data constructors) since they
are used to construct concrete types such as List[Int] or List[Double]. List
is similar to other collection types such as Vector or Option. Vectors and
lists are collections of zero or more elements. Whereas an Option is a
container for turning a partial function into a total function by explicitly
representing failure cases using None. Option can be thought of as a
collection which can be either empty (None) or contain exactly one element.

The implementations of List and Vector are different and both have distinct
performance characteristics. A singly linked list provides $O(1)$ access to
the head of the list but $O(n)$ access time to a specific element. Whereas
Vector in Scala is implemented as a trie (Fredkin, 1960) with a branching
factor of 32. This means that random access to elements in a Vector is faster
than in a List. Consequently, concrete collections should be appropriately
chosen for a given purpose.

These collections (and other type constructors) share common abstractions
which can be used to unify them and then write code which is polymorphic not
just in the type, but in the type constructor. This allows for greater code
re-use. Higher kinded types provide powerful abstractions leading to cleaner,
more elegant code. Very few mainstream programming languages include the
necessary language support for using higher kinded types to safely abstract
over generic type constructors, with Scala and Haskell arguably being the best
known examples of languages which do.

###  3.5 Higher Order Functions

In programming languages which have functions as first class citizens,
functions can be passed as an argument to other functions. Functions which
accept other functions as arguments or return functions are called _higher
order functions_.

A common operation performed on collections using a for-loop is to apply a
function, f, to each element in the collection. The map function can be used
to apply a function to a value in a context, and a context, as in this case,
can be a collection. The signature of the map function is

⬇

def map[F[_], A, B](fa: F[A])(f: A => B): F[B]

map is a higher order function which accepts a function, f: A => B and applies
it to each element of F. The term F is used to stand for any construct which
can be mapped over using a function, this construct is known as a Functor (see
section 3.6). The relationship with the category theoretic functor is made
clearer if the map is restated such that it transforms the function f: A => B
into g: F[A] => F[B].

Another common operation for collections is a reduction, and this can be
achieved by applying a binary function to each pair of elements of the
collection in order from either the left or right. A higher order function
which reduces values of the list from the left starting with an initial seed
value is foldLeft, the signature of foldLeft is

⬇

def foldLeft[F[_], A, B](fa: F[A])(b: B)(f: (B, A) => B): B

map and foldLeft are implemented using recursion, but in order to avoid stack
overflows they are implemented using tail-calls which are straightforward for
the compiler to optimise. In contrast to imperative for-loops these higher
order functions require no auxiliary variables such as counters. They are
typically shorter and express the intent of the programmer more clearly.

###  3.6 Functors

A functor represents a type constructor which can be mapped over. A functor
can be used to lift a function between two types into a new function between
two types in a context, where the context refers to the type constructor. The
signature of the map function for each collection List, Option or Vector can
be written as follows:

⬇

def map[A, B](fa: Option[A])(f: A => B): Option[B]

def map[A, B](fa: Vector[A])(f: A => B): Vector[B]

The function f: A => B is lifted into the context of List, Option and Vector
respectively. The higher order function map allows f to operate on an element
in a context. The map functions defined above are identical apart from the
type constructor, the type constructor can thus be abstracted over into the
functor type class:

⬇

def map[A, B](fa: F[A])(f: A => B): F[B]

}

The type constructor, denoted as F[_], has been abstracted over using a
_higher kinded type_ , here a functor. Now functions can be defined with a
functor context bound on the type constructor:

⬇

map(fa)(a => a + 1)

}

This function can be applied to any data constructor which has a functor
instance and which contains Numeric values. In order to use the function,
addOneGen with a concrete collection type, a functor instance must be defined
by defining the map function. The Scala library
Cats111https://github.com/typelevel/cats provides type class instances for
many of the built in Scala types constructors.

###  3.7 Applicatives

The data constructors considered so far are also applicatives (McBride and
Paterson, 2008). An applicative is a special functor which has a pure method
which lifts a value into the context of the Applicative and a ap which can be
used to apply a function to a value in the Applicative context (the F[_])

⬇

def pure[A](a: A): F[A]

def ap[A, B](fa: F[A])(f: F[A => B]): F[B]

}

To define pure, use the type constructors for the collections:

⬇

def pure[A](x: A): Option[A] = Some(x)

def pure[A](x: A): Vector[A] = Vector(x)

Where :: is used as an infix operator which can be used to construct a singly-
linked list by prepending a single element on the left to a list on the right.
In Cartesian closed categories (CCCs), like Set, ap can be defined in terms of
zip and map. In category theory an applicative defined in terms of ap is a
_lax closed functor_ and an applicative defined in terms of zip is a _lax
monoidal functor_. In CCCs, these are equivalent. As such, the applicative
typeclass can be defined in terms of the functions zip, map and pure. The type
signature of zip is given by

⬇

def zip[A, B](fa: F[A], fb: F[B]): F[(A, B)]

Then ap can be derived

⬇

map(zip(f, fa)){ case (fi, a) => fi(a) }

The definition of zip for the List collection can be defined recursively as
follows,

⬇

case (Nil, lb) => Nil

case (la, Nil) => Nil

case (a :: ta, b :: tb) => (a, b) :: zip(ta, tb)

}

If the lists are different lengths, zip will return a list of tuples the same
length as the shortest list.

###  3.8 Monads

The collections considered so far are actually monads. A monad represents a
value in a context which can be manipulated in a consistent way (Wadler,
1995). Monads are integral to safe purely functional programming and are used
to encapsulate unsafe program behaviour such as IO (input output), async
programs, random number generators and partial functions. Suitable monads used
in Scala for these unsafe behaviours include the IO monad (implemented in Cats
Effect222https://typelevel.org/cats-effect/), Future, Rand and Option
respectively. These monads have unsafe methods which can be used to access the
contents of the monad and hence perform the side effect. The unsafe access
method is only performed at the end of the program when the side effect is
desired, until then combinations of map and flatMap are used to manipulate the
values inside of the monadic context.

To define a monad in Scala, simply define the two natural transformations
required. The natural transformation corresponding to
$\eta:\textrm{Id}_{C}\rightarrow T$ is pure, inherited from the Applicative
type class. The second natural transformation, $\mu:T\circ T\rightarrow T$, is
referred to as join

⬇

def join[A](fa: F[F[A]]): F[A]

}

Typically, in functional programming monads are defined in terms flatMap

⬇

join(map(fa)(f))

which is a map followed by a join or flatten. The monad can equivalently be
defined in terms of flatMap and pure

⬇

def flatMap[A, B](fa: F[A])(f: A => F[B]): F[B]

def pure[A](x: A): F[A]

}

Note that map can be written in terms of flatMap and pure and ap can be
written in terms of flatMap and map. This means when a monad instance is
defined for a given type, only flatMap and pure must be defined and the
methods from functor and applicative are available. Hence all monads in CCCs
are applicative functors, though the converse is not true. The ubiquity of
monads in functional programming derives from the power of flatMap, which
enables composition of functions returning values in a monadic context.
Explicitly, functions f: X => F[Y] and g: Y => F[Z] don’t compose directly,
since the types don’t align, but for a monad F[_], they can be composed using
flatMap to give a function X => F[Z]. The ability to chain together such
monadic functions turns out to be exceptionally useful.

Option is a monad which can be used to turn a partial function into a total-
function by capturing error states in the type information. Consider several
functions which could fail:

⬇

if (x < 0) None else Some(math.sqrt(x))

def log(x: Double): Option[Double] =

if (x <= 0) None else Some(math.log(x))

In order to compose these two functions, a flatMap can be used:

⬇

sqrt(x).flatMap(y => log(y))

The same principal can be extended to more than two applications:

⬇

sqrt(x).flatMap(y => log(y)).flatMap(z => log(z))

This can be visually improved using a for-comprehension which is syntactic-
sugar for applications of flatMap and map:

⬇

y <- sqrt(x)

z <- log(y)

res <- log(z)

} yield res

The functions returning an Option (or any monad) are placed on the right side
of the arrow <- extracting the value from the monadic context. The values on
the left represent the value inside of the monadic context and can be
referenced in the same for-comprehension. The function composeFor is
equivalent to composeMore. The for-comprehension is “desugared” by the
compiler into a chain of flatMaps and maps in a very early compiler pass. The
equivalent to a for-comprehension in Haskell is known as `do`-notation.

In general, _different_ monads do not compose (Jones and Duponcheel, 1993),
hence functions with different monadic contexts can not be mixed inside of a
for-comprehension. However, different contexts are often required, and for
certain combinations, monad transformers can be used which combine monadic
contexts (Liang et al., 1995). Alternatively, extensible effects and the Freer
monad can be used for combining effects (Kiselyov et al., 2013).

###  3.9 Monadic Collections

Modern hardware is parallel, with most laptops containing CPUs with multiple
cores. The sequential collection types considered so far, List and Vector do
not take advantage of multithreading. In Scala there are parallel collections
which share the same monadic interface, and this enables functions to be
written in a generic way and easily parallelised to take advantage of
environments with multiple CPU cores.

Any collection can be transformed into a parallel collection, by calling .par.
Higher order functions defined on the parallel collections will run in
parallel. map and filter are naturally parallel and are evaluated by splitting
the collection into roughly equal sized chunks and evaluating a sequential map
or filter on each partition using different threads. Reduction operations need
more careful consideration. foldLeft and foldRight are inherently sequential,
however reduce can be implemented using an associative binary reduction
function and then the reduction can be performed in parallel using a parallel
tree reduction.

Consider evaluating the likelihood of a linear model, the specification of the
model is presented in example (8.1). The likelihood can be written as:

$p(y|\psi)=(2\pi\sigma^{2})^{-\frac{N}{2}}\exp\left\\{-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-(\beta_{0}+\beta_{1}x_{i}))^{2}\right\\}$

where $\psi=\\{\beta_{0},\beta_{1},\sigma\\}$. The log-likelihood can then be
written as:

$\log
p(y|\psi)=-\frac{N}{2}\log(2\pi\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-(\beta_{0}+\beta_{1}x_{i}))^{2}$

This log-likelihood can be evaluated in Scala using a map followed by a reduce

⬇

map { case (x, y) => Gaussian(beta0 + beta1 * x, sigma).logPdf(y) }.

reduce(_ + _)

Where ys: F[(Double, Double)] is a collection of tuples containing the
dependent and independent variables. This code is completely agnostic to
whether the collection being operated on is serial or parallel. If ys is a
serial collection, then map and reduce will execute serially on a single core.
However, if ys is a parallel collection, then map and reduce will execute on
the collection in parallel on all available cores. In practice, the parallel
reduction will only speed up processing significantly if the collection of
data is large, or the individual likelihood functions are expensive to
evaluate.

If the data contained in the collection is too large to fit into the memory of
a single computer, then higher order functions such as map and reduce can be
used for distributed computing (Dean and Ghemawat, 2008). This pattern was
named MapReduce after the familiar higher order functions. MapReduce is
commonly used when a large amount of input data which does not fit in the
memory of a single node requires processing. The large dataset is partitioned
onto several worker nodes, orchestrated by a single master node. During the
map stage, each worker node then performs a function independently on each
element of the data. In the reduce stage a binary operator is used to combine
the elements. The map and reduce steps can be performed by multiple nodes in
parallel, meaning computationally intensive tasks can be sped up by increasing
the number of worker nodes. In practice the speed of computation does not
increase linearly, because of the additional cost of messaging between nodes.

MapReduce is a general pattern which can be used to implement many algorithms,
including the page rank algorithm used by Google to rank search results (Brin
and Page, 1998). MapReduce aims to make distributed programming easier by
removing the possibility of race-conditions. Listing LABEL:lst:map-reduce
shows a word count implementation written in the MapReduce style. The map
stage assigns the number one to each word in the corpus of documents. groupBy
is used to group the words into key-value pairs where the key is a single word
and the value is a collection of tuples containing the same word and the
number one. The reduce stage is implemented using map and reduce, the map
takes each tuple containing each word and collection and performs the binary
reduction function + to add the numbers associated with each word. The reduce
function is a fold without an initial value, as such it can only be performed
on a non-empty collection.

⬇

map(w => (w -> 1)).

toList.

groupBy(_._1).

map { case (w, ws) => (w, ws.map(_._2).reduce(_ + _)) }

Listing 1: Illustrative MapReduce program written using Scala

Apache Spark (Apache, 2018b) is a sophisticated platform for distributed
computing for big data. Spark utilises a high level functional interface
written in Scala allowing parallel processing for large datasets and streaming
data. Spark features APIs (application programming interfaces) for Scala,
Python and Java, and each API utilises higher order functions which can be
used to manipulate data across a cluster of compute nodes. Spark primarily
stores data in main memory of compute nodes, whereas Hadoop (Apache, 2018a) (a
map reduce implementation) stores the data on the hard disk drive (HDD) or
solid-state drive (SSD) of the worker nodes. This means that Spark is
typically orders of magnitude faster for the same tasks. Apache Spark has
additional capabilities, including batch processing of streaming data and a
library for machine learning called MLlib. The main collection type in Apache
spark is the resilient distributed dataset (RDD) which is an immutable
distributed monadic collection partitioned onto different nodes in the Spark
cluster. The RDD can be operated on using the higher order functions familiar
to users of the built in Scala collections.

###  3.10 Purely Functional Pseudo Random Number Generator

The State monad can be used to keep track of internal state which is
cumbersome to pass around as an argument. The state monad is a constructor
with a function from the current state to a tuple containing a value and an
updated state:

⬇

case class State[S, A](run: S => (S, A))

Since it is a monad, there are implementations of flatMap and pure:

⬇

val (st, a) = run(x)

f(a).run(st)

}

def pure[B](a: B) = State(st => (st, a))

The State monad can be used to implement a purely-functional pseudo-random
number generator (PRNG) by threading an PRNG state through multiple function
calls.

Sequences of apparently random numbers are generated by mutating an internal
PRNG state each time a function returning a random number is evaluated. Each
time a pseudo-random number is asked for a different value is returned and the
internal PRNG state is mutated

⬇

randNaive

// 0.26167644649515154

randNaive

// 0.4732297985296451

This function violates referential transparency (note that a def is needed
here rather than a val since val is eagerly evaluated) since each time
randNaive returns different values at each invocation. In order to develop a
purely functional random number generator the internal state can be passed
around explicitly. One of the simplest implementations of a pseudo random
number generator is a linear congruential generator. This generator is defined
by a recursive relation

$X_{n+1}=(aX_{n}+c)\operatorname{mod}m$ (1)

$m>0$ is the modulus, $0\leq c<m$ is the increment and $0<a<m$ is the
multiplier and the initial seed is $0\leq X_{0}<m$. It is challenging to
choose the parameters $a,c$ and $m$ to produce a long apparently random chain,
a problem which is not considered here. The generator can be implemented in
Scala as follows:

⬇

a * state + c

def toDouble(x: Long) =

(x >>> 11) * math.pow(2, -53)

def myGenerator(state: Long): (Long, Double) = {

val newState = linearGenerator(6364136223846793005L,

1442695040888963407L, state)

(newState, toDouble(newState))

}

The linear generator is initialised without a modulus, which is equivalent to
the modulus being $2^{64}$ (since a Long is 64 bits in Scala). The function
toDouble performs a logical right bit-shift by 11 places, then divides by
$2^{53}$ (since $53=64-11$). The right shift of the Long value by 11 bits and
the subsequent division ensures there is no numeric overflow, and avoids sign
issues that would arise if the first bit was set. This has the effect of
transforming the seed to be a Double between zero and one. In order to
generate a new number myGenerator must be given an initial value, which is the
initial state of the of the random number generator:

⬇

// (3730281199248434349,0.2022189490103492)

This function returns the next state on the left of the tuple and a “random”
number between zero and one. The function can be initialised from a “random”
or chosen seed. In practice, pseudo-random number generators are often more
complex than this; a well used PRNG with a long non-repeating chain is the
Mersenne Twister algorithm (Matsumoto and Nishimura, 1998).

A function to draw a pseudo-random number can be now be written as:

⬇

def randDouble: Rand[Double] =

State(myGenerator)

Then in order to sample several numbers, flatMap can be used (hidden in a for-
comprehension):

⬇

firstNumber <- randDouble

nextNumber <- randDouble

} yield (firstNumber, nextNumber)

println(res.run(13515670))

// (2060103227662993080,(0.2022189490103492,0.11167841974883075))

The State monad passes along the state to the next invocation of randDouble.
In order to extract the pseudo-random numbers from the state monad, the
function run is called with a starting seed. The result is on the right of the
tuple with the next state on the left. This construct forms a referentially
transparent PRNG. Note also, that the same initial seed is used as in the
previous example and hence the same pseudo-random number is produced. The next
number however is different, showing that the PRNG state is passed along by
the State monad when randDouble is called for the second time.

Since State is a functor, map can be used to transform the output of the
function randDouble in order to build other random generators:

⬇

randDouble.map(x => (x * (to - from) + from).toInt)

def randBool: Rand[Boolean] =

randDouble.map(_ > 0.5)

The State monad can be used to define a Rand monad using other algorithms
which generate uniformly distributed random variables, and these can be
transformed appropriately to generate values distributed according to other
useful distributions. For example the Box-Muller transform is an efficient
method of generating standard Normal random numbers (Box and Muller, 1958).

##  4 Probability Monads

In order to define a probabilistic program using the abstractions introduced
so far, it remains to show that probability measures form a monad. Monads are
well understood and ubiquitous in functional programming and programming
language theory, hence have well supported language features such as the
`for`-comprehension in Scala and `do`-notation in Haskell. Their uniformity
enables programmers to manipulate monadic values in a consistent way using map
and flatMap independent of the monadic context. This makes defining
probabilistic programs straightforward using powerful, well-developed
functional programming languages with support for monads.

The foundations of the probability monad were first developed by Lawvere
(1962), later extended in Giry (1982). A review of the Giry monad is presented
in Ramsey and Pfeffer (2002). The Giry monad is defined on the category of
measurable spaces, Meas. The Giry monad maps each measurable space $X$ to the
space of of all probability measures on $X$, which we will denote $P(X)$. The
monad is characterised by the two natural transformations $\eta$ and $\mu$.
The unit $\eta_{X}:X\rightarrow P(X)$ must map each $x\in X$ to a probability
measure on $X$, and does so in the obvious way by mapping to the Dirac
measure, $\eta_{X}(x)=\delta_{x}$. The multiplication
$\mu_{X}:P(P(X))\rightarrow P(X)$ must flatten a probability measure over
probability measures on $X$, down to a probability measure on $X$, and does so
in the obvious way, by marginalisation (via Lebesgue integration). Then the
monadic bind operation (`flatMap` in Scala) can be interpreted as the law of
total probability. Explicitly, for probability kernels (conditional
distributions) $f:X\rightarrow P(Y)$ and $g:Y\rightarrow P(Z)$ we bind them as
$g\bullet f:X\rightarrow P(Z)$ via

$g\bullet f=\int_{Y}g(y)\,\textrm{d}f(x).$

One issue with the Giry monad is that it is defined on Meas, which is not
Cartesian-closed. This is not ideal for any kind of higher-order probabilistic
programming language, since it prevents convenient use of probability
distributions over functions. It is also awkward for formalising the semantics
of monadic probabilistic programming languages developed as embedded DSLs in
functional programming languages with support for monads, since there will
typically be an implicit assumption that all monads are defined over CCCs.

Measurable spaces are very rich, and include many pathological cases which are
not of significant practical interest, yet cause substantial technical
difficulties. So one approach to “fixing” the problem is to restrict attention
to spaces of more direct relevance, in a principled way. This has led to the
development of a probability monad based on Quasi-Borel spaces, which form a
CCC, as explored in Heunen et al. (2017). Other approaches make a more radical
departure from conventional measure-theoretic probability, such as the
Kantorovich monad, defined on metric spaces (Fritz and Perrone, 2017; Jacobs,
2018). The use of probability monads as a foundation for Bayesian inference is
explored in Culbertson and Sturtz (2014).

This theoretical underpinning paves the way for developing monadic
probabilistic programming DSLs in functional programming languages with
minimal effort beyond specifying the two natural transformations required for
the distribution to form a monad (see section 3.8). Although the formal
connection with the theory of probability monads is important for establishing
rigorous semantics for a given probabilistic programming language (Staton,
2017; Ścibior et al., 2018b), the practical implications associated with
picking a specific probability monad formalism for the development of actual
implementations are limited (Ścibior et al., 2015, 2018a).

##  5 Hamiltonian Monte Carlo

One of the aims of probabilistic programming is to separate the details of
inference algorithms from the specification of statistical models. This allows
for rapid exploration of different models without having to devise new
inference algorithms. Metropolis-Hastings (MH) is an MCMC algorithm with
convergence guarantees and if run for long enough will give samples from the
stationary distribution corresponding to the posterior distribution of
interest in a Bayesian statistical model (Metropolis et al., 1953; Hastings,
1970). The MH algorithm uses a parameter proposal distribution which requires
tuning in order to achieve the optimal acceptance rate of around 0.234
(Roberts et al., 1997). Adaptive methods can be used to tune the proposal
distribution automatically (Atchadé et al., 2005). Even optimal tuning of the
proposal distribution does not take into account the geometry of the posterior
distribution, which is contained in the gradient of the posterior. A more
efficient inference algorithm is desirable for general probabilistic
programming and Bayesian inference.

Hamiltonian Monte Carlo utilises the gradient of the un-normalised log-
posterior to more efficiently explore the posterior distribution. The
intuition for HMC is developed from Hamiltonian dynamics. Hamilton’s equations
are written as:

$\displaystyle\frac{\mathrm{d}\boldsymbol{p}}{\mathrm{d}t}$
$\displaystyle=-\frac{\partial\mathcal{H}}{\partial\boldsymbol{q}},$
$\displaystyle\frac{\mathrm{d}\boldsymbol{q}}{\mathrm{d}t}$
$\displaystyle=+\frac{\partial\mathcal{H}}{\partial\boldsymbol{p}},$

where $\boldsymbol{q}$ represents a particle position and $\boldsymbol{p}$ is
the momentum of the particle. The Hamiltonian of a physical system can be
written as the sum of the kinetic and potential energy:

$\mathcal{H}(\boldsymbol{q},\boldsymbol{p})=T(\boldsymbol{q})+V(\boldsymbol{p})$
(2)

The HMC algorithm uses a combination of Gibbs sampling, Hamiltonian Dynamics
and a Metropolis-Hastings step; hence, it is sometimes called Hybrid Monte
Carlo. First the posterior distribution is augmented with an additional
momentum parameter $\phi$. This is an auxiliary parameter which is not of
direct interest when calculating the parameter posterior distribution. The
parameters, $\psi$ correspond to the position in Hamilton’s equations. We take
the kinetic energy to be $T(\phi)=\frac{1}{2}\phi^{T}\phi$, assuming a unit
particle mass. The joint density of the position and momentum can then be
written as:

$p(\psi,\phi)\propto\exp\left\\{\log
p(\psi|y)-\frac{1}{2}\phi^{T}\phi\right\\},$ (3)

where $\log p(\psi|y)$ is the log of the target posterior distribution written
up to a normalising constant and represents the negative potential energy. The
kinetic energy is the kernel of a standard multivariate Normal distribution,
and the identity covariance matrix can replaced by a tuning parameter,
$\Sigma$ termed the mass matrix, from analytical mechanics (Betancourt, 2017).
Hamilton’s equations are discretised in order to update the values of the
static parameters, $\psi$ (the position in Hamilton’s equations) and the
momentum $\phi$. A special discretisation of Hamilton’s equations is used
called a leapfrog step:

$\displaystyle\phi_{t+\varepsilon/2}$
$\displaystyle=\phi_{t}+\frac{\varepsilon}{2}\nabla_{\psi}\log p(\psi_{t}|y),$
$\displaystyle\psi_{t+\varepsilon}$
$\displaystyle=\psi_{t}+\varepsilon\phi_{t+\varepsilon/2},$
$\displaystyle\phi_{t+\varepsilon}$
$\displaystyle=\phi_{t+\varepsilon/2}+\frac{\varepsilon}{2}\nabla_{\psi}\log
p(\psi_{t+\varepsilon}|y).$ (4)

$\nabla_{\psi}\log p(\psi|y)$ is the gradient of the un-normalised log-
posterior distribution with respect to the parameters, $\psi$. $\varepsilon$
is a tuning parameter in the HMC algorithm and represents the step-size of a
leapfrog step. This leapfrog update is more accurate than a naive Euler
discretisation, primarily due to the fact that it is volume-preserving (since
each step is a shear). The volume-preservation is important for HMC, since it
facilitates reversibility, and avoids the need to track Jacobians. The steps
required to perform HMC are summarised in Algorithm 1. The function
$\operatorname{leapfrogs}$ is a recursive function which applies
$\operatorname{leapfrog\\_step}$ $L$ times without mutating state and returns
the updated position and momentum.

Result: Return $\psi_{i},i=1,\dots,M$

1 Given the prior distribution $p(\psi)$;

2 The un-normalised the log-posterior distribution $\log p(\psi|y)$ and its
derivative $\nabla_{\psi}\log p(\psi|y)$;

3 The tuning parameters $\varepsilon$, $L$;

4 Initialise the parameters by sampling $\psi_{0}\sim p(\psi)$;

5 for  _i in 1 to M_ do

6 Sample $\phi\sim\textrm{MVN}(0,\Sigma)$.;

7
$\psi^{\star},\phi^{\star}=\operatorname{leapfrogs}(\psi^{\star},\phi^{\star},L,\varepsilon)$;

8 Calculate $\alpha=\log p(\psi^{\star}(L)|y)-\frac{1}{2}\phi^{\star
T}(L)\Sigma^{-1}\phi^{\star}(L)-\log
p(\psi_{i-1}|y)+\frac{1}{2}\phi^{T}\Sigma^{-1}\phi$;

9 Sample $u\sim U(0,1)$;

10 if  _$\log(u) <\alpha$_ then

11 Set $\psi_{i}:=\psi^{\star}$

12 else

13 Set $\psi_{i}:=\psi_{i-1}$

14

15 end for

16 function  _leapfrog_step(_$\psi$ , $\phi$, $\varepsilon$_)_:

17 $\tilde{\phi}=\phi+\frac{\varepsilon}{2}\nabla_{\psi}\log p(\psi|y)$;

18 $\tilde{\psi}=\psi+\varepsilon\Sigma^{-1}\tilde{\phi}$;

_19_ _ _ return
_$\tilde{\psi},\tilde{\phi}+\frac{\varepsilon}{2}\nabla_{\psi}\log
p(\tilde{\psi}|y)$ ;_

20

21

22 function  _leapfrogs(_$\psi$ , $\phi$, $\varepsilon$, $L$_)_:

23 if  _$L=0$_ then

24 return _$\psi,\phi$_

25 else

26
$\tilde{\psi},\tilde{\phi}=\operatorname{leapfrog\\_step}(\psi,\phi,\varepsilon)$;

27 $\operatorname{leapfrogs}(\tilde{\psi},\tilde{\phi},\varepsilon,L-1)$

28

29

Algorithm 1 Hamiltonian Monte Carlo

For differentiable targets, HMC can be shown to exactly preserve the required
distribution, and is typically more efficient then random walk Metropolis-
Hastings schemes. The optimal acceptance rate for HMC is approximately 0.65
(Neal et al., 2011). Selecting the tuning parameters corresponding to the
leapfrog step size $\varepsilon$ and the number of leapfrog steps $L$ is
typically done using short pilot runs of the chain and targeting an acceptance
rate of 0.65. Efforts have been made to automate selection of tuning
parameters in the HMC algorithm and have led to the No-U-turn sampler which
performs the optimum number of leapfrog steps (Hoffman and Gelman, 2014) and
empirical HMC (Wu et al., 2018). In both cases the optimum number of leapfrog
steps are chosen such that the discretised steps in the posterior distribution
does not make an U-turn by first heading away from the previously accepted
parameter value, then turning back in the direction of the previous parameter
when the value of the gradient changes.

##  6 Automatic Differentiation

The need to calculate gradients in many inference and optimisation algorithms
such as Hamiltonian Monte Carlo (HMC) and variational inference (Kucukelbir et
al., 2017) has led to a renewed interest in automatic differentiation (AD). AD
is a way of calculating derivatives of functions whilst at the same time
evaluating them. This is not numerical differentiation or symbolic
differentiation but rather exact differentiation which returns the value of a
derivative at a point (Wang et al., 2018). It is straightforward to implement
forward mode automatic differentiation in an FP language. However the number
of computations performed using forward mode AD depends on the dimension of
the input space, and hence does not scale well to parameter inference in
models involving a large number of parameters. On the other hand the number of
computations required for reverse mode automatic differentiation scales with
the dimension of the output dimension. In Bayesian inference and HMC in
particular the gradient of the log-posterior with respect to the free-
parameters is required, this is a function from
$\mathbb{R}^{n}\rightarrow\mathbb{R}$, where $n$ is the dimension of the
parameter space, hence reverse mode AD is typically more efficient.

Hamiltonian Monte Carlo was introduced in section 5, along with some
implementation difficulties, such as choosing the leapfrog step size, number
of leapfrog steps and deriving the partial derivatives required for the
proposal. Automatic differentiation can be used to calculate the exact
derivatives needed when performing the leapfrog step of HMC.

Dual numbers can be used in order to calculate derivatives (of univariate
functions) automatically and exactly. Each real number has a corresponding
dual number, which is the number, $y\in\mathbb{R}$ plus a small innovation,
$\varepsilon$ such that $\varepsilon^{2}=0$.

Derivatives can be calculated by evaluating functions using the dual number,
for instance the function $f:\mathbb{R}\rightarrow\mathbb{R}$ defined by
$f(x)=x^{2}+2x+5$, the derivative is $f^{\prime}(x)=2x+2$. In order to
calculate the derivative automatically using dual numbers, the function is
evaluated using the dual number equivalent to a chosen value of $x$ for
instance $x=5$, has the dual number $5+\varepsilon$ then

$\displaystyle f(5+\varepsilon)$
$\displaystyle=(5+\varepsilon)^{2}+2(5+\varepsilon)+5,$
$\displaystyle=25+10\varepsilon+\varepsilon^{2}+10+2\varepsilon+5,$
$\displaystyle=40+12\varepsilon.$

In this way, the evaluation of $f$ using the dual number has resulted in the
evaluation of $f(x)$, and $f^{\prime}(x)$ simultaneously with $f^{\prime}(x)$
being given by the coefficient of $\varepsilon$. This is the essence of
forward mode automatic differentiation.

Automatic differentiation using dual numbers is equivalent to applications of
the chain rule:

$(f\circ g)^{\prime}=(f^{\prime}\circ g)g^{\prime}$ (5)

Combinations of primitive functions can be differentiated using repeated
applications of the chain rule. flatMap can be defined for the Dual class
which encapsulates the chain rule. This presentation of forward mode automatic
differentiation was first outlined by (Welsh, 2018).

⬇

def flatMap[B: Numeric](f: A => Dual[B]): Dual[B] =

val x = f(real)

val nextEps = eps * x.eps

Dual(x.real, nextEps)

}

The Numeric context bound is required since the primitive functions, such as
multiplication, are defined explicitly in terms of their real result and
derivative; hence an eps: Double has to be multiplied by a real: A when
applying the product rule for differentiation:

⬇

Dual[A](x.real * y.real,

x.eps * ev.toType[Double](y.real) +

y.eps * ev.toType[Double](x.real))

The derivatives of special functions can be written as:

⬇

Dual(math.sin(real), math.cos(real) * eps)

def log[A: Numeric](a: Dual[A])(implicit ev: A =:= Double) =

Dual(math.log(a.real), a.eps / a.real)

The derivatives of the special functions can also be written using flatMap

⬇

a.flatMap(x => Dual(math.sin(x), math.cos(x)))

def log(a: Dual[A])(implicit ev: A =:= Double) =

a.flatMap(x => Dual(math.log(real), 1 / real))

The implementation becomes more complex as functions of multiple arguments are
considered, since a separate $\varepsilon$ is required for the derivative with
respect to each argument (Manzyuk, 2012).

Figure 3 is a representation of a function with multiple arguments
$f(x_{1},x_{2})=x_{1}x_{2}+x_{1}^{2}$. The arrows represent function
application, intermediate nodes represent the result of each primitive
function and are labelled $x_{3}$ and $x_{4}$.

Figure 3: Computational graph showing the decomposition of
$f(x_{1},x_{2})=x_{1}x_{2}+x_{1}^{2}$ into primitive functions

Differentiating the function with respect to all variables in the graph using
forward mode AD requires multiple traversals of the graph. The derivative of a
node, $x_{k}$ with parent nodes $Pa(x_{k})$ with respect to the input variable
$x_{1}$ is derived from the chain rule

$\frac{\partial x_{k}}{\partial x_{1}}=\sum_{j\in Pa(x_{k})}\frac{\partial
x_{k}}{\partial x_{j}}\frac{\partial x_{j}}{\partial x_{1}}.$ (6)

Figure 4 shows the forward pass required to calculate the gradient
$\frac{\partial f}{\partial x_{1}}\rvert_{x_{1}=1,x_{2}=2}$. This single
forward pass evaluates the function and the derivative with respect to $x_{1}$
given the known derivatives of primitive functions such as product and power.
A second pass has to be computed in order to determine the derivative with
respect to $x_{2}$. The gradient calculated by hand is, $\frac{\partial
f}{\partial x_{1}}=x_{2}+2x_{1}$.

Figure 4: Computational graph showing the applications of the chain rule
required to calculate the partial derivative $\frac{\partial f}{\partial
x_{1}}$.

A method often used to implement AD for functions of multiple arguments
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is to track a vector of derivatives.

A minimal implementation of forward mode AD along with a test suite verifying
the algorithm using property based testing and a selection of different
functions is available at 333https://git.io/probprog.

###  6.1 Reverse Mode Automatic Differentiation

Reverse mode automatic differentiation is typically faster than forward mode
AD when functions have a larger input space than output space, ie.
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ where $n>m$. This is true for
Bayesian inference algorithms which require the gradient of the un-normalised
log-posterior such as Hamiltonian Monte Carlo or Metropolis-adjusted Langevin
algorithm (Roberts and Rosenthal, 1998). The un-normalised log-posterior is a
function from the parameters of dimension $n$ to a single real number.

In order to differentiate a function using reverse mode automatic
differentiation the steps are similar to forward mode differentiation, the
derivatives of simple functions are defined. The function is then decomposed
into its constituent primitive functions. The derivative of with respect to
each node of the adjoint graph is calculated at $x_{1}=2$, $x_{2}=3$ in figure
5.

Figure 5: Graph with derivatives of each input variable calculated on the
nodes of the adjoint graph

Then proceeding in reverse through the graph to calculate the derivative of
the function with respect to the input variables:

$\frac{\partial f}{\partial x_{k}}=\sum_{j\in Ch(x_{k})}\frac{\partial
x_{j}}{\partial x_{k}}\frac{\partial f}{\partial x_{j}}$

Where $Ch(x_{k})$ represent the set of children of node $x_{k}$. Figure 6
shows the reverse sweep through the computation graph required to calculate
the derivative with respect to all of the arguments of the function.

Figure 6: Reverse pass through the graph to calculate the derivative of $f$
with respect to all the input variables

Automatic differentiation can be a part of a monadic probabilistic program
which allows distributions to be composed into hierarchical models. Reverse
mode automatic differentiation can be implemented as a monad using
continuations 444https://na.scaladays.org/schedule/differentiable-functional-
programming.

##  7 Architecture of a Probabilistic Program

The Scala library Rainier (Bryant, 2018) is an example of a monadic
probabilistic programming language. It can be used to build a compute graph
over the parameters of a model. The compute graph is essentially a variadic
function from many input parameters to an output. This function is a
representation of the un-normalised log-posterior function and returns both
the evaluation of the function and the gradient. The compute graph can be
compiled to a function with signature Array[Double] => Double to evaluate the
un-normalised log-posterior and Array[Double] => Array[Double] for the
corresponding gradient.

The casual user is not aware of the compute graph and simply constructs a
model by combining primitive probabilistic functions. A distribution is
defined as a trait:

⬇

def logPdf(a: A): Double

def gen: Rand[A]

def param: RandomVariable[Real]

def fit(ys: A*): RandomVariable[A]

}

Listing 2: Abstract interfaces for a continuous distribution

A continuous distribution can evaluate the log probability density function,
logPdf and be transformed into a RandomVariable monad for use as a parameter
in a composed model using param. Discrete distributions can be defined
similarly, except the log probability mass function is defined: logPmf and
discrete distributions can not be used as a parameter in HMC since it can not
be differentiated and hence the param method is not defined for discrete
distributions. The fit method accepts as an argument a single observation, or
multiple and allows the distribution to be used as the likelihood for the
observation. Each interface has a method, gen which can transform the
distribution to a sampling based distribution, called a Generator in Rainier
and sometimes referred to as Rand. Rand is also a monad:

⬇

def draw(): A

def map[B](f: A => B): Rand[B] = new Rand[B] {

override def draw = f(self.draw)

}

def flatMap[B](f: A => Rand[B]): Rand[B] = new Rand[B] {

override def draw = f(self.draw).draw

}

}

and has similar semantics to the Rand monad discussed earlier in the context
of pure functional random number generation. Then the only function which
remains abstract and must be defined for a new Rand is draw, hence this monad
has no information about the logPmf or logPdf functions. The functions flatMap
and pure are the two functions required to define a monad, if the natural
transformations defined on Rand satisfy the monad laws then Rand is a monad.

A model to infer the posterior distribution for the probability of heads on a
coin from ten coin flips can be written as

⬇

p <- Beta(3, 3).param

_ <- Binomial(p, 10).fit(6)

} yield p

Listing 3: Rainier probabilistic program expressing a coin flip example with a
beta prior distribution for the probability of heads

Using .param to indicate the bounded parameter $p$ has a beta prior
distribution with bounded support on $[0,1]$, then the discrete binomial
distribution is used as the likelihood by calling the method fit with the
total trials $n=10$ and observed number of heads $y=6$. Each model defined
using this monadic syntax can be sampled from using a Sampler. The Sampler
depends on the compute graph which can derive the un-normalised log-posterior
and its gradient.

##  8 Example Probabilistic Programs

This section contains some example programs written using the Rainier DSL. The
models are intentionally simple, as they are intended to present just the most
basic concepts. Some areas in which assumptions of these simple models can be
relaxed are highlighted — in order to improve model fit. This emphasises the
flexibility of probabilistic programming.

For each of the examples the HMC algorithm with dual averaging and $L=5$
leapfrog steps was used with 10,000 iterations for the warm-up run during
which the initial leapfrog step size is determined followed by 50,000 sampling
iterations with a thinning factor of 5.

###  8.1 Linear Model

Linear models are statistical models specifying a relationship between
covariates, $X\in\mathbb{R}^{n\times p}$ and a univariate outcome
$Y\in\mathbb{R}^{n}$ for each example via the coefficient $\beta$, a
$(p+1)$-vector. Each row of the outcome matrix $Y$ is denoted as $y_{i}$ and
is related to each row of the covariate matrix, $x_{i}$ by the coefficient
matrix $\beta$. The covariate matrix has a column vector of ones prepended as
the first column which represents the intercept. The observations have
independent normally distributed noise with equal variance:

$\displaystyle y_{i}$
$\displaystyle=\beta^{T}x_{i}+\varepsilon_{i},\qquad\varepsilon_{i}\sim\mathcal{N}(0,\sigma),$
(7) $\displaystyle\beta$
$\displaystyle\sim\textrm{MVN}(\mu_{\beta},\Sigma_{\beta}^{2}),$
$\displaystyle\sigma$
$\displaystyle\sim\textrm{Exponential}(\lambda_{\sigma}).$

The model assumptions can be relaxed, but then model interpretation can become
more challenging. 1,000 observations from this model are simulated with
$\beta_{0}=4.0,\beta_{1}=-1.5,\sigma=0.5$ and the covariates simulated from
the standard Normal distribution. The bivariate relationships of the simulated
data is plotted in Figure 7 (a). The figures were produced using R (R Core
Team, 2019) and ggplot2 (Wickham, 2016).

(a)  (b)

Figure 7: (a) Simulated data from a linear model with
$\beta_{0}=4.0,\beta_{1}=-1.5,\sigma=0.5$ and covariates simulated
independently from a standard Normal distribution. (b) Parameter posterior
diagnostic plots for the linear model using the HMC kernel and simulated data,
the actual values of the parameters are plotted using the dashed lines (Top)
Traceplots (Middle) Autocorrelation plots (Bottom) Marginal Densities

Listing 4 shows a probabilistic program for a linear model with covariates x
and outcome y. Each parameter is added to the model using the .param notation
and the $p=1$ covariates are related to each scalar outcome using
Predictor.from. The linear model is expressed as a function which can be re-
used in more complex hierarchical models, such as the random effects model
presented in section 8.3.

⬇

1 sigma: Real, data: Vector[(Double, Double)]) = for {

2 alpha <- Normal(pa, pas).param

3 beta <- Normal(pb, pbs).param

4 _ <- Predictor[Double].from { x =>

5 Normal(alpha + beta * x, sigma)

6 }

7 .fit(data)

8} yield Map("alpha" -> alpha, "beta" -> beta, "sigma" -> sigma)

9

10val model = for {

11 sigma <- Exponential(3.0).param

12 params <- linearModel(0.0, 0.01, 0.0, 0.01, sigma, data)

13} yield params

Listing 4: Simple linear model with one covariate and an intercept with a
Normal observation model and fully specified prior distributions

Figure 7 (b) shows posterior inferences using simulated data with
$\beta_{0}=4.0,\beta_{1}=-1.5,\sigma=0.5$ and the covariates simulated from a
standard Normal distribution.

It is straightforward to change the observation distribution of this model,
for instance to a generalised linear model with a Poisson observation
distribution by changing line 7 to:

⬇

Poisson((alpha + beta * x).exp)

In addition the standard deviation parameter sigma is no longer required when
using a Poisson observation distribution. Any suitable prior distributions can
be specified for the parameters; they do not have to be conditionally
conjugate as in Gibbs sampling.

###  8.2 Mixture Model

In a mixture model the observed data is assumed to arise from a finite mixture
of independent distributions. The mixture model considered here is a mixture
of $m=3$ Normal distributions with different mean values,
$\mu_{k},k=1,\dots,m$ and a common variance.

$\displaystyle y_{i}$ $\displaystyle\sim\mathcal{N}(\mu_{k},\sigma),$
$\displaystyle k$ $\displaystyle\sim\mathcal{F}(\theta),$
$\displaystyle\theta$
$\displaystyle\sim\textrm{Dir}(\boldsymbol{\alpha}_{\theta}),$
$\displaystyle\mu_{k}$
$\displaystyle\sim\mathcal{N}(\mu_{\mu_{k}},\sigma^{2}_{\mu_{k}}),$
$\displaystyle\sigma$
$\displaystyle\sim\textrm{Exponential}(\alpha_{\sigma},\beta_{\sigma}).$

where each index $k$ is drawn from a discrete distribution $\mathcal{F}$ with
probabilities $\theta$. The likelihood of the mixture model can be written as:

$p(y|\theta,\mu,\sigma)=\sum_{k=1}^{m}\theta_{k}\mathcal{N}(y;\mu_{k},\sigma).$
(8)

This mixture distribution has a smooth log-density and hence can be
differentiated and inference can be performed using HMC. In addition the
mixture distribution is already implemented in Rainier. The mixture model can
be defined as:

⬇

val total = alphas.reduce(_ + _)

alphas.map(a => a / total)

}

val model = for {

unnormThetas <- RandomVariable.traverse(Vector.fill(3)(Gamma(3.0, 1.0).param))

thetas = normalise(unnormThetas)

mus <- RandomVariable.traverse(Vector.fill(3)(Normal(0, 1).param))

sigma <- Exponential(3.0).param

dists = mus.

map(mu => Normal(mu, sigma)

components: Map[Continuous, Real] = dists.zip(thetas).toMap

_ <- Mixture(components).fit(ys)

} yield thetas ++ mus ++ sigma

Listing 5: Rainier model for a Gaussian mixture model with $m=3$ components.

The probability of each mixing component is $\theta_{i}$, $i=1,2,3$. The
mixing components must be greater than zero and sum-to-one, they are drawn
from a Dirichlet distribution by drawing each component from a Gamma
distribution with scale, $\theta=1$ then normalising the values. Figure 8 (a)
shows a simulation from the mixture model and (b) shows the posterior
diagnostics for the means and mixing components.

(a)  (b)

Figure 8: (a) 10,000 simulations from a three component Normal mixture model
with mean values $\mu=\\{-2.0,1.0,3.0\\}$, common variance $\sigma=0.5$ and
proportions $\theta=\\{0.3,0.2,0.5\\}$ (b) Parameter posterior marginal
density plots for the mixture model parameters, the parameter values used to
simulate the data are vertical dashed lines

###  8.3 Random Effects Model

Hierarchical models give an opportunity to emphasise the manipulation of
probabilistic programs as values in the embedded DSL. This model is a random
effects model with $k$ classes each with $n$ observations. Each class is
modelled using an independent linear regression.

$\displaystyle Y_{ij}$
$\displaystyle=\alpha_{i}+\beta_{i}x_{j}+\varepsilon_{ij},$
$\displaystyle\varepsilon_{ij}\sim\mathcal{N}(0,\sigma^{2}),$
$\displaystyle\alpha_{i}$
$\displaystyle\sim\mathcal{N}(\alpha_{c},\sigma_{a}^{2}),$
$\displaystyle\beta_{i}$
$\displaystyle\sim\mathcal{N}(\beta_{c},\sigma_{b}^{2}).$

Each of the $i=1,\dots,k$ classes have the same covariates $x_{j},j=1,\dots,n$
and have a linear relationship with independent coefficients, $\alpha_{i}$ and
$\beta_{i}$. The standard deviation of the observation, $\sigma$ is assumed to
be the same for each class and observation. The coefficients of the covariates
for all models are assumed to be Normally distributed with the same mean and
variance, this induces correlation between the $k$ linear models.
$\alpha_{c},\sigma_{a},\beta_{c},\sigma_{b}$ and $\sigma$ are given the
following prior distributions:

$\displaystyle\alpha_{c}$ $\displaystyle\sim\mathcal{N}(0,10^{6}),$
$\displaystyle\sigma_{a}$ $\displaystyle\sim\textrm{Gamma}(10^{-3},10^{3}),$
$\displaystyle\beta_{c}$ $\displaystyle\sim\mathcal{N}(0,10^{6}),$
$\displaystyle\sigma_{b}$ $\displaystyle\sim\textrm{Gamma}(10^{-3},10^{3}),$
$\displaystyle\sigma$ $\displaystyle\sim\textrm{Exponential}(3).$

The Gamma distribution is parameterised here using shape, $k$ and scale,
$\theta$ such that the mean is $k\theta$. Each observation is represented by a
class containing the index of the class, the covariate, $x$, and the
observation, $y$:

⬇

id: Int,

x: Double,

y: Double

)

Then the prior distributions are translated to Rainier, note that the Gamma
distribution is parameterised in terms of shape, $k$ and scale, $\theta$ and
the Normal distribution by its mean, $\mu$ and standard deviation $\sigma$.

⬇

alphaC <- Normal(0.0, 1000).param

alphaSigma <- Gamma(0.001, 1000).param

betaC <- Normal(0.0, 1000).param

betaSigma <- Gamma(0.001, 1000).param

sigma <- Exponential(3).param

} yield (alphaC, alphaSigma, betaC, betaSigma, sigma)

A function to fit a single linear regression for the $i^{\textrm{th}}$ class
can be re-used from the function linearModel presented in section 8.1.

linearModel returns a RandomVariable monad which is its own probabilistic
program, but this program only specifies the model for a single class. The
collection of observations, x and y for a single class can be grouped by the
id of the class using groupBy and a linear regression can be fit for each
class using map. When the function groupBy is called on the collection of
observations, it returns a Map[Int, Iterable[Observation]], the key of the map
is an integer representing the id of each class and the Iterable in the value
of the Map corresponds to the observations for each class. The application of
linearModel to each element of the Map returns an
Iterable[RandomVariable[(Double, Double)] which is converted to a Vector
before applying the function RandomVariable.traverse. traverse “reverses” the
order of the effects to return a RandomVariable[Vector[(Double, Double)]]
which is a model and can be sampled from.

⬇

(ac, sa, bc, sb, sigma) <- prior

_ <- RandomVariable.traverse(

observations.

groupBy(_.id).

map { case (id, obs) =>

linearModel(ac, sa, bc, sb, sigma, obs.map(ys => (ys.x, ys.y)))

}.toVector)

} yield (ac, sa, bc, sb, sigma)

This implementation emphasises the compositionality of probabilistic programs
embedded in a host language. The higher order functions, groupBy and map are
familiar to any functional programmer and can be used to combine simple models
into complex hierarchical models.

The diagnostics of the draws from the posterior distribution obtained using
HMC on the random effects hierarchical example are presented in figure 9. The
actual values used to simulate the data are plotted using dashed lines.

Figure 9: Diagnostic plots for the hierarchical random effects model with the
values used to simulate the overlaid using dashed lines (Top) Traceplots
(Middle) Autocorrelation (Bottom) Empirical marginal densities

##  9 Conclusion

Probabilistic programming aims to unify Bayesian inference with general
purpose programming languages in order to simplify model exploration and
inference. This allows practitioners to develop novel statistical models with
minimal consideration for the underlying inference algorithms. In addition,
multiple inference algorithms can be compared without re-writing the model
code.

It has been shown that functional programming in Scala is a suitable, powerful
language for developing a probabilistic programming language as an embedded
DSL. Static typing reduces difficult to find runtime errors and ensures that
refactoring code is straightforward with the help of the compiler. This means
code is not over-engineered from the start of project as changes can be made
to any size code base as and when the specification changes. Static typing
does not increase the size of the code base since well defined abstractions
from category theory allow for polymorphic code with type safety.
Additionally, Scala has local type inference, allowing for return types which
are obvious to be omitted from the source code and inferred by the compiler,
which is especially helpful when defining anonymous functions.

In order to maximise the utility of the statistical models, the DSL is
embedded in an existing general purpose programming language. Functional
programming languages have powerful, theoretically sound abstractions such as
the free monad and tagless final encodings which can be used to develop
embedded DSLs for probabilistic programming. This means the same model code is
used in production as in development, reducing the possibility of bugs in
production code leading to incorrect decisions being made. The Rainier DSL
provides a powerful modelling language embedded in the Scala language. The
‘for’-comprehension used to build models from probability distributions is
familiar to Scala programmers since many common data constructors form a
monad. The monadic syntax provided by ‘for’ is reminiscent of the syntax in
BUGS, JAGS and Stan allowing experts in these languages to easily transition
to using the embedded DSL in Scala. Embedding the probabilistic programming
DSL in a host functional language not only renders trivial the problem of
“compiling” the language to an executable sampling algorithm, but also allows
probabilistic programs to be manipulated as regular values in the host
language. This maximises the opportunities for model composition, modular
model development and code reuse. This is exemplified in the random effects
example in section 8.3 .

The Rainier DSL provides a powerful modelling language embedded in the Scala
language. The for comprehension used to build models from probability
distributions is familiar to Scala programmers since many common data
constructors form a monad. The monadic syntax provided by for is reminiscent
of the syntax in BUGS, Jags and Stan allowing experts in these languages to
easily transition to using the embedded DSL in Scala. This syntax is available
with minimal effort when a monad is defined by specifying the flatMap and pure
functions.

## Acknowledgements

JL is supported by the Engineering and Physical Sciences Research Council,
Centre for Doctoral Training in Cloud Computing for Big Data (grant number
EP/L015358/1) and Digital Catapult Teaching Grant Award (KH153326). DJW would
like to thank the Isaac Newton Institute for Mathematical Sciences for support
and hospitality during the programme Statistical Scalability when some work on
this paper was undertaken; EPSRC grant EP/R014604/1. This work was also
supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1.

## References

  * Abadi et al. (2015) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Apache (2018a) Apache.  Apache Hadoop, 2018a.  URL https://hadoop.apache.org/.
  * Apache (2018b) Apache.  Apache Spark - Unified Analytics Engine for Big Data, 2018b.  URL http://spark.apache.org.
  * Atchadé et al. (2005) Y. F. Atchadé, J. S. Rosenthal, et al.  On adaptive Markov chain Monte Carlo algorithms.  _Bernoulli_ , 11(5):815–828, 2005.
  * Awodey (2010) S. Awodey.  _Category theory_.  Oxford University Press, New York, second edition, 2010.
  * Barr and Wells (1990) M. Barr and C. Wells.  _Category theory for computing science_ , volume 49.  Prentice Hall New York, 1990.
  * Betancourt (2017) M. Betancourt.  A conceptual introduction to Hamiltonian Monte Carlo.  _arXiv preprint arXiv:1701.02434_ , 2017.
  * Box and Muller (1958) G. E. P. Box and M. E. Muller.  A note on the generation of random normal deviates.  _Ann. Math. Statist._ , 29(2):610–611, 06 1958\.  doi: 10.1214/aoms/1177706645.  URL https://doi.org/10.1214/aoms/1177706645.
  * Brin and Page (1998) S. Brin and L. Page.  The anatomy of a large-scale hypertextual web search engine.  _Computer networks and ISDN systems_ , 30(1-7):107–117, 1998.
  * Bryant (2018) A. Bryant.  Rainier: Bayesian inference in Scala, 2018.  URL https://github.com/stripe/rainier/.
  * Carpenter et al. (2016) B. Carpenter, A. Gelman, M. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. A. Brubaker, J. Guo, P. Li, A. Riddell, et al.  Stan: A probabilistic programming language.  _Journal of Statistical Software_ , 20(2):1–37, 2016.
  * Chiusano and Bjarnason (2014) P. Chiusano and R. Bjarnason.  _Functional programming in Scala_.  Manning Publications Co., 2014.
  * Culbertson and Sturtz (2014) J. Culbertson and K. Sturtz.  A categorical foundation for Bayesian probability.  _Applied Categorical Structures_ , 22(4):647–662, 2014.
  * Dean and Ghemawat (2008) J. Dean and S. Ghemawat.  MapReduce: simplified data processing on large clusters.  _Communications of the ACM_ , 51(1):107–113, 2008\.
  * Duane et al. (1987) S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth.  Hybrid Monte Carlo.  _Physics letters B_ , 195(2):216–222, 1987.
  * Fredkin (1960) E. Fredkin.  Trie memory.  _Commun. ACM_ , 3(9):490–499, Sept. 1960.  ISSN 0001-0782.  doi: 10.1145/367390.367400.  URL http://doi.acm.org/10.1145/367390.367400.
  * Fritz and Perrone (2017) T. Fritz and P. Perrone.  A probability monad as the colimit of spaces of finite samples.  Technical Report arXiv:1712.05363, arXiv, 2017.
  * Giry (1982) M. Giry.  A categorical approach to probability theory.  In _Categorical aspects of topology and analysis_ , pages 68–85. Springer, 1982.
  * Hastings (1970) W. K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.
  * Heunen et al. (2017) C. Heunen, O. Kammar, S. Staton, and H. Yang.  A convenient category for higher-order probability theory.  In _2017 32nd Annual ACM/IEEE Symposium on Logic in Computer Science (LICS)_ , pages 1–12. IEEE, 2017.
  * Hoffman and Gelman (2014) M. D. Hoffman and A. Gelman.  The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15(1):1593–1623, 2014.
  * Huang et al. (2017) G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger.  Densely connected convolutional networks.  In _CVPR_ , volume 1, page 3, 2017.
  * Jacobs (2018) B. Jacobs.  From probability monads to commutative effectuses.  _Journal of Logical and Algebraic Methods in Programming_ , 94:200–237, January 2018.
  * Jones and Duponcheel (1993) M. P. Jones and L. Duponcheel.  Composing monads.  Technical report, Technical Report YALEU/DCS/RR-1004, Department of Computer Science. Yale, 1993.
  * Khomtchouk et al. (2016) B. B. Khomtchouk, E. Weitz, P. D. Karp, and C. Wahlestedt.  How the strengths of lisp-family languages facilitate building complex and flexible bioinformatics applications.  _Briefings in bioinformatics_ , 19(3):537–543, 2016.
  * Kiselyov et al. (2013) O. Kiselyov, A. Sabry, and C. Swords.  Extensible effects: an alternative to monad transformers.  In _ACM SIGPLAN Notices_ , volume 48, pages 59–70. ACM, 2013.
  * Kucukelbir et al. (2017) A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei.  Automatic differentiation variational inference.  _The Journal of Machine Learning Research_ , 18(1):430–474, 2017.
  * Lawvere (1962) F. W. Lawvere.  The category of probabilistic mappings, 1962.  URL https://ncatlab.org/nlab/files/lawvereprobability1962.pdf.
  * Liang et al. (1995) S. Liang, P. Hudak, and M. Jones.  Monad transformers and modular interpreters.  In _Proceedings of the 22nd ACM SIGPLAN-SIGACT symposium on Principles of programming languages_ , pages 333–343. ACM, 1995.
  * Lunn et al. (2000) D. J. Lunn, A. Thomas, N. Best, and D. Spiegelhalter.  WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility.  _Statistics and computing_ , 10(4):325–337, 2000\.
  * Manzyuk (2012) O. Manzyuk.  A simply typed $\lambda$-calculus of forward automatic differentiation.  _Electronic Notes in Theoretical Computer Science_ , 286:257–272, 2012.
  * Matsumoto and Nishimura (1998) M. Matsumoto and T. Nishimura.  Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator.  _ACM Transactions on Modeling and Computer Simulation (TOMACS)_ , 8(1):3–30, 1998.
  * McBride and Paterson (2008) C. McBride and R. Paterson.  Applicative programming with effects.  _Journal of functional programming_ , 18(1):1–13, 2008.
  * Metropolis et al. (1953) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  _The journal of chemical physics_ , 21(6):1087–1092, 1953.
  * Milewski (2018) B. Milewski.  _Category theory for programmers_.  Blurb, 2018.
  * Narayanan et al. (2016) P. Narayanan, J. Carette, W. Romano, C. Shan, and R. Zinkov.  Probabilistic inference by program transformation in hakaru (system description).  In _International Symposium on Functional and Logic Programming \- 13th International Symposium, FLOPS 2016, Kochi, Japan, March 4-6, 2016, Proceedings_ , pages 62–79. Springer, 2016.  doi: 10.1007/978-3-319-29604-3˙5.  URL http://dx.doi.org/10.1007/978-3-319-29604-3_5.
  * Neal et al. (2011) R. M. Neal et al.  MCMC using Hamiltonian dynamics.  _Handbook of Markov Chain Monte Carlo_ , 2(11), 2011.
  * Odersky et al. (2004) M. Odersky, P. Altherr, V. Cremet, B. Emir, S. Maneth, S. Micheloud, N. Mihaylov, M. Schinz, E. Stenman, and M. Zenger.  An overview of the Scala programming language.  Technical Report IC/2004/64, EPFL Lausanne, Switzerland, 2004.
  * Park et al. (2005) S. Park, F. Pfenning, and S. Thrun.  _A probabilistic language based upon sampling functions_ , volume 40.  ACM, 2005.
  * Paszke et al. (2017) A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in PyTorch.  In _NIPS-W_ , 2017.
  * Plummer et al. (2003) M. Plummer et al.  JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling.  In _Proceedings of the 3rd international workshop on distributed statistical computing_ , volume 124. Vienna, Austria, 2003.
  * R Core Team (2019) R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2019.  URL https://www.R-project.org/.
  * Ramsey and Pfeffer (2002) N. Ramsey and A. Pfeffer.  Stochastic lambda calculus and monads of probability distributions.  In _ACM SIGPLAN Notices_ , volume 37, pages 154–165. ACM, 2002.
  * Roberts and Rosenthal (1998) G. O. Roberts and J. S. Rosenthal.  Optimal scaling of discrete approximations to Langevin diffusions.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 60(1):255–268, 1998.
  * Roberts et al. (1997) G. O. Roberts, A. Gelman, W. R. Gilks, et al.  Weak convergence and optimal scaling of random walk Metropolis algorithms.  _The annals of applied probability_ , 7(1):110–120, 1997.
  * Ścibior et al. (2015) A. Ścibior, Z. Ghahramani, and A. D. Gordon.  Practical probabilistic programming with monads.  In _ACM SIGPLAN Notices_ , volume 50, pages 165–176. ACM, 2015.
  * Ścibior et al. (2018a) A. Ścibior, O. Kammar, and Z. Ghahramani.  Functional programming for modular Bayesian inference.  _Proceedings of the ACM on Programming Languages_ , 2:83, 2018a.
  * Ścibior et al. (2018b) A. Ścibior, O. Kammar, M. Vakar, S. Staton, H. Yang, Y. Cai, K. Ostermann, S. K. Moss, C. Heunen, and Z. Ghahramani.  Denotational validation of higher-order Bayesian inference.  _Proc. ACM Prog. Lang._ , 2(60), 2018b.
  * Staton (2017) S. Staton.  Commutative semantics for probabilistic programming.  In H. Yang, editor, _Programming languages and systems_ , volume 10201 of _ESOP 2017. Lecture notes in computer science_ , Berlin, Heidelberg, 2017. Springer.
  * Swierstra (2008) W. Swierstra.  Data types à la carte.  _Journal of functional programming_ , 18(4):423–436, 2008.
  * TensorFlow (2018) TensorFlow.  Tensorflow probability, 2018.  URL https://www.tensorflow.org/probability/.
  * Uber (2018) Uber.  Pyro: Deep universal probabilistic programming with Python and PyTorch, 2018.  URL https://github.com/uber/pyro.
  * Wadler (1995) P. Wadler.  Monads for functional programming.  In _International School on Advanced Functional Programming_ , pages 24–52. Springer, 1995.
  * Wang et al. (2018) F. Wang, X. Wu, G. Essertel, J. Decker, and T. Rompf.  Demystifying differentiable programming: Shift/reset the penultimate backpropagator.  _arXiv preprint arXiv:1803.10228_ , 2018.
  * Welsh (2018) N. Welsh.  Differentiable Functional Programming, 2018.  URL https://slideslive.ch/38908798/differentiable-functional-programming.
  * Wickham (2016) H. Wickham.  _ggplot2: Elegant Graphics for Data Analysis_.  Springer-Verlag New York, 2016.  ISBN 978-3-319-24277-4.  URL https://ggplot2.tidyverse.org.
  * Wu et al. (2018) C. Wu, J. Stoehr, and C. P. Robert.  Faster Hamiltonian Monte Carlo by Learning Leapfrog Scale.  _arXiv preprint arXiv:1810.04449_ , 2018.