Lightweight Implementations of Probabilistic Programming Languages Via
Transformational Compilation

David Wingate

Andreas Stuhlmüller

Noah D. Goodman

wingated@mit.edu

ast@mit.edu

ngoodman@stanford.edu

Brain and Cognitive Sciences
MIT

Brain and Cognitive Sciences
MIT

Department of Psychology
Stanford

Abstract
We describe a general method of transforming
arbitrary programming languages into probabilistic programming languages with straightforward MCMC inference engines. Random
choices in the program are “named” with information about their position in an execution
trace; these names are used in conjunction witha database holding values of random variables
to implement MCMC inference in the space ofexecution traces. We encode naming information using lightweight source-to-source compilers. Our method enables us to reuse existing
infrastructure (compilers, profilers, etc.) withminimal additional code, implying fast models
with low development overhead. We illustrate
the technique on two languages, one functionaland one imperative: Bher, a compiled versionof the Church language which eliminates interpretive overhead of the original MIT-Church
implementation, and Stochastic Matlab, a new
open-source language.

1 INTRODUCTION
Probabilistic programming languages simplify the development of probabilistic models by allowing programmers to specify a stochastic process using syntax that
resembles modern programming languages. These languages allow programmers to freely mix deterministic
and stochastic elements, resulting in tremendous modeling flexibility. The resulting programs define prior distributions: running the (unconditional) program forward many times results in a distribution over execuAppearing in Proceedings of the 14t h International Conference on Artificial Intelligence and Statistics (AISTATS) 2011, Fort
Lauderdale, FL, USA. Volume 15 of JMLR: W&CP 15. Copyright2011 by the authors. Revision 3. Updated February 8, 2014.

tion traces, with each trace generating a sample of datafrom the prior. The goal of inference in such programs
is to reason about the posterior distribution over execution traces conditioned on a particular program output.
Examples include BLOG [7], PRISM [15], Bayesian Logic
Programs [5] Stochastic Logic Programs [9], Markov Logic[14], Independent Choice Logic [12], IBAL [11], Probabilistic Scheme [13], Λ◦ [10], and HANSEI [6].
We present a general technique for turning any programming language into a probabilistic programming language with an accompanying universal Markov chain
Monte Carlo inference engine. Our method allows thefull use of all language constructs permitted by the original (non-probabilistic) language by leveraging native
compilers (or interpreters). The key technical idea is to
give each random choice of a fixed program a unique“name” that depends on its position in a given execution
trace. We then convert stochastic functions into deterministic ones which use these names to look up their return values in a database that stores the values of all random choices. By controlling the values in this database,execution traces of the program can be controlled, whichallows us to construct the key operations needed for the
Metropolis-Hastings algorithm.
Different techniques for naming random choices resultin different MCMC dynamics. Our central contribution
is to name random choices based on their structural position in the execution trace, which allows fine-grained
sharing of random choice values between subsequent
MCMC states. We describe how naming information canbe provided by a side computation that is constructed bya source-to-source transformation at compile time: theoriginal probabilistic program can be augmented withnaming information to create a new program that makes
the names available. The resulting new program can be
executed at full speed, with only minimal overhead to
track names and control the database, and no changes tothe compiler or interpreter for the underlying language.
We illustrate our technique on two languages, one functional and one imperative. The first is named “Bher,”Lightweight Implementations of Probabilistic Programming Languages Via Transformational Compilationand is a compiler for the Church probabilistic programming language [4]. The current MIT-Church implementation uses a custom interpreter to reason about the distribution of execution traces, but since our method does
not incur interpretive overhead, the resulting programs
can run orders of magnitude faster. We also illustrate
the technique on an imperative language by creating
Stochastic Matlab, an entirely new language.

2 OVERVIEW OF THE METHOD
We begin by outlining our setup. We define an unconditioned probabilistic program to be a parameterless function f with an arbitrary mix of stochastic and deterministic elements (hereafter, we will use the term function
and program interchangeably). The function f may be
written in any language, but our running example will be
Matlab. We allow the function to be arbitrarily complex
inside, using any additional functions, recursion, language constructs or external libraries it wishes. The only
constraint is that the function must be self-contained,with no external side-effects which would impact the execution of the function from one run to another.
The stochastic elements of f must come from a set of
known, fixed elementary random primitives, or ERPs.
Complex distributions are constructed compositionally,using ERPs as building blocks. In Matlab, ERPs may befunctions such as rand (sample uniformly from [0,1]) or
randn (sample from a standard normal). Higher-order
random primitives, such as nonparametric distributions,
may also be defined, but must be fixed ahead of time.
Formally, let T be the set of ERP types. We assume thateach type t ∈ T is a parametric family of distributions
p t (x|θt ), where θt are the parameters of the distribution.
As f is executed, it encounters a series of ERPs. Algorithm 1 shows an example of a simple f written in Matlab.
Here, there are three syntactic ERPs: rand, randn, andgammarnd. During execution, depending on the return
value of each call to rand, different paths will be taken
through the program, and different ERPs will be encountered. We call this path an execution trace. A total of 2000
random choices will be made when executing this f .
Let f k|x1 ,···,xk−1 be the k’th ERP encountered while executing f , and let x k be the value it returns. Note that theparameters passed to the k’th ERP may change depending on previous x k ’s (indeed, its type may also change, as
well as the total number of ERPs). We denote by x all of
the random choices which are made by f , meaning that fdefines the probability distribution p(x). In our example,x ∈ R2000 . The probability of x is therefore the product ofthe probability of all of the ERP choices made:p(x) =

KY
k=1

p tk (x k |θtk , x 1 , · · · , x k−1 )Algorithm 1 A simple Stochastic Matlab program representing a Gaussian-Gamma mixture model.
1: for i=1:1000
2:
if ( rand > 0.5 )3:X(i) = randn;
4:
else
5:
X(i) = gammarnd;
6:end;
7: end;again noting explicitly that types and parameters may depend arbitrarily on previous random choices. To simplify
notation, we will omit the conditioning on the values ofprevious ERPs, but again wish to emphasize that thesedependencies are critical and cannot be ignored. By f k , itshould therefore be understood that we mean f k|x1 ,···,xk−1 ,and by p tk (x k |θtk ) we mean p tk (x k |θtk , x 1 , · · · , x k−1 ).
Generative functions as described above are, of course,
easy to write. A much harder problem, and our goal inthis paper, is to reason about the posterior conditional
distribution p(x \c |x c ), where we define x c to be a subset of random choices which we condition on and x \c
to be the remaining random choices. For example, we
may condition our example program on the values of
the X(i)’s, and reason about the sequence of rand’s most
likely to generate the X(i)’s. Rejection sampling is one
possible inference method: we may run f forward until
x c is generated, somehow recording x \c . Of course, suchan approach is intractable for all but the simplest models.
We propose a different approach, based on the following intuition. Consider repeatedly running an unconditioned f . Every time an f k is encountered, the program
samples an x k and continues execution. Importantly, if
two runs of f sample exactly the same values, their traces
will be the same. This suggests the following procedure
for controlling the execution of f :1. Give each f k in every possible execution trace a“name.” This name does not have to be unique, andcan depend on previous x k ’s or other information
about program state.
2. Rewrite the source code of f to generate f 0 , replacing random functions f k with deterministic functions f k0 . When encountered in the execution trace,these functions f k0 deterministically use their nameto look up a current value x k in a database and return it (and if no value exists, they sample x k ∼p tk (|θtk )., store it in the database, and then returnit). Behind the scenes, these functions also accumulate the likelihood of each random choice.
Thus, the execution trace of f 0 can be controlled by manipulating the values in the database. This makes MCMC

David Wingate, Andreas Stuhlmüller, Noah D. Goodman

Algorithm 2 An MCMC trace sampler.

Algorithm 3 Function trace_update( D )1: Initialize: [ D, l l , , ] = trace_update( ; )2: Repeat forever:3:
Select a random f k via its name n
4:Look up its current value (t , x, l , θd b ) = D(n).
5:
Propose a new value x 0 ∼ K t (|x, θd b )6:Compute F = log K t (x 0 |x, θd b )7:Compute R = log K t (x|x 0 , θd b )8:Compute l 0 = log p t (x 0 |θd b )9:
Let D0 = D
10:Set D0 (n) = (t , x 0 , l 0 , θd b )11:[ D0 , l l 0 , l l fresh , l l stale ] = trace_update( D0 )12:
α = l l 0 − l l + R − F + log |D| − log |D0 | + l l stale − l l fresh
13:
if ( log(rand) < α )14:// accept
15:
D = D016:ll = ll0
17:
else
18:
// reject; discard D0

1: Initialize l l , l l fresh , and l l stale to 0
2: Execute f :3: For all random choices k:
4:
Run computation until choice k,inference possible by enabling proposals, scoring, andthe ability to accept or reject proposals. We make aproposal to an execution trace by picking an x k in the
database, modifying it according to a proposal distribution, and then re-executing f 0 and computing the likelihood of the new trace. As f 0 is executed, some random
choices may be reused, new randomness may be triggered (which is sampled from the prior), and randomness
that was previously used may no longer be needed. In our
example, we can make a proposal to, say, the 500th rand.
This will change the execution trace, switching between
the Gaussian branch and the gamma branch.
We now address the details. Section 2.1 discusses theoverall MCMC algorithm, Section 2.2 discusses how random variables are named, and Section 2.3 discusses conditional inference.
2.1 MCMC IN TRACE SPACE
We now describe our MCMC algorithm. Here, we assume
that some naming scheme has been defined, and that amechanism for computing those names has been implemented; we will discuss how to do this in the next section.
Let a database D be defined as a mapping N → T × X ×
L × θ, where N is the name of a random choice, T is its
ERP type, X is its value, θ are ERP parameters, and L is
this random value’s likelihood. We allow missing entries.
We can control the execution of f 0 by controlling the values in D. When f 0 encounters an f k0 , it computes its
name n ∈ N, its parameters θc (where ‘c’ is a mnemonicfor “current”), and its current type t c ∈ T . It then looksup (t , x, l , θd b ) = D(n). If a value is found and the types

determining n, t c , θc for k.
5:
Look up (t , x, l , θd b ) = D(n)6:
if a value is found in the database and t = t c
7:if parameters match (ie, θc == θd b )8:ll = ll +l
9:
else
10:
// rescore ERP with new parameters
11:
Compute l = log p tc (x|θc )12:Store D(n) = (t c , x, l , θc )13:
ll = ll +l
14:
else
15:
// sample new randomness
16:
Sample x ∼ p t (.|θc )17:Compute l = log p t (x|θc )18:Store D(n) = (t c , x, l , θc )19:
ll = ll +l
20:
l l fresh = l l fresh + l
21:
Set return value of f k to x
22: For all random choices k in D:
23:
if k has not been active in this update:24:
Look up (t , x, l , θd b ) = D(k)25:l l stale = l l stale + log p t (x|θd b )26:
Remove k from D
27: return [ D, l l , l l fresh , l l stale ]match, we set the return value of f k to be x. Otherwise, the value is sampled from the appropriate ERP x ∼p tc (|θc ), its likelihood is computed, and the corresponding entry in the database is updated. We accumulate thelog-probability of new randomness, l l fresh , and the logprobability of randomness that is no longer used, l l stale .
This is formalized in the trace_update procedure, defined in Algorithm 3.
We use trace_update to define our overall MCMC algorithm, shown in Algorithm 2. Given a current trace x
and score p(x), we proceed by reconsidering one random
choice x k . We equip each ERP type with a proposal kernel
K t (x 0 |x, θ), which we use to generate proposals to x k . After a proposal, we call trace_update to generate a new
trace x 0 , and compute its likelihood p(x 0 ) and the probabilities of new and stale randomness, p fresh and p stale .
This gives us an overall score, which is used as an MH accept ratio¾½
p(x 0 ) K t (x|x 0 , θ) |D| p stale
α = min 1,
p(x) K t (x 0 |x, θ) |D0 | p fresh
There are a few subtleties to the implementation. Even
if an x for an f k0 is found in D, we must compare the pa-

1: m = poissrnd();
2: for i=1:m
3:
X(i) = gammarnd();
4: end;5: for i=m+1:2*m
6:
X(i) = randn;
7: end;rameters θc and θd b to make sure they are equal; if they
are not, we may still be able to re-use x, but must rescore it under the new parameters θc . We also remove unused (“stale”) random choices from the database to ensure that subsequent new randomness is sampled fromthe prior, and allows us to avoid reversible jump corrections. Finally, note that many of the probabilities in α
cancel when proposals are made from the local prior.

NaivenamingAlgorithm 4 An illustration how the number and type ofrandom choices can change.

Trace 1

1

2

3

4

5

Trace 2

1

2

3

4

5

Structural
namingLightweight Implementations of Probabilistic Programming Languages Via Transformational Compilation

Trace 1

1-1

3-1

3-2

3-3

Trace 2

1-1

3-1

3-2

(m=3)

(m=2)

(m=3)

(m=2)

poisson

gamma

6

7

6-1

6-2

6-3

6-1

6-2
gaussianFigure 1: An illustration of two naming schemes. Variables are boxes with names inside. Top: the problem
with naive variable naming. Notice the type mismatchat time 4 as well as the general sequence misalignment.
Bottom: structural naming yields the desired behavior,reusing the maximum number of random choices.

2.2 NAMING RANDOM VARIABLES
What properties do we desire in a naming scheme? We
build around a simple idea: when we propose a new x k0
and update the trace, we wish to reuse as many x k ’s aspossible. This will give us a new execution trace “close” tothe original, which is important to ensure a high acceptance rate in our MCMC sampler (recall that if an f k0 cannot reuse a value from D, one is sampled from the prior,which is unlikely to be good).
How might different naming schemes affect reuse of randomness? Algorithm 4 shows a simple program f which
we will use to illustrate different approaches. As f executes, it first encounters a Poisson, and then a sequenceof gamma and Gaussian variables. One might be tempted
to name random variables sequentially, according to theorder in which they are encountered during execution—that is, use k itself as the name. Figure 1 (top) illustrates
why this is a bad idea. Two traces are shown schematically, illustrating the f k0 ’s encountered. Suppose that intrace #1, the value of m is 3. A proposal is made, changing
m to 2, yielding trace #2. If variables are named sequentially, then x 1 and x 2 can be reused without problem, butwhen f 30 looks itself up in D, it discovers that its type haschanged, and must sample x 3 from the prior. The rest of
the sequence is now “misaligned,” and must probably bediscarded. More sophisticated versions of this problem
complicate this naming scheme to the point of unworkability, especially with complex models.
What is an alternative? Examining Alg. 4, we see thatwould like to name variables according to the way they
are used at a more semantic level. For example, we mightwant to name a variable “the third gamma variable in theloop on lines 2-4.”

Our key technical contribution formalizes this intuition.
We name random choices according to their structural
position in the execution trace, which we define in a way
roughly analogous to a stack address: a variable’s nameis defined as list of the functions, their line numbers, andtheir loop indices that precede it in the call stack. Thus, in
Alg. 4, variables are named “1” (the initial Poisson), “3-1,”,“3-2,” ... (the gamma variables sampled on line 3), and“6-1,”, “6-2,” ... (the Gaussian variables sampled on line
6). Figure 1(bottom) illustrates this on Alg. 4: we see thatthe maximum number of random variables are reused.
We now discuss naming schemes for imperative andfunctional languages.
Imperative.
For an imperative language, such as
Matlab, we define the structural position with an abstract stack trace that is augmented with a combination of function identifiers, line numbers, and loop iteration numbers. The formal specification is inductive,
as shown in Figure 2. In an interpreted imperative language, such as Python, Matlab or Ruby, this naming information could potentially be acquired via dynamic stackexamination, assuming suitable libraries exist. Alternatively, this information can be provided by augmenting
the source code for f , as we shortly discuss.
Functional. For a functional language, such as Lisp, we
do not need to worry about loops or line numbers: everything is a function, so we may define the structural position with a stack trace that is augmented with function
identifiers. The specification is shown in Figure 3.
There is no guarantee that this naming scheme is optimal
in any sense, but it does capture our desiderata, is simple
and fast, and incurs minimal overhead.

David Wingate, Andreas Stuhlmüller, Noah D. Goodman´´

A t op E  = ((lambda (addr) A E ) (top))A (lambda ( I in=1 ) E bod y ) = (lambda (addr . I in=1 ) A E bod y )A (mem E ) = ((lambda (maddr f) (lambda (addr . args) (apply f (cons args maddr) args))) addr A E )A (begin E in=1 ) = (begin A E i ni=1 )A (letrec (( I i E i )ni=1 ) E bod y ) = (letrec (( I i A E i )ni=1 ) A E bod y )A (if E t E c E a ) = (if A E t  A E c  A E a )A (define I E ) = (define I A E )A (quote E ) = (quote E )A (E op E in=1 ) = (A E op  (cons S addr) A E i ni=1 ) with S being a globally unique symbol.
A E  = E , otherwise.

Figure 3: The naming specification for Bher is given by a syntactic transform A t op .

• When entering a new function,◦ push a unique function id on the function stack.
◦ push a 0 on the line stack.
• When moving to a new line, increment the last valueon the line stack.
• When starting a loop, push a 0 on the loop stack.
• When iterating through a loop, increment the lastvalue on the loop stack.
• When exiting a loop, pop the loop stack.
• When exiting a function, pop the function stack andthe line stack.
Figure 2: The imperative naming specification. The
name of an f k is the state of the three stacks when f k is
encountered.

L (lambda ( I in=1 ) E bod y ) = (lambda ( I in=1 ) L E bod y )L (letrec (( I i E i )ni=1 ) E ) = (letrec (( I i D E i )ni=1 ) L E )L (begin E in=1 ) = (begin (force L E i )n−1
L E n )
i =1
L (mem E ) = (mem L E )L (quote E ) = (quote E )L (if E t E c E a ) = (if (force L E t ) D E c  D E a )L (E op E in=1 ) = ((force L E op ) D E i ni=1 )L E  = E , otherwise.
D (lambda I E ) = L (lambda I E )D (mem (lambda I E )) = L (mem (lambda I E ))
D E  = (list delayed (mem (lambda () L E ))), otherwise.
´• Begin executing f with empty function, line, and loop
stacks.

(define (force address val)(if (and (pair? val) (eq? (car val) ’delayed))(force address ((cadr val) address)))val)

Figure 4: Laziness transform L and definition of the function force referenced in the transform.

3 IMPLEMENTATION AND EXPERIMENTS
2.3 CONDITIONAL INFERENCE
In principle, a program can be conditioned on any expression that evaluates to a Boolean value, not just on a
fixed value assignment for a subset of the random variables of the program. The meaning of such a conditioned
program is defined by rejection sampling: Sample an execution trace from the unconditioned program. If thevalue of the conditioning expression is true, return thistrace. Otherwise, retry.
In practice, a MCMC implementation that allows general conditions requires that for any valid condition,
we can find an initial assignment of random variables
to values that satisfies this condition. If conditioning
is not restricted to random variables (Stochastic Matlab), then this is a search problem that needs to besolved. Possible approaches include constraint propagation (MIT-Church), annealing (Bher), and rejection sampling (Bher). We leave the analysis of this problem to future research.

We now present two implementations of our technique.
Bher is a compiled implementation of Church, andStochastic Matlab is a new, open-source language. Bothillustrate the virtues of our approach: small code bases,
high speed, and language flexibility.
3.1 EXAMPLE: BHER
Bher is an implementation of the Church probabilistic
programming language that compiles Church code to native Scheme code. Like MIT-Church, the original implementation of the Church language, Bher is written in
R6RS Scheme. However, MIT-Church is an interpreter—
not based on the lightweight implementation principle—
and thus provides a useful point of comparison.
In the following, we first discuss the addressing transform that names random choices in Bher, then compare
performance between Bher and MIT-Church on Hidden
Markov Models and Latent Dirichlet Allocation.

Lightweight Implementations of Probabilistic Programming Languages Via Transformational Compilation(begin(define geometric(lambda (p)(if (flip p)1(+ 1 (geometric p)))))(geometric .7))

3.1.1 Syntactic Transformations
The syntactic transform that provides addresses is shownin Figure 3. The address that is built at runtime consistsof a list of symbols. The transform is designed to change
the program code such that a new symbol is added to thefront of this list whenever the stack is extended.
The transform has two central parts: First, change eachfunction definition such that the function takes anotherargument, the address. Second, change each function
application such that (1) the current address is extended
with a symbol that uniquely identifies this (syntactic) application within the program and (2) pass on this modified address as an argument when the function is applied.

´

On the top level of the program, the address variable
addr is initialized to (top). Except for memoization(mem) and quote, all remaining syntax cases are handled
by applying the addressing transform to all proper subexpressions. Primitive functions like + and cons need tobe redefined to take (and ignore) the address argument,
e.g. to (lambda (addr . args) (apply prim args)).
In Bher, all globally free variables are treated as primitives
and redefined in this way.
In deterministic programs, memoization is an optimization technique that does not affect semantics: When amemoized function returns, it stores the return value forits current set of arguments in a cache. When the function is called a second time with the same arguments, thevalue is not computed again, but instead the value foundin the cache is returned.
In stochastic programs, we can distinguish mem as a semantic construct from mem as an optimization technique[4]. The distribution defined by (eq? (flip) (flip))
is very different from the distribution defined by(let ([mf (mem flip)]) (eq? (mf) (mf))). The former is true with probability .5, false with probability .5;
the latter is true with probability 1.
When addresses are available, semantic memoization isa simple syntactic transform. In a program that has beentransformed to provide addresses, random variables get
their values by looking up the value that is stored in the
global database at their address. Thus, to enforce that´

Algorithm 6 Transformed version of the Bher program
shown in Algorithm 5. a1 to a4 are the names generated
for function applications.
´

Algorithm 5 A Church program that samples from a geometric distribution.

((lambda (addr)(begin(define geometric(lambda (addr p)(if (flip (cons ’a1 addr) p)
1(+ (cons ’a2 addr)
1(geometric (cons ’a3 addr) p)))))(geometric (cons ’a4 addr) 0.7)))’(top))

two random choices always have the same value, we
make sure that the address they receive is the same. To
achieve this, the mem case of the addressing transformbuilds a function that captures within its closure the address of the mem creation and uses this address (extended
by any function arguments) instead of the address provided to the function when it is applied.
The Church program shown in Algorithm 5 and its transformed version (Algorithm 6) illustrate the address computation that is induced by the syntactic transform. The
program defines a simple geometric distribution with apotentially unbounded number of random choices, eachof which gets a unique address. When we propose tochange the value of the flip that controls the recursion depth, we want to reuse as much of the existing
trace as possible. An example of such a trace is shownin Figure 5. Note that apply nodes (shown in green) are
the only place where addresses change and that nodes
that refer to the same syntactic place (e.g. the two blue
(flip p) nodes) have different addresses when they occur in different places in the trace. A change to the return value of such a flip does not affect the addresses of
the choices “above” it. Thus, any randomness stored forthese choices persists across the change.
A major virtue of source-to-source transformations is
simplicity and compositionality. To illustrate this point,
we now sketch a complementary transformation, eager
style to lazy style, that can be used to avoid unnecessary
computation. As shown in Figure 4, we delay expressions
that must not necessarily get evaluated—for example, arguments to compound functions—and force expressions
that are always evaluated—operators, for example. To
delay an expression E, we wrap it in (mem (lambda ()E)). Thus, we ensure that it always returns the same result when evaluated multiple times. To force a value, we
check whether it is delayed and if it is, we simply apply
the memoized function stored in the value and recurse.
This difference in simplicity manifests itself in the codebase of the two projects. While the MIT-Church imple-

David Wingate, Andreas Stuhlmüller, Noah D. Goodman(geometric 0.3)
→ 2 (p=-1.561)'(top)apply→ 2 (p=-1.561)'(a5 top)(geometric p)→ 1 (p=-1.204)'(a5 top)

geometric→ procedure (p=0.0)'(top)apply→ 1 (p=-1.204)'(a3 a5 top)0.3
→ 0.3 (p=0.0)'(top)(if (flip p) 1 (+ 1 (geometric p)))→ 2 (p=-1.561)'(a5 top)(flip p)→ #f (p=-0.357)'(a5 top)apply→ #f (p=-0.357)'(a1 a5 top)flip
→ erp (p=0.0)'(a5 top)p→ 0.3 (p=0.0)'(a5 top)+→ primitive (p=0.0)'(a5 top)p→ 0.3 (p=0.0)'(a5 top)(if (flip p) 1 (+ 1 (geometric p)))→ 1 (p=-1.204)'(a3 a5 top)(+ 1 (geometric p))→ 2 (p=-1.204)'(a5 top)apply→ 2 (p=0.0)'(a2 a5 top)

geometric→ procedure (p=0.0)'(a5 top)1→ 1 (p=0.0)'(a5 top)

apply→ #t (p=-1.204)'(a1 a3 a5 top)

(flip p)→ #t (p=-1.204)'(a3 a5 top)1
→ 1(p=0.0)'(a3 a5 top)flip
→ erp (p=0.0)'(a3 a5 top)p
→ 0.3 (p=0.0)'(a3 a5 top)

Figure 5: Example of an execution trace for the geometric sampler defined by Algorithm 5. Each node is annotated with
its expression, value, log probability, and address. Applications (blue nodes) have subtraces for operator, operands, andfor the application of operator to operands (apply, green nodes; dashed border indicates a random choice).
Hidden Markov ModelOn the HMM model, Bher outperforms the MIT-Churchinterpreter. On the LDA model, Bher shows better performance at smaller problem sizes whereas MIT-Church
scales better to larger instances. This difference in scaling is due to (1) reuse of cached computation in mem and(2) computational short-circuiting in MIT-Church. In theprocess of updating a trace after a proposal has been
made, MIT-Church can detect that the information flowing into a function application has not changed; in thiscase, it does not recompute the return value. Since mem
is transformed away by the compiler, perhaps there isa source transformation that results in similar reuse ofcomputation. This is an area of future research.
3.2 EXAMPLE: STOCHASTIC MATLAB
We now illustrate our technique to create a new,imperative probabilistic programming language called“Stochastic Matlab,” which can take advantage of Matlab’s rich native libraries and toolboxes, and its parallelization, profiling and debugging support.
To rewrite f , we used the imperative ERP naming specification as outlined in Fig. 2. We implemented the rewriterusing a parser developed with code from the open-source120
20 40 60 80

Bher
MIT−Church

0Inference time (s)Bher
MIT−Church

20 40 60 80

We have tested inference performance in MIT-Church
and Bher on Hidden Markov Models and Latent Dirichlet
Allocation (Figure 6). In both cases, we conditioned on an
observation sampled from the model. For each HMM, we
generated a sentence of length 15. For LDA, we used 10
topics, 10 words per document, and 30 possible words,
varying the number of observed documents. All experiments use the optimizing compiler Ikarus [3].

Latent Dirichlet Allocation

0

3.1.2 Experiments

120

mentation has about 11, 119 lines of Scheme code, Bher
consists of only 1, 658 lines of code.

0

20

40

60

80

Number of latent states

100

0

10

20

30

40

Number of documents

Figure 6: Inference performance for the HMM (10,000
samples) and LDA model (1,000 samples).

Octave project [2]. ERP names were rewritten to deterministic functions using a MEX-based database. Remaining code was implemented in pure Matlab.
We now present a case study that relies on a deterministic GPU-based MEX function. The point is twofold: first,
to show how easy it is to mix arbitrary deterministic elements with stochastic ones—reusing complex external
libraries—and second, to show how easy it can be to create complex models with a few simple commands.
Our example is a small (but real) application of inverse
geological modeling. Layers of rock and sediment are deposited over millions of years to create a volumetric cubeof rock with porosity that varies according to the geometries of the layers. Given this forward model, we wish tocondition on a given volume, and recover a descriptionof the individual layers that built up the volume.
Figure 7 shows the model, which generates a rock volume by first sampling parametric descriptions for eachlayer (the details are not shown; layers are described asa weighted set of basis functions), and then “rendering”
them into a porosity volume (with Gaussian noise added
at each voxel). The render_layers function encapsu-

Lightweight Implementations of Probabilistic Programming Languages Via Transformational Compilation

A layer model of rock porosity
function dp_render_rock( data )num_layers = poissrnd();dp_sample_layer = dpmem( @sample_layer, 1.0 );for i=1:num_layers
layer(i,:) = dp_sample_layer();end;
rock_volume = render_layers( layers );data = rock_volume + randn( size(rock_volume) );return;6

Inferred porosity0
Log−likelihood

True porosity

x 10

−2
−4

Flat
DP

−6
0

200
400
Steps of inference

Figure 7: Results on modeling rock porosity.

lates the MEX/GPU-based renderer. Given this model, we
perform inference using tempered MCMC, conditioning
on the data shown in Fig. 7. A sample from the posteriorshows a reconstructed volume (also shown in Fig. 7) thatclosely matches the true volume—the program has successfully inferred a parametric form for each layer.
The program calls a function named dpmem, a stochastic memoizer similar to that found in Church. dpmem
accepts a function representing a base measure and aconcentration parameter, and returns a new function
which samples according to a Dirichlet process (DP). We
use this to create a DP in “layer space,” where parameters of layers can be shared. Fig. 7 (lower left) shows
a comparison of inference with and without the dpmem:
shared layers enable the inference engine to find higherlikelihood models faster. Hierarchical DPs are equally
easy: we could (for example) let f=dpmem(@randn,1.0),
and g=dpmem(@f,1.0) to create an HDP with a Gaussianbase measure. We consider such flexibility in defining
models—combined with the rich libraries of Matlab—tobe the major virtue of our approach.

4 RELATED WORK
Our MCMC algorithm is a new implementation of the
Church MCMC algorithm [4], and is similar to the algorithm used in BLOG [7]. However, we provide a simpler
and more uniform approach to constructing such algorithms for new languages, with small, clean code bases.
The closest related work is that of Kiselyov and Shan [6],who share the goal of transforming standard languages

into probabilistic versions with little interpretive overhead. They use delimited continuations to control random choices without explicit choice naming. However,their technique is only applicable (prima facie) to enumeration or importance sampling, not MCMC.
Other related approaches include BUGS [16], HBC [1],and Infer.NET [8]. BUGS is not transformational (it simply loads a data structure into memory from a config file),
but the implementations of both Infer.NET and HBC arecomparable to our approach. Infer.NET uses a complicated transformation to compile csoft to C#. HBC compiles from a config file to C, but the specification language used by HBC is not a full-fledged language. We
consider our approach simpler and more general.
With respect to modeling flexibility, our approach provides a simple dynamic knowledge-based model construction, which allows unbounded models, including
complex recursions, to be handled efficiently. Models
such as the naive implementation of the geometric distribution in Bher (Alg. 5), the gamma/randn example in
Stochastic Matlab (Alg. 4), and nonparametric distributions such as the Dirichlet Process (Alg. 7) cannot be
specified in these other languages.
A performance comparison between these languages isimportant, but beyond the scope of this paper. A careful comparison must account for issues such as mixing
rates, initialization values, differences in parameters ofproposal kernels, etc. and so we have explicitly limited
this paper to an explanation of the technique, with comparisons reserved for a longer future paper.

5 CONCLUSION
We have described a lightweight implementation technique for probabilistic programming languages. The
method is simple and fast, and the resulting languages
permit a great deal of flexibility in the specification of
distributions by mixing stochastic and deterministic elements with arbitrary language features (such as objects,
inheritance, operator overloading, closures, recursion, libraries, etc.). Our example implementations have compact code bases and reasonable inference speed.
The main directions for improvement are better mixing
and faster inference. Because source-to-source transformations are often compositional, more transformations
could be applied to enhance performance. These could
reduce redundant computation between traces, generatecompound proposals, or implement constraint propagation for initializing conditioners. Future work will investigate these issues, as well as the possibility of new languages, such as stochastic Python.

David Wingate, Andreas Stuhlmüller, Noah D. Goodman

ACKNOWLEDGMENTS
The authors gratefully acknowledge the support of ONR,
the James S. McDonnell Foundation and AFOSR grantFA9550-07-1-0075.

References[1] H. Daume III. Hbc: Hierarchical bayes compiler,
2007.
[2] J. W. Eaton. GNU Octave Manual. Network Theory
Limited, 2002.
[3] A. Ghuloum. Ikarus scheme user’s guide, 2008.
[4] N. Goodman, V. Mansinghka, D. Roy, K. Bonawitz,
and J. Tenenbaum. Church: a language for generative models. In Uncertainty in Artificial Intelligence(UAI), 2008.
[5] K. Kersting and L. D. Raedt. Bayesian logic programming: Theory and tool. In L. Getoor and
B. Taskar, editors, An Introduction to Statistical Relational Learning. MIT Press, 2007.
[6] O. Kiselyov and C. Shan. Embedded probabilistic programming. In Domain-Specific Languages,
pages 360–384, 2009.
[7] B. Milch, B. Marthi, S. Russell, D. Sontag, D. L. Ong,
and A. Kolobov. Blog: Probabilistic models with unknown objects. In International Joint Conference on
Artificial Intelligence (IJCAI), pages 1352–1359, 2005.
[8] T. Minka, J. Winn, J. Guiver, and D. Knowles. Infer.NET 2.4, 2010. Microsoft Research Cambridge.
http://research.microsoft.com/infernet.
[9] S. Muggleton. Stochastic logic programs. In New
Generation Computing. Academic Press, 1996.
[10] S. Park, F. Pfenning, and S. Thrun. A probabilistic
language based on sampling functions. ACM Trans.
Program. Lang. Syst., 31(1):1–46, 2008.
[11] A. Pfeffer. IBAL: A probabilistic rational programming language. In Proc. of the 17th Int. Joint Conf. on
Artificial Intelligence, pages 733–740. Morgan Kaufmann Publ., 2001.
[12] D. Poole. The independent choice logic and beyond.
pages 222–243, 2008.
[13] A. Radul. Report on the probabilistic language
scheme. Technical Report MIT-CSAIL-TR-2007-059,
Massachusetts Institute of Technology, 2007.
[14] M. Richardson and P. Domingos. Markov logic networks. Machine Learning, 62:107–136, 2006.
[15] T. Sato and Y. Kameya. Prism: A symbolic-statisticalmodeling language. In International Joint Conference on Artificial Intelligence (IJCAI), 1997.

[16] D. Spiegelhalter, A. Thomas, N. Best, and W. Gilks.
Bugs - bayesian inference using gibbs sampling version 0.50, 1995.