If you have used core.logic, you already know half of probabilistic programming. A logic program states constraints and lets the system search for the values that satisfy them. A probabilistic program does the same, except that a constraint can be soft: a value that fits the data poorly is not ruled out, only made less likely. Two dice show the idea below. A small forecasting model then runs under several inference methods in Foerster, an MIT-licensed probabilistic programming library for Clojure, and every result is checked against an exact answer.
Foerster continues the Anglican line of probabilistic programming in Clojure. Christian Weilbach, who wrote it, did his PhD in Frank Wood’s group at the University of British Columbia, where Anglican was built.
What is a probabilistic program? #
Take two dice and the constraint that they sum to 7. Without any library, the plainest way to find the pairs that meet it is to generate pairs at random and keep the ones that pass:
The random choices form a generative model, a small simulator of two dice. The test conditions the simulator on what we know. What survives is the posterior: the pairs that are still possible, and how often each occurs.
A probabilistic program states the same model without fixing how to search it:
sample introduces an unknown, as a fresh logic variable does in core.logic, but with a probability on each value: here a die whose six faces are equally likely. factor adds to the run’s log weight. Zero leaves the run unchanged, and minus infinity makes it impossible, as a failing goal does. sp/spin turns the body into a program that an inference method can , copy and resume at each sample and factor. With {:method :enumerate}, Foerster runs all 36 pairs and returns the six that hold, each with probability 1/6. That is the answer run* gives, with weights.
Now soften the constraint. We do not know the sum, but a sensor that measures it with noise of about one showed 7:
observe weights each run by how likely the reading is under it, here a normal density centred on the run’s sum. No pair is ruled out, and the arithmetic can be done by hand. Each of the six pairs that sum to 7 gets weight e⁰ = 1, each of the five pairs that sum to 6 or to 8 gets e^−½ ≈ 0.61, and the weights keep falling with the distance from 7. Normalized, a sum of 7 has probability 0.45, sums of 6 and 8 have 0.23 each, sums of 5 and 9 have 0.04 each, and the rest almost nothing. Foerster’s enumeration gives the same numbers.
| core.logic | Foerster |
|---|---|
| a fresh logic variable | sample , an unknown with a probability on each value |
a goal such as == , which holds or fails |
observe orfactor , a weight; failing is weight zero |
conde , a choice between branches |
a discrete sample , with a probability on each branch |
run* , the set of answers |
the posterior, the answers with their weights |
| the search strategy | the inference method |
Two more ideas carry over. The total weight before normalizing, called the evidence, says how well the model as a whole explains the data. A hard constraint that no run meets leaves nothing, and a soft one leaves a number that can be compared between models. And as core.logic leaves the search to the system, a probabilistic program leaves the inference method open. ProbLog joins the two traditions as probabilistic logic programming.
How does the system find the answer? #
Enumeration works for 36 pairs. It fails for many choices, or for continuous ones, because no list of all runs exists. The other methods search the space of runs, and each resembles a familiar strategy:
- Importance sampling runs the program many times and weights each run instead of discarding it. It is generate and test with soft scores, and it fails when good runs are rare.
- Sequential Monte Carlo (SMC) runs many copies side by side. At each observation it copies the runs that fit and drops the ones that do not, like a search that prunes partial solutions as soon as a constraint can be checked.
- Markov chain Monte Carlo (MCMC) changes one choice of a single run at a time and keeps the change by a rule that sometimes accepts a worse run. This local search visits runs in proportion to their weight.
Does the choice of method matter? #
A shop sells around a hundred units of a product each week. The weekly figure is noisy, and the underlying demand drifts. How high is demand now, and how certain is that? The model says that demand starts near 100, moves by a random amount each week, and that each week’s figure is demand plus noise:
As with the dice, sample marks the unknowns, here twenty-one weekly demand levels, and observe scores each run against the twenty sales figures. The data are an ordinary argument, weeks. The :id names each choice, so a method can refer to the level in week 7 across runs. The pictures article explains priors and posteriors if those terms are new. Each method is one option of the same call:
The model is linear and Gaussian, so a Kalman filter computes the exact answer, which makes it a fair test. The data are twenty simulated weeks, with a true drift of 4 per week and observation noise of 8. Each particle method ran ten times with different seeds, so the table shows how much an answer varies from run to run. All numbers come from Foerster’s main branch on 15 October 2026.
| Method | Final level (mean) | Log evidence | Check |
|---|---|---|---|
| Exact (Kalman filter) | 102.89 | −77.49 | |
| Sequential Monte Carlo, 1000 particles | 102.92 ± 0.23 | −77.44 ± 0.14 | about 120 histories survive to week 0 |
| Importance sampling, 1000 samples | 102.9 ± 1.8 | −77.7 ± 0.5 | 2 to 10 effective samples |
| Metropolis–Hastings, 4 chains | 103.4 | R-hat 1.03, 114 effective draws | |
| Particle Gibbs, 4 chains | 102.93 | R-hat 1.002, about 1500 effective draws |
SMC lands within 0.36 of the exact final level in every run, and its log evidence averages −77.44 against the exact −77.49. It pays for its pruning on the early weeks: of a thousand particles at the end, only 114 to 134 still descend from different starting points. Importance sampling never prunes, so after twenty weeks a handful of runs carry nearly all the weight, and the estimate moves by almost two units between runs.
Metropolis–Hastings (MH), the simplest MCMC, fails here, and its diagnostics say so before any exact answer is consulted. Its four chains disagree: their rank-normalized R-hat is 1.03, above the usual threshold of 1.01, and their 12,000 draws are worth about 114 independent ones. Foerster computes these numbers as Stan and ArviZ do. The cause lies in the model. Neighbouring weeks are strongly tied, so changing one week’s level at a time barely moves the run. Andrew Gelman’s folk theorem of statistical computing, restated in the Bayesian workflow paper, says as much: when the computation struggles, look at the model.
The fix is a method that changes all weeks together. Particle Gibbs runs a whole SMC through the twenty weeks at every step, holding one surviving run fixed, and so redraws all twenty-one levels at once. Four chains of 500 steps after warm-up, with 50 particles each, reached an R-hat of 1.002 and about 1500 effective draws of the final level, counting each step’s weighted average as one draw. Their mean was 102.93 against the exact 102.89. The first week’s level, which SMC sees through only a few surviving histories, reached an R-hat of 1.002 and about 660 effective draws. Particle Gibbs took four minutes against half a minute for MH, but produced about twice as many effective draws per second. Nothing in the model changed, only the search.
What if a parameter is unknown too? #
So far the drift was given as 4. In practice nobody measures how fast demand moves, so the drift becomes one more unknown with a prior:
Particle marginal Metropolis–Hastings (PMMH) searches over the drift alone. For each proposed drift it runs SMC over the weekly levels and accepts or rejects the drift by the evidence that SMC reports. That evidence estimate is unbiased, so the chain targets the exact posterior of the drift even though every step uses an estimate.
The chain’s mean drift is 2.82 against an exact 2.85, and its 90% interval runs from 0.61 to 6.09 against an exact 0.61 to 6.28. The data were simulated with a drift of 4, inside the interval but not at its centre, because twenty noisy weeks cannot pin the drift down. A forecast with one fitted drift would hide that width. The 2000 iterations took five and a half minutes on a laptop. For this linear Gaussian model, Stan or a Kalman filter inside a one-parameter chain would take seconds. PMMH is for models without a gradient to follow, models that branch, count or call a simulator.
How do I run it? #
Foerster is on Clojars. With the Clojure CLI:
{:deps {org.replikativ/foerster {:mvn/version "0.1.113"}}}
A model runs inside a Spindel world, an execution context, and inference returns a spin that the REPL dereferences:
Build the model inside with-context, as here: outside a world, the program has nowhere to run. The scripts that produced every number above are in the Foerster repository.
What can Foerster not do yet? #
Foerster is beta software. Beyond the methods shown here, it implements tempered and waste-free SMC, SMC², involutive MCMC with reversible jump, parallel tempering, the No-U-Turn sampler on numerical blocks, MAP estimation, the Laplace approximation, Pathfinder, and variational inference. The methods notebook describes each one with its papers, and the gallery works through classic models against exact answers.
For a static model with continuous parameters, Stan, PyMC and NumPyro remain the mature choice. Foerster’s No-U-Turn sampler works on numerical blocks whose gradient the model supplies or Raster compiles. On one hierarchical model with 89 parameters, it matched PyMC’s warm run with a gradient compiled by Raster, and ran 1.5 times faster with a hand-written one. Its general methods over program traces are slower: single-site MH takes about twelve times as long per move as Gen, because each move replays the rest of the program, and its particle filter is about four and a half times slower than Gen’s Unfold combinator. The comparison has the details, and coming from other systems maps Stan, PyMC, Turing, Gen, Anglican and WebPPL concepts onto Foerster, including what is missing.
What comes next? #
Every model so far is a pure function. Useful models often are not: a forecast that reads current stock, a simulation that writes into a scenario database, an agent that calls tools. If a thousand particles write to the same database, inference corrupts the state it reasons about, which is why Gen requires that a model change no external state.
Foerster runs each particle in a Spindel world, whose state can be copied cheaply, as the copy-on-write article describes. When SMC copies a particle, it copies the particle’s world, and when it drops a particle, it discards the world. With :world-policy :fork, every particle starts from a copy of the caller’s world, so a model can read and change the systems registered there, and nothing a particle writes reaches the caller. This costs about three times as much as running in fresh worlds. Systems that cannot be copied safely, such as live network handles, are refused before the model runs, and Markov chains run in fresh worlds only. The same mechanism makes counterfactuals direct: intervene on one choice in the world where the posterior was computed (Pearl’s do operator, as in the causal graph article) and replay the program with everything else held to what was inferred. The worlds guide covers both, and a later article shows them on a model that books orders into a forked database.
The same machinery steers processes that are not statistical models. Foerster’s steer namespace runs SMC over the steps of an agent working on a task, following twisted SMC (Zhao et al., 2024) and the SMC steering of language models (Lew et al., 2023). On short spreadsheet tasks it matched picking the best of four attempts with 17% fewer tokens, but solved nothing more. It should pay on long tasks whose steps can be checked cheaply.
For Simmis, inference over forked worlds is research rather than a product feature. It is how a governed organizational record could later answer a question such as “what happens to our cash if this order slips” with a distribution instead of a single number.