An Embedded Probabilistic Programming Language

A probabilistic program describes a generative model in ordinary code: it draws latent parameters from priors, then ties those parameters to observed data through likelihoods. A probabilistic programming language (PPL) separates the model from the inference. You write the model once; the language supplies a generic engine that returns the posterior distribution of the parameters given the data. The same model can be fit by several different engines without rewriting a line.

This chapter builds a small but complete PPL in Common Lisp. It provides a macro DSL for declaring models and three inference engines written from scratch: random-walk Metropolis-Hastings, Hamiltonian Monte Carlo with automatic differentiation, and mean-field variational inference. A terminal REPL and text plots let you inspect the results.

The example program for this chapter is in the file 11_probabilistic_dsl.lisp.

Chapter 9 introduced Bayesian updating for a single conjugate pair, where the posterior fell out as a closed-form Beta. Real models rarely have that luxury. A probabilistic program encodes the joint distribution over latent parameters and data . The program draws each latent from a prior with sample and ties latents to data through a likelihood with observe. Inference targets the posterior

where the marginal likelihood is usually an intractable integral over the parameter space.

The key observation that makes all three engines in this chapter work is that we never need . Every acceptance ratio and every gradient involves a ratio or a difference of log densities, and the constant cancels. So each engine needs only the unnormalized log posterior

The program computes for any . The engine’s job is to explore the shape of and report where its mass sits. Sampling engines return a set of draws whose empirical distribution approximates the posterior; the variational engine returns a simple approximating distribution tuned to match it.

A prior often lives on a bounded set. A standard deviation is positive; a probability lies in . The samplers in this chapter move on the whole real line, so each latent is reparameterized. We work in an unconstrained space and map back to the natural scale with a transform : for positive support, for the unit interval, and a scaled sigmoid for a bounded interval.

A change of variables bends the density. For the transformed variable the density picks up the Jacobian:

For the Jacobian term is . For the sigmoid it is . The DSL applies these corrections automatically, so the model author never writes a transform or a Jacobian by hand. This keeps the model readable and keeps the math correct.

There are three standard computational routes from to the posterior, and this chapter implements all three.

Metropolis-Hastings (MCMC). Start at a point . Propose a new point . Accept it with probability

If accepted, the chain moves to ; otherwise it stays. Under mild conditions the chain’s stationary distribution is the posterior. The method is simple and needs no gradients, but it mixes slowly in high dimensions because random proposals wander.

Hamiltonian Monte Carlo. Treat as a potential energy and add a momentum variable with kinetic energy . The total energy is

Simulate the Hamiltonian dynamics with the leapfrog integrator for several steps, which slides the state along the posterior’s level sets. Accept the endpoint with probability . Long, low-rejection moves need the gradient of , which we get exactly from automatic differentiation rather than finite differences. HMC mixes far better than random-walk Metropolis in moderate dimensions.

Variational inference. Instead of simulating the posterior, fit a simple distribution to it. We use a factorized Normal . Fitting means maximizing the Evidence Lower Bound

where is the entropy of . The reparameterization trick writes a sample as with , so the expectation becomes a differentiable function of and . The same automatic differentiation supplies the gradient, and Adam ascends the ELBO. Variational inference is fast and deterministic, but it approximates the posterior with a shape (independent Gaussians) that may not match the truth.

The next sections show how the file builds the pieces these engines share: automatic differentiation, a distribution library, support transforms, and the model DSL.

HMC and variational inference both need the gradient of . The file implements forward-mode automatic differentiation with dual numbers. A dual number carries a value together with the full gradient vector of that value with respect to the model parameters. Every primitive operation propagates the gradient by the chain rule, so a single evaluation of the log density yields both its value and its exact gradient.

For a unary function with derivative , the rule is . For a binary function the partials combine by the chain rule. Constants stay plain doubles, so Metropolis-Hastings, which needs no gradient, pays no AD overhead.

On top of these two combinators, the file defines differentiable versions of the arithmetic operators. Each one passes the right partials through d-binary or d-unary:

The entry point seeds coordinate with the unit dual (value everywhere except partial set to ), runs the function once, and reads the value and full gradient off the result. One pass gives the whole gradient vector:

The cost is one evaluation of per parameter dimension, which is fine for the small models in this chapter. (Reverse-mode AD would be cheaper for functions with many inputs and one output, but forward mode is simple to build and exact.) The numerical guards in and matter: a divergent HMC trajectory can push a parameter to an extreme value, and clamping keeps the gradient finite so the trajectory can be rejected rather than crashing the sampler.

Each distribution is a small record holding three things: a log-density function written with the differentiable operators, a plain sampler used to initialize chains, and the support of the distribution (:real, :positive, :unit, :bounded, or :discrete). The support fixes the transform to unconstrained space.

The log-density is built from g-, g*, g/, gsquare, and glog, so it carries gradients when its inputs are duals and computes plainly when they are doubles. The parameters and may themselves be duals, which is what makes hierarchical models work: a likelihood mean can depend on other latents, and the gradient still flows through.

The library supplies Normal, Half-Normal, Exponential, Gamma, Beta, Uniform, Bernoulli, and Poisson. The Beta and Gamma log densities use a Lanczos approximation to , and the Poisson likelihood uses . These are the same special functions that appeared in Chapter 5 (continuous distributions) and Chapter 9 (the Beta-Bernoulli update), restated here so the file stands alone.

Two small functions handle the mapping between constrained and unconstrained space. constrain-support maps an unconstrained real back to the distribution’s support, and returns the log Jacobian of that map. Both are written with the differentiable operators so the Jacobian term enters the gradient correctly.

For support the Jacobian term is just , because and . For support it is , the log of the sigmoid derivative. The sampler side uses a plain and to invert the transform when placing initial values.

The macro defmodel defines a function that, when called with data, returns a model object. The model body runs under a dynamic context in one of two modes:

• :init discovers the latents in order and picks dispersed starting values.
• :density reads a parameter vector and accumulates the unconstrained log posterior (prior plus Jacobian plus likelihood).

The body never mentions transforms or gradients. sample and observe handle all of it.

sample declares a latent. In mode it records the latent’s support and returns a starting value. In mode it reads the next parameter from the vector, transforms it to the natural scale, and adds the prior log density plus the Jacobian to the running total:

The same body serves both modes because sample and observe dispatch on ctx-mode. Running the body once in mode builds the latent list and the start vector; running it in mode with a parameter vector builds the log posterior. Two thin wrappers package these:

The dynamic variable carries the mode and the accumulator, so the model body reads like a direct description of the generative process. This is the whole point of an embedded DSL: the model is just a Lisp function, and the host language’s machinery does the rest.

The file ships three data generators, one per built-in model. Before fitting, it helps to see what the data look like.

The three formats are:

• Coin: a list of and , here ones and zeros: (1 1 1 0 1 ... 1 0 0).
• Normal: a vector of doubles drawn from , for example #(4.83 6.21 3.40 5.57 ...).
• Regression: two parallel vectors and , where each with .

Because the global random state is freshly seeded on every load, the exact data (and so the exact posterior) differ slightly each session. The true parameters stay fixed: coin bias near , normal mean and standard deviation , regression intercept , slope , noise .

The simplest engine is a random-walk Metropolis sampler in unconstrained space. At each iteration it adds isotropic Gaussian noise to the current point, evaluates at the proposal, and accepts with the Metropolis ratio. Divergent proposals (infinite or NaN log density) are rejected. The engine records the acceptance rate and returns posterior draws transformed back to the natural scale.

The proposal scale controls the tradeoff between acceptance rate and move size. Too small and nearly every proposal is accepted but the chain crawls; too large and nearly every proposal is rejected and the chain stalls. A rate around to usually mixes well for random-walk Metropolis.

HMC treats as a potential energy. Each iteration draws a fresh momentum , then simulates the Hamiltonian dynamics with the leapfrog integrator. The integrator alternates half-step momentum updates with full-step position updates, using the gradient of the potential at each step. The gradient comes from ad-gradient, so it is exact.

Two details matter in practice. The path length is jittered each iteration (steps = 1 + random(steps0)) so the integrator cannot resonate with the target’s curvature, which would cripple mixing. And a divergent trajectory (the potential becomes infinite) breaks out of the leapfrog loop early and is rejected, the same safety net as in Metropolis. During warmup a Nesterov dual-averaging scheme adapts the step size toward a target acceptance rate of , then freezes it to the running average for the sampling phase. This is the adaptation recipe from Hoffman and Gelman’s NUTS paper, applied to plain HMC.

The variational engine fits a factorized Normal by maximizing the ELBO. It stores the mean and the log standard deviation (so the standard deviation stays positive for free). Each iteration draws samples from , and for each sample it computes the gradient of at . The reparameterization trick moves the gradient of the expectation onto and :

The on the log-standard-deviation gradient is the entropy term : it widens unless the data pulls the mean toward higher . Adam ascends both parameter sets.

After Adam converges, the engine draws a large sample from the fitted and returns those draws, transformed back to the natural scale, as its posterior summary. Variational inference is fast and deterministic, but its answer is only as good as the factorized Normal approximation. For a unimodal posterior on the unconstrained scale it does well; for a posterior with strong correlations or multiple modes it can miss.

Sampling only approximates the posterior, so the file reports two standard diagnostics alongside each fit.

The effective sample size (ESS) measures how many independent draws the autocorrelated chain is worth. If consecutive draws are correlated, the chain carries less information than its length suggests. Geyer’s initial positive sequence estimator sums the autocorrelations until they turn negative:

where is the lag- autocorrelation. An ESS near means the draws are nearly independent; an ESS far below means the chain mixes poorly and you should run longer or switch engines.

The Gelman-Rubin (R-hat) compares the variance within each chain to the variance between chains. Run several chains from dispersed starts. If they all sample the same posterior, the between-chain variance should match the within-chain variance and . If the chains disagree, rises above , signalling that the chains have not converged. The statistic is

with the between-chain variance and the average within-chain variance. A common rule of thumb is to trust the fit when for every parameter. R-hat needs at least two chains; with one chain the file reports ESS only.

Load the file and it drops you into the PPL REPL:

Type demo at the ppl> prompt for a guided tour, or fit a model directly. The (demo) function fits the coin and the Normal models, plots their posteriors, and compares all three engines against the exact Beta answer. Because the random state is freshly seeded each session, your numbers will differ slightly from those below, but the patterns hold.

The coin model puts a prior (uniform on ) on the bias and scores flips, of them heads, with a Bernoulli likelihood:

With a prior and heads, tails, the exact posterior is with mean (Chapter 9’s conjugate update). HMC with two chains reproduces it:

The posterior mean sits next to the exact . The 95% credible interval is what the conjugate gives. The histogram is the shape of that Beta density: roughly symmetric, peaked near , tapering to both sides. ESS equals the full draws and R-hat is , so the two chains agree and mix almost perfectly. The acceptance rate is right at the HMC target of .

The demo fits the same coin with each engine and checks the posterior mean against the exact :

All three land within a couple of percentage points of . Metropolis and HMC, both sampling the posterior directly, are essentially exact up to Monte Carlo noise. Variational inference is less accurate here because the factorized Normal is fit in unconstrained space and then pushed through the sigmoid; the mean of a transformed Normal is not the transform of the mean, so the bias shifts slightly. This is the usual cost of variational inference: speed and determinism in exchange for a small, controllable bias.

The mean-var model puts a prior on the mean, a Half-Normal prior on the standard deviation, and scores draws from with a Normal likelihood:

Two chains recover both parameters:

The posterior mean of is , bracketing the true . The posterior mean of is : with only draws the sample standard deviation has upward scatter, and the half-normal prior adds a gentle pull, so the estimate sits a little above the true . The credible intervals cover both true values. R-hat is for both and ESS is high, so the chains agree and mix well.

The regression model puts Normal priors on the intercept and slope and a Half-Normal prior on the noise, then scores points. The mean is built with the differentiable operators and , so HMC and VI get gradients through it:

Two chains recover the true line with noise :

The intercept and slope both cover their true values of and , and the noise covers the true . This is the first model with a linear predictor: the likelihood mean depends on two latents at once. The AD pass threads the gradient through and without the model author doing anything special, which is what makes the DSL usable for models more complex than a single prior.

The same coin fit by random-walk Metropolis shows why diagnostics matter. A standalone run gives:

The posterior mean is right on the exact , but note the cost. Metropolis used iterations to get ESS , while HMC got ESS from only draws. Each Metropolis draw is worth roughly a fifth of an independent draw; each HMC draw is worth nearly one. The autocorrelation plot makes the cause visible:

The autocorrelation decays slowly: even at lag it is still positive. The area under this curve is the in the ESS denominator, and that area is large, so ESS is far below . HMC’s leapfrog moves slide along the posterior and decorrelate quickly, which is why its autocorrelation drops to near zero within a few lags and its ESS sits near . R-hat is for both engines, so both have converged; the difference is purely efficiency. This is the practical case for HMC over random-walk Metropolis in even modest dimensions.

This chapter assembled a working probabilistic programming language from three ingredients: a tiny DSL (defmodel, sample, observe) that lets you state a model as plain Lisp; a layer of forward-mode automatic differentiation that turns a log-density function into its gradient; and three inference engines that share the same model and the same unnormalized log posterior. Metropolis-Hastings is the simple, gradient-free baseline. HMC adds a gradient and a momentum and mixes far better. Variational inference trades exact sampling for a fast, deterministic approximation.

The whole system fits in one file of roughly lines and uses only Common Lisp primitives plus a few SBCL-specific float-trap masks. It connects the earlier chapters into a single tool: the Normal, Beta, Gamma, and Poisson distributions from Chapters 5 and 9, the Monte Carlo idea from Chapter 8, and the Markov chain theory from Chapter 10 all reappear here. The Metropolis-Hastings chain is a Markov chain whose stationary distribution is the posterior, and the convergence diagnostics R-hat and ESS measure exactly the mixing behavior Chapter 10 analyzed.

The DSL is extensible. Adding a new distribution means writing one call with a log density and a sampler. Adding a new model means writing one form. To go further you might add a Student-t likelihood for robust regression, a log-normal prior, a No-U-Turn Sampler to tune HMC’s path length automatically, or a reparameterization that handles posteriors with strong correlations. The framework is already there; the model and the engine are finally separate.

Problem 11.1. For the coin model with a prior and heads in flips, the exact posterior is with mean and variance . Compute the exact posterior standard deviation and compare it with the column the HMC summary reports (about ).

Problem 11.2 (Two modes). Trace through the coin model body in both and modes. In mode, what value does sample return, and what gets pushed onto the context? In mode with parameter vector , write down the three terms that sample adds to ctx-logdens (the Beta prior log density, the Jacobian, and the Bernoulli likelihood is added later by observe). Why does observe do nothing in mode?

Problem 11.3 (ESS and efficiency). The Metropolis run used iterations to reach ESS , while HMC used draws to reach ESS . How many Metropolis iterations would you need to match HMC’s effective draws, assuming the same efficiency? What does this say about the per-draw cost of random-walk Metropolis versus HMC?

Problem 11.4 (R-hat). Run the coin model with and read the R-hat column. Now deliberately break mixing by setting the HMC step size very large (:step 2.0d0) with few leapfrog steps. What happens to the acceptance rate, the ESS, and R-hat? Explain how R-hat detects the problem even when ESS alone might look acceptable.

Problem 11.5 (The Jacobian). A standard deviation has support, so the DSL uses . Show that the log Jacobian is , as the code claims. Now suppose a probability has support with . Derive the log Jacobian and confirm it matches log-jac-support. What would go wrong if the sampler moved in the natural space and ignored the Jacobian?

Problem 11.6 (Adding a distribution). Add a distribution to the library. Its density on the positive axis is for . Write the call using the differentiable operators, choosing the right :support. What is the Jacobian contribution, and why is it the same as for any latent?

Problem 11.7 (A new model). Using the existing distributions and the DSL, write a for a Poisson rate: prior, and a vector of count data with (poisson lambda). Generate counts from a Poisson with true rate , fit the model with HMC and two chains, and check that the posterior mean of covers . Compare with the exact Gamma posterior from Chapter 9’s conjugacy (Gamma-Poisson).

Problem 11.8 (Comparing engines). Fit the regression model with all three engines (:mh, :hmc, :vi), two chains each where supported. For each engine, report the posterior mean of , the ESS, and the wall-clock time (use around the call). Rank the engines by accuracy and by speed. Which engine gives the best tradeoff for this three-parameter model?

Problem 11.9 (Variational bias). The three-engine table shows variational inference with a larger error than the two samplers. The variational guide is a factorized Normal in unconstrained space, transformed by the sigmoid. Explain in your own words why the mean of with is not , and why this introduces a bias that grows with the variance . Suggest one way to reduce the bias (for example, fitting on the natural scale, or using a Beta guide).

Problem 11.10 (Coding exercise). The HMC engine uses a fixed, jittered path length. Replace it with a simple doubling scheme: start with one leapfrog step, and double the path length while the proposed Hamiltonian energy decreases, stopping at a maximum or when the energy starts to rise. This is a simplified version of the No-U-Turn condition. Test your modified engine on the regression model and report the effect on ESS and acceptance rate compared with the fixed-path-length version.