Learning Stochastic Inverses for
Adaptive Inference in Probabilistic Programs

Andreas Stuhlmüller
Department of Brain and Cognitive Sciences
Massachusetts Institute of Technology

Noah D. Goodman
Department of Psychology
Stanford UniversityAbstract
We describe an algorithm for adaptive inference in probabilistic programs. During sampling, the algorithm accumulates information about the local probability
distributions that compose the program’s overall distribution. We use this information to construct targeted samples: given a value for an intermediate expression,
we stochastically invert each of the steps giving rise to this value, sampling backwards in order to assign values to random choices such that we get a likely parseof the intermediate value. We propose this algorithm for importance sampling and
as a means of constructing blocked proposals for a Metropolis-Hastings sampler.

1Introduction

Probabilistic programming is a recent merger between programming languages and Bayesian statistics that enables rapid development of complexly structured probabilistic models. Given a probabilistic program, the problem of inference is typically formulated in terms of conditioning: theprogram samples from a joint distribution p(x, y), and our goal is to sample from a conditional
distribution p(x|y). Implementations of probabilistic programming languages provide generic algorithms that attempt to solve this problem.
Two considerations guide us in the design of new inference algorithms: first, we aim to provide
algorithms that efficiently solve common inference problems, and, second, we aim to model humanreasoning. Taken together, these goals suggest two desiderata for algorithms: We want inference tobe adaptive in the sense that our algorithms should improve at solving a given inference problem
over time, inspired by the way people learn to reason about a problem type. We further want inference to be compositional in the sense that the algorithms decompose a big, challenging query into
smaller queries that are easier to solve.
In the following, we describe how to extend existing importance sampling and MCMC algorithms
such that they more strongly exhibit these properties. We outline the dimensions along which adaptive inference algorithms vary and explain where our proposed algorithm is located in this space. In
our analysis, we focus on functional probabilistic languages, using Church (Goodman et al., 2008)
and its foundation, the stochastic lambda calculus, as representatives for this class.

2

Learning stochastic inverses

The goal of adaptive inference is to improve the quality and/or efficiency of sampling over time,ideally converging to an efficient source of exact samples from the posterior distribution. Adaptivealgorithms vary along a number of dimensions:
What are the units we learn? There are many choices, including parameters that depend on what
underlying inference algorithm we use, e.g., the step size and number of leapfrog steps for Hamiltonian Monte Carlo. We will focus on learning about the distribution of the model.
1(define (noisify value)(gaussian value 2))(define reflectance(gaussian 1 1))(define illumination(+ (gamma 9 .5)(gaussian 1 .5)))(define luminance(* reflectanceillumination))(noisify luminance)(gaussian value 2)

apply

*

noisify(gaussian 1 1)

apply

(noisify luminance)(* reflectance illumination)+(+ (gamma 0 .5) (gaussian 1 .5))

(gamma 0 .5)(gaussian 1 .5)

Figure 1: Left: A probabilistic program modeling brightness constancy in visual perception. Right:
an execution trace corresponding to the program’s evaluation.
In functional probabilistic programs, this distribution exhibits a clear independence structure. The
evaluation of a functional program induces a tree of nested function applications—an execution
trace. Figure 1 shows an example. For any application (A1 A2 . . . An ), the probability of returninga value v isZZ
p(v) = · · · p(a1 )p(a2 ) . . . p(an )p(a1 (a2 , . . . , an ) = v) da1 . . . dan .
This suggests that candidates for learning include priors p(a1 ), p(a2 ), . . . , p(an ), conditional distributions p(a1 , a2 , . . . , an |v), and joint distributions p(a1 , a2 , . . . , an , v). The dynamic programming algorithm presented in Stuhlmüller and Goodman (2012) can be viewed as an approach tolearning about priors. Adaptive importance sampling algorithms such as AIS-BN learn about conditional distributions p(r|c) where r is a primitive random choice and c the global query condition(Cheng and Druzdzel, 2000). We are interested in learning about local conditional distributionsp(a1 , a2 , . . . , an |v), which we call stochastic inverses. Taken together, these inverses provide analternative factorization of the program’s joint distribution. For example, for the single application shown above, p(a1 )p(a2 ) · · · p(an )p(v|a1 , . . . , an ) is replaced with the inverse factorizationp(v)p(a1 , a2 , . . . , an |v).
How do we learn the units given program executions? In other words, how do we turn theinformation we gather during program evaluation into estimates for the units of interest? In general,many machine learning methods are of potential relevance, including discriminative methods.
We take a local density approximation approach.1 Whenever we encounter an application(A1 A2 . . . An ) during evaluation, we store the values ~a := a1 , a2 , . . . , an together with application return value v. Figure 2 shows the local joint distributions for the program in Figure 1. We
can view these joint distribution samples as an implicit representation of the approximate conditional distribution: In order to sample from the conditional distribution p(~a|v) given a new value v ∗ ,we find the subset of points {(~a(i) , v (i) ), . . . , (~a(j) , v (j) )} containing nearest neighbors to v ∗ . We
construct a kernel density estimate of the preimage distribution from these neighbors, and sample apreimage ~a∗ from this density (Algorithm 1).
For these estimates to be accurate, our sampled program executions need to provide unbiased estimates of the distributions of interest. This can be achieved by computing the estimates using samples
from the unconditioned program. In the discussion, we consider the case where the program executions are generated in the process of approximate inference.
How do we use the units to change the inference process? This depends on the underlying
inference algorithm. We consider two cases, importance sampling and MCMC.
In order to use stochastic inverses for importance sampling, we simply sample according to thealternative factorization described above. Unconditioned, this process starts by sampling from the1Density estimation techniques such as Gaussian Processes and infinite Gaussian mixtures are alternatives.
However, quality of estimation needs to be balanced against computational cost. It is also unclear whether thepriors expressed by such nonparametric techniques mesh with the distributions we actually encounter.

2

−100

10203046reflectanceillumination

0 2 4 6 8 10 1201080 2 4 6 8 10 120246108(gaussian 1 0.5)26illumination

luminance−20−10 0 10 20 3030
20100(noisify luminance)−10−4 −2 024(gamma 9 .5)

luminance

Figure 2: Local joint distributions for the brightness constancy program shown in Figure 1.
Algorithm 1: Given value v for expression e, sample immediate subexpressions.
function S TOCHASTIC I NVERT(e, v ∗ )find nearest neighbors for value v ∗ in expression data {(~a(1) , v (1) ), (~a(2) , v (2) ), . . . , (~a(m) , v (m) )}construct density estimate pK on nearest-neighbor subexpression values {~a(i) , . . . , ~a(j) }draw ~a∗ from density estimate
return ~a∗ , pK (~a∗ )marginal distribution on the return value of the program. Without loss of generality, we assume that
this value represents the Boolean value of the conditioning expression. We set this value to trueand sample backwards using our inverses. This process bottoms out at elementary random choices,which allows us to read off the prior probability p of the sampled execution trace and to computethe importance weight p/q of our sample (Algorithm 2).
For Metropolis-Hastings, we build on the Church MCMC algorithm that executes a random walkon program execution traces (Goodman et al., 2008). To make proposals, this algorithm selects arandom choice, proposes a new value, and propagates the changes through the execution trace. We
augment this algorithm by introducing a new type of proposal: select a random application node in
the execution trace, then regenerate the subtree below this node using stochastic inverses. This sets
all random choices in this tree, creating a blocked proposal. Once the random choices are set, we
propagate the changes analogous to the original Church MCMC algorithm.

3Examples: XOR, brightness constancy

Consider the program (xor (flip)(flip)) conditioned to return true. Single-site MCMC cannot switch between the two valid assignments [1 0] and [0 1]. In contrast, our algorithm can constructblocked proposals that swap both assignments simultaneously: given knowledge of the joint distribution, our algorithm will resample the arguments to xor in a way that leaves the condition satisfied.
This could be extremely useful if the xor expression occurred within a larger program.
As a second example, consider the brightness constancy program shown in Figure 1. Figure 2 shows
the corresponding local joint distributions. To sample from the program under the condition that wenoisily observed a luminance of 3, we proceed as follows (left to right in the Figure): conditioned on(noisify luminance) being 3, we sample a luminance l. Next, conditioned on luminance l, wesample a reflectance r and illumination i. Finally, conditioned on illuminance i, we sample values
for the gamma and Gaussian random primitives. Figure 3 shows that the empirical variance of theimportance weights for a sampler that follows this procedure goes to 0 as we refine the inverses
using forward samples (percentiles computed across 50 runs). If the stochastic inverses are exact,
this produces a posterior sample in a way that satisfies the goals we set in the introduction: we havedecomposed the overall problem into smaller problems of sampling from local conditional distributions (stochastic inverses), and adaptive learning of these inverses gradually improves inference.
31000.0
1.0 10.0

procedure I NVERSE I MPORTANCE(e, v)~a, q0 ← S TOCHASTIC I NVERT(e, v)
for i = 1 to |e| doif ei is a primitive random choice thenstore value ~ai for choice eipi , qi ← pei (~ai ), 1else
pi , qi ← I NVERSE I MPORTANCE(ei , ai )return p1 · · · pn , q0 q1 · · · qn

median
75th
90th

0.1

Variance of weights

Algorithm 2: Recursively importance-sample
expression e conditioned on return value v.

500

1000

1500

2000

2500

3000

Number of samples

Figure 3: On the right: Results for adaptive importance sampling applied to the program in Figure
1. As the number of samples used to estimate stochastic inverses increases, the empirical varianceof the importance weights approaches 0 (implying convergence to a perfect posterior sampler).

4

Related work and discussion

There exists prior work on both adaptive importance sampling and adaptive MCMC. A number ofadaptive importance sampling algorithms have been proposed for Bayesian networks, e.g., Shachterand Peot (1989), Cheng and Druzdzel (2000), and Ortiz and Kaelbling (2000). These techniques
typically learn importance distributions with the same independence structure as the prior, whichin general excludes exact representation of the posterior. Similarly, adaptive MCMC techniquessuch as those presented in Roberts and Rosenthal (2009) and Haario et al. (2006) are used to tune
parameters of the MCMC algorithms, but do not allow arbitrarily close adaptation to the posterior,
whereas our method is designed to allow such close approximation. Future work will be required toevaluate whether our methods for importance sampling and MCMC provide better results in practicethan these other approaches.
A key question for our method is whether we can use samples that are produced from a distributionthat is not the prior and still guarantee eventual convergence of the estimated inverse distributions.
Samples generated from the prior may be sparse in regions that have high probability under theposterior, resulting in slow convergence of our stochastic inverses for the problem of interest. We canencourage higher coverage by temporarily increasing the entropy of our random variables, or target
the region of interest by using the samples generated during MCMC, however these alternatives bothresult in biased inverses. When can we correct for this bias, or rely on bias washing out over time?
Acknowledgments
We thank Thomas Icard, Jacob Steinhardt and Lewis Fishgold for useful discussions, and ONR grant
N00014-09-0124 for funding.

ReferencesJ. Cheng and M. Druzdzel. AIS-BN: An adaptive importance sampling algorithm for evidential reasoning inlarge bayesian networks. Journal of Artificial Intelligence Research, 2000.
N. D. Goodman, V. Mansinghka, D. M. Roy, K. Bonawitz, and J. B. Tenenbaum. Church: a language forgenerative models. Proceedings of the 24th Conference in Uncertainty in Artificial Intelligence, pages 220–229, 2008.
H. Haario, M. Laine, A. Mira, and E. Saksman. DRAM: efficient adaptive MCMC. Statistics and Computing,16(4):339–354, 2006.
L. E. Ortiz and L. P. Kaelbling. Adaptive importance sampling for estimation in structured domains. In Proc. ofthe 16th Ann. Conf. on Uncertainty in A.I. (UAI-00), pages 446–454. Morgan Kaufmann Publishers, 2000.
G. Roberts and J. Rosenthal. Examples of adaptive MCMC. Journal of Computational and Graphical Statistics,
18(2):349–367, 2009.
R. D. Shachter and M. A. Peot. Simulation approaches to general probabilistic inference on belief networks. InProc. of the 5th Ann. Conf. on Uncertainty in A.I. (UAI-89), pages 311–318, New York, NY, 1989. ElsevierScience.
A. Stuhlmüller and N. D. Goodman. A dynamic programming algorithm for inference in recursive probabilisticprograms. Second Statistical Relational AI workshop at UAI 2012 (StaRAI-12), 2012.

4