Conceptual background: Models as programs

A Bayesian hierarchical model can be expressed as a program that samples random variables from their conditional prior distributions and observes data. Inference samples from the posterior distribution over the program's execution traces. This page explains how BAli-Phy performs that inference and updates model computations efficiently.

See also the model-program examples.

Model programs and execution traces

Running a model program samples its random variables and evaluates the likelihood of the observations. In a hierarchical model, the distribution of one random variable can depend on values sampled earlier. Together, these sampling distributions define the joint prior.

For example, the program below samples a mean and standard deviation, then observes two heights from normal distributions with those parameters.

Two executions might produce:

The sequence of random choices made during an execution is its execution trace. For a fixed program and fixed inputs, these choices determine the execution. Replaying a trace uses its recorded choices to reproduce the execution and its probability densities. Modifying choices in the trace can change the subsequent execution.

Random choices can determine which conditional branch is evaluated, how far recursion proceeds, and which additional random variables are created. The program thus defines a distribution over executions whose structures, as well as parameter values, may differ.

BAli-Phy implements a universal probabilistic programming language, supporting random control flow and changing collections of random variables. A fixed graphical model provides a useful special case: values change, but the organization of random variables and their dependencies remains fixed.

The term probabilistic programming language is also used for more restricted systems. The relevant distinction here is whether random choices can change the structure of the probabilistic computation. For a discussion in evolutionary biology, see Ronquist et al. (2021).

Conditioning on observations

The prior density of a trace is the product of its conditional sampling densities; the random variables need not have independent priors. Calls to observe contribute likelihood factors. Together these define the posterior over traces:

posterior density ∝ prior density × likelihood

One approach to posterior inference is to generate independent traces from the prior and weight them by their likelihoods. This is importance sampling with the prior as the proposal distribution. It can be inefficient when most prior draws have negligible likelihood compared with a small fraction of draws, which then carry most of the weight. BAli-Phy instead uses Markov chain Monte Carlo (MCMC).

MCMC: modifying and replaying traces

The host environment executes the model program and performs inference over its traces. The model program specifies distributions and observations; it need not implement the MCMC loop. In BAli-Phy, the host environment includes the interpreter and inference runtime and is implemented mainly in C++.

A Metropolis–Hastings update involves:

1. Proposing a change to one or more random choices in the current trace.
2. Replaying the modified trace and evaluating the resulting prior and likelihood.
3. Accepting or rejecting the proposal using the Metropolis–Hastings acceptance probability, accounting for both the target density and the proposal mechanism.

A proposal can change control flow as well as parameter values. Random choices on the new execution path may need to be generated, while choices on the old path may no longer be used. The proposed trace must correspond to a valid execution of the model program. For example, changing a tree topology can change which parent node a node's state depends on.

Model programs can also register specialized MCMC moves, as the custom codon example illustrates.

Replay describes these updates conceptually. BAli-Phy implements them by updating a retained execution dependency graph, avoiding full re-execution after each proposal.

Incremental computation and dependency tracking

BAli-Phy uses incremental computation to update a previous execution when designated modifiable variables, such as random choices, receive new values. Its execution dependency graph records how computations depend on these variables and on one another. These dependencies identify computations that need reevaluation and allow unaffected results to be reused.

Nodes in this graph represent computations and values; edges record dependencies of results on their inputs. It includes deterministic intermediate calculations, not just random variables and probability factors.

Model programs are written in Haskell. A pure computation can reuse its result when the values on which it depends remain unchanged. BAli-Phy's interpreter records these dependencies as it evaluates the program, including dependencies on control flow: the result of a conditional expression depends on which branch is selected as well as on the computation within that branch.

Changes in control flow can change the computations performed and hence the dependency graph itself. Unlike a fixed probabilistic graphical model, this graph describes the computation of a particular execution.

The virtual machine

The interpreter implements a virtual machine modeled on Sestoft's (1997) Deriving a lazy abstract machine, extended to support incremental computation. Heap entries hold closures: expressions together with environments that bind their free variables to other heap entries. Evaluation is call-by-need: expressions are evaluated when demanded, and their results are shared. The machine reduces expressions to weak head normal form, leaving the components of data structures unevaluated until needed.

Ordinary lazy evaluation can overwrite an expression with its result. Incremental evaluation must preserve enough information to recompute it when an input changes. Modifiable variables are changed directly; computations that depend on them are changeable. When a reduction step does not depend on modifiable values, the machine can replace the closure with the result of that step. Otherwise, it retains the original closure and separately records both the result of the reduction step and the eventual result in weak head normal form. This preserves the expression for reevaluation while allowing direct access to its cached result.

Dependency edges record which values a reduction uses and which computation it produces. Reverse edges let the machine propagate invalidation from a changed input to affected reductions and cached results. A cached final result may become invalid while the reduction step remains valid; the machine can then reevaluate the computation produced by that step without repeating the step itself. Since reductions can allocate closures and select different branches, this mechanism supports changes in the structure of the computation as well as changes in numeric values.

Continue with the model-program examples.