A reference for every inference entry point in org.replikativ.foerster.core
(infer below). Each returns a spin resolving an EmpiricalMeasure (see
posteriors); at the REPL deref it inside
sp/with-context, inside a spin await it. For which to use when, see the
algorithms notebook.
infer/infer runs any method by name, as Anglican's doquery did:
(infer/infer (model) {:method :smc :particles 1000})
(infer/infer (model) {:method :mh :iterations 4000 :chains 4 :burn 1000})
(infer/infer (model) {:method :pmmh :particles 100 :iterations 2000 :params #{:sigma}})
Methods: :enumerate (exact, see below); :importance, :smc, :tempered (:particles);
:parallel-tempering (:chains, one block site); :pimh,
:pgibbs, :pgas, :ipmcmc, :bbvi (:particles, :iterations); :mh,
:rmh, :nuts (:iterations per chain, :chains, :burn; :step-size for
:rmh and :nuts; :target-accept, :max-depth, :metric,
:inverse-metric, :jitter for :nuts); :kernel
(:kernel, :chains; the kernel carries its own :burn and output
options); :pmmh (:particles, :iterations, :burn, :params). A size
option a method does not take is refused; other options go to the method's
own function, described below. :seed makes a run reproducible without
touching the process generator, for the particle methods (importance
sampling, SMC, particle MCMC, BBVI) and the Markov chains alike: chain i is seeded by (derive-seed seed ::chain i),
so several sessions sharing a JVM do not disturb each other's draws.
The papers behind each method are in literature, and the
inference methods notebook
runs each on a model with a known answer.
Which to start with:
:smc with as many particles as you can afford; it also
estimates the evidence.:nuts on a
numerical block (:rmh on ordinary sites), several chains, and check
R-hat and ESS
(posteriors).:pmmh, or SMC² online;
a multimodal posterior: :tempered, or :parallel-tempering for a
model that is one block.optimize/pathfinder or optimize/laplace on a block site.:enumerate, the reference
to check the others against.SMC and tempered SMC spread their particles over the executor's threads (about one worker per processor); the results per seed do not depend on how many threads ran them.
{:method :enumerate} (org.replikativ.foerster.inference.enumerate) computes the
exact posterior of a program whose latent choices all have finite support
(dist/support: Bernoulli, flip, discrete, categorical, binomial,
uniform-discrete): at every latent site the particle's world forks once per
value, each branch weighted by that value's probability. Branches share the
program up to where they part, so the cost is the number of complete
assignments (:max-branches, default 100 000, bounds it). Its measure has
one particle per assignment and its m/log-marginal is the exact evidence:
the reference for checking the other methods on small discrete models. A
latent site with infinite or continuous support is refused.
smc-infer, importance-sampling, pimh-infer, pgibbs-infer,
pgas-infer, ipmcmc-infer, bbvi-infer, and kernel-infer with a
non-Markov-chain kernel:
| Option | Default | Meaning |
|---|---|---|
:world-policy | :fresh | :fork runs the particles in canonical forks of the caller's world (worlds) |
:world-opts | {} | fork options for :fork: :systems (:all, :none or a set of system ids), :rights, :snapshots |
:authority, :grant | — | budgets under :fork (worlds) |
:executor | spindel's default | the executor the worlds run on |
:resample-threshold | 0.5 | SMC resamples when the effective sample size is below this fraction of the particles |
:policy | the prior | a foerster.runtime.trace/policy deciding the sites: constraints, interventions, custom proposals (extending) |
:anchors, :rejuvenate | — | resample-move (below); smc-infer and pimh-infer, fresh worlds |
:resampling | :systematic | :stratified, :residual or :multinomial (Douc, Cappé & Moulines 2005); smc-infer, pimh-infer |
:genealogy? | false | record every resampling's ancestor indices in the measure's :history; smc-infer, pimh-infer |
:batch | — | B: resample in batches of B in arrival order instead of at a barrier for all (below); smc-infer, fresh worlds |
:smcp3 | — | {:forward K :backward L}: SMCP3 move-reweight steps at each observation (below); needs :anchors; smc-infer, pimh-infer, fresh worlds |
An option a method does not take is refused (::infer/unknown-options), and
so are :world-opts, :authority and :grant under :world-policy :fresh
(::infer/fork-only-options): an option that would be ignored is a mistake.
importance-sampling takes no :resample-threshold, BBVI no :policy or
:resample-threshold (it decides both), and the Markov-chain kernels only
:executor, :seed, :world-policy and :world-opts.
A particle method draws latent sites from the prior unless given a
:proposal (proposals).
A program can also name a proposal at a site:
(sample (dist/normal 0 1) :id :mu :proposal guide) draws a fresh value
from guide — a distribution the program computes, for instance from its
data, which is how an amortized guide plugs in — and weighs it by
log p − log q, so every method's target stays the program's posterior.
Replays, moves and a policy's :draw take precedence.
| Function | What it does |
|---|---|
(smc-infer model n opts) | Sequential Monte Carlo with n particles, resampling at observations. The measure's m/log-marginal estimates the evidence. |
(importance-sampling model n opts) | n weighted runs, never resampled. |
(pimh-infer model n iterations opts) | Particle independent MH: iterations SMC sweeps of n particles, each accepted on its evidence estimate; every sweep's particles pooled. |
(pgibbs-infer model n iterations opts) | Particle Gibbs: iterated conditional SMC, each sweep keeping a trajectory drawn from the previous one; sweeps pooled. |
(pgas-infer model n iterations opts) | PGAS: particle Gibbs whose retained particle redraws its ancestor at every observation. Costs a replay of the retained future per particle per observation. |
(ipmcmc-infer model n iterations opts) | Interacting particle MCMC: :num-nodes (default 8) SMC nodes, :num-csmc-nodes (default half) of them conditional, exchanged by Gibbs updates on their evidence; :all-particles? (default true) emits every node's particles, false one per node. |
(bbvi-infer model n iterations opts) | Black-box variational inference: learns a mean-field q in iterations updates from n samples each. Each site's q starts as its prior when first reached, then moves by gradient only (it does not follow a prior that depends on other latents). :base-lr (1.0), :robbins-monro (0.0; step size base-lr / (t+1)^robbins-monro), :adagrad (true: AdaGrad; γ ∈ (0,1): RMSprop; false: plain steps). The measure holds n importance samples from the final q; (infer/get-variational-dists measure) returns q, {address distribution}. |
The pooled particle-MCMC measures are MCMC estimates: their weights are
normalized per sweep, and m/log-marginal of the pool is not an evidence
estimate.
(tempered-infer model n opts) (org.replikativ.foerster.inference.tempering) moves a
population of complete program runs from the prior to the posterior through
the targets p(x)·L(x)^β, L the observations' and factors' likelihood, with β
rising from 0 to 1 (an SMC sampler; Del Moral, Doucet & Jasra 2006). Each
step chooses the next β so that the conditional ESS of the reweighting stays
at :ess-target·N (Zhou, Johansen & Aston 2016), resamples, and moves every
particle by Metropolis-Hastings at the new temperature — random walks scaled
by the population's spread on continuous sites, prior proposals elsewhere.
Use it for posteriors with separated modes, which a Markov chain does not cross and SMC in program order only reaches if its early particles happen to:
(infer/tempered-infer (model) 1000) ; β schedule in :temperatures
(infer/tempered-infer (model) 1000 {:waste-free 10}) ; Dau & Chopin 2022
| Option | Default | |
|---|---|---|
:ess-target | 0.5 | higher: more, smaller temperature steps |
:moves | one sweep | single-site MH moves per particle per step |
:scale | 2.38 | random-walk scale in population standard deviations |
:waste-free | — | P: resample N/P particles and keep every state of their P-step chains |
:max-steps | 1000 |
m/log-marginal estimates the evidence. Every move re-runs the program from
the moved site, so a step costs N sweeps of the model; fresh worlds only.
Block sites are tempered too. A block keeps its observations inside its
target, so a block site with a draw density q (its :sample-log-density)
follows the geometric path q^(1−β)·t^β from q to its target t, and is moved
by random walks on its θ vector scaled by the population's spread at the
step's temperature (foerster.moves.block-mh). A block without
:sample-log-density is refused (::block/no-sample-density). A model
that is one block site and nothing else never makes a world per
particle: the population is a set of θ vectors that are reweighted,
resampled and moved directly (foerster.population.theta), and only the distinct
final θ are written back to traces, by one replay each. On golf putting as
one block, 1000 particles take about a second and estimate log Z within
0.2 of the exact −191.93.
:kernel moves the particles with a Markov kernel of your choice instead
of the built-in moves (one iteration per particle per step): any fixed
kernel that targets temperatures below one — single-site MH (:only to
restrict it), finite and block Gibbs, their cycle and mixture, or a
kernel of your own declaring :tempering?. NUTS and HMC target β = 1 only
and are refused, as are kernels that adapt. With a kernel, a one-block model
runs on traces too.
For data that arrive one observation at a time — IBIS (Chopin 2002) — use SMC with resample-move on the static parameters, below.
(parallel-tempering-infer model chains opts)
(org.replikativ.foerster.inference.parallel-tempering) runs non-reversible parallel
tempering (Syed, Bouchard-Côté, Deligiannidis & Doucet 2022) on a model that
is one block site with a :sample and :sample-log-density. chains
Markov chains sit on the path q^(1−β)·t^β from that draw density q to the
posterior t; every scan moves each chain at its own β (an exact draw at
β = 0, random-walk Metropolis elsewhere) and then offers swaps to neighbours,
even pairs and odd pairs in turn, so states travel through the temperatures
in one direction instead of diffusing. The schedule is tuned over :rounds
rounds of doubling length, spacing the temperatures evenly in the
communication barrier Λ the swap rejections measure.
(infer/infer (model) {:method :parallel-tempering :chains 10 :rounds 11})
The last round's states of the β = 1 chain are the sample (2^(rounds−1) of
them, written back to traces once per distinct θ); m/log-marginal is the
stepping-stone evidence. :barrier reports Λ: about 2Λ chains is enough,
and :round-trips counts the states that went from the reference to the
posterior and back. Where tempered SMC spends its budget on many particles,
parallel tempering spends it on long chains, and its result improves with
more rounds without rerunning.
(scmc/scmc model opts) (org.replikativ.foerster.inference.scmc) is SMC over data
that arrive in batches, for a model given as a prior and a per-datum
log-likelihood on θ vectors (Beronov, Weilbach, Wood & Campbell 2021).
Resample-move SMC over all past data (IBIS) must score every past
observation in each move, so its step t costs O(t) and the run O(T²).
SCMC keeps at most :memory past points with weights, a core-set whose
weighted log-likelihood approximates that of all data so far. The particles
move against the core-set. After each step, the core-set plus the new batch
is recompressed with GIGA (Campbell & Broderick 2018), fitted where the
current particles are. GIGA's weights are then refitted by nonnegative
least squares on the points it chose, the fully corrective update of
Campbell & Beronov (2019) (:refit :nnls, the default; :refit :none keeps
GIGA's weights, as in the UAI 2021 algorithm):
(def posterior
(scmc/scmc {:prior-sample (fn [] [a b])
:prior-log-density (fn [[a b]] …)
:log-lik (fn [[a b] [x y]] …) ; one datum
:batches (partition 20 data)}
{:particles 500 :memory 50}))
(mapv (comp m/get-value first) (m/get-particles posterior)) ; θ vectors
(m/log-marginal posterior) ; evidence estimate
(:coreset posterior) ; [[w datum] …]
The synchronous numeric API returns the same EmpiricalMeasure of Samples
as the program-based inference APIs, with equal particle weights. Its
:inference-contract, :coreset, :steps, :max-retained-coreset-size
(at most :memory), :max-transient-coreset-size (core-set plus batch;
:max-coreset-size is its old name), :moves, :accepted and
:last-update (ESS before resampling and likelihood evaluations by phase)
fields describe approximation and execution.
Each batch reweights the particles by its exact likelihood, but the
particles were moved against the core-set before the last compression, and
nothing corrects for the switch to the compressed one (as in the paper's
Algorithm 1); the moves then target the compressed core-set. The result
therefore carries the core-set's approximation error, which is small when
the log-likelihoods compress well (few parameters, exchangeable data). With
:memory at least the amount of data, nothing is compressed and SCMC is
resample-move SMC.
A step evaluates each likelihood once per state: the batch at every
particle for the weights, and the expanded core-set at each proposal of the
moves. The values at the particles are kept in the state, so a move's start
state and the compression cost no new evaluations; they are used only with
the log-likelihood, core-set and particles they were computed for, and
re-evaluated otherwise. :projection-inflation k
is a research option, off by default and a deviation from the published
algorithm: it fits the compression at the particles spread by k about their
mean, μ + k·(θ − μ), at N·|E| more evaluations per compression.
For streaming use, scmc/start, scmc/append-batch and scmc/result keep
the state as a persistent value. start records its seed (drawn when none
is given); (scmc/initial state) and (scmc/options state) restart the
same run, so a replay after a historical correction is reproducible.
On a two-parameter Bayesian regression in batches of 20, with 500 particles:
| data | :memory 50 | all data | posterior error (:memory 50) |
|---|---|---|---|
| 1000 | 1.2 s | 3.5 s | means within 0.08 sd, sds within 2% |
| 2000 | 2.4 s | 12.6 s | within 0.03 sd and 4% |
| 4000 | 4.7 s | 48.7 s | within 0.02 sd and 3% |
scmc/giga is available on its own: nonnegative, at most m-sparse weights
over a set of vectors that approximate a target vector.
SMC resamples at each observation once every particle has arrived, so a
population whose steps take uneven time — a model call, a tool, a simulation
that sometimes runs long — waits for its slowest particle at every
observation. With :batch B, the particles arriving at an observation form
batches of B in arrival order; a full batch resamples within itself at once
(B children, each weighted by the batch's mean weight), and the last batch
of an observation closes once no other particle can still reach it.
(infer/smc-infer (model) 1000 {:batch 32})
Weights are never reset, so the evidence is the final mean weight. Each
batch's resampling keeps its total weight in expectation, whatever decided
which particles are in it, so the estimate stays unbiased even when the
arrival order depends on the particles' states (Paige et al. 2014's
condition for the particle cascade, here per batch; an island model with
islands formed by arrival). B = N is ordinary SMC; smaller B resamples
sooner but within fewer particles. The measure's :history lists every
batch. Batch membership follows arrival order, so on a multi-threaded
executor runs are not reproducible from their seed. Not with stream sites,
:retained, :anchors or :smcp3.
(cascade/cascade model n opts) (org.replikativ.foerster.inference.cascade; Paige,
Wood, Doucet & Teh 2014) runs every particle on its own. At an observation a
particle compares its weight W with the running mean W̄ of the weights that
arrived there so far. Below the mean it survives with probability W/W̄, as
one copy weighted W̄; above it, it branches into ⌈W/W̄⌉ copies (forks of its
world) while the stage has produced no more children than min(n, arrivals),
⌊W/W̄⌋ otherwise, each weighted W/M — the paper's Eq. 14, whose feedback
keeps the population near n. Nobody waits for anybody; the evidence
estimate (1/n) Σ W over the particles that reach the end is unbiased
whatever order they arrive in.
| Option | Default | |
|---|---|---|
:cap | 4·n | most particles alive at once; copies beyond it fold into a multiplicity on the particle |
:concurrency | 32 | particles running between observations at once; the rest wait and resume in random order, which keeps the population stable (the paper's §4) |
:policy, :executor |
It ends when every particle has finished, never by wall clock, which would
favour the fast ones (Murray, Singh & Lee 2021). Each particle has one child
per observation on average, so the number that reach the end varies a lot
from run to run (a critical branching process); its :cascade reports
:launched, :finished, :peak and :collapsed. Runs on a
multi-threaded executor are not reproducible from their seed. Where the
population should stay fixed, :batch (above) is the steadier choice.
SMC resamples by copying good particles and dropping bad ones, so after a few observations most particles share their early choices (a static parameter sampled once is the extreme case). Resample-move (Gilks & Berzuini 2001) gives every particle Metropolis-Hastings moves after each resampling, which restores diversity without changing the target or the weights:
(infer/smc-infer (model) 1000
{:anchors #{:mu} ; where moves may start
:rejuvenate {:moves 2 ; per particle, per resampling
:propose (proposal/random-walk 0.2)}})
Any fixed Markov kernel can be the move instead (:rejuvenate {:moves m :kernel k}): a cycle of MH and random-walk moves, HMC on a block, a
kernel of your own. Kernels that tune themselves (:burn > 0) are refused
(warm-up is a separate chain), and so is :kernel with :retained
(particle Gibbs draws its retained state backwards through the move, which
needs a reversible move; use :propose). The evidence stays :unbiased
with built-in kernels and is :unverified with a kernel of your own.
A move replays the program from an anchor — a kept fork of the world at a latent site — up to the observation the particle is parked at, and accepts on the MH ratio of the two partial traces. Which sites keep anchors is the trade-off:
:anchors | Moves may target | Cost |
|---|---|---|
| a set of addresses | those sites (static parameters) | one anchor per particle per site; a move replays from the site to the current observation, O(t) |
{:lag L} | latent sites of the last L observations | anchors of older sites are released; a move replays O(L) steps — fixed-lag rejuvenation, and the choice for smc/stream |
:all | every latent site | an anchor per latent site, shared between particles of one ancestor |
(fn [sp]) | sites the predicate names |
Moves happen only when the population is resampled, so with the default
:resample-threshold they become rare as a static posterior concentrates
(Chopin 2002). The measure's :rejuvenation reports :moves, :accepted
and :max-anchors. Resample-move runs in fresh worlds only. A model whose
moves reach effects (a model call) runs them again on every replay.
In particle Gibbs (:retained, pgibbs-infer with :anchors and
:rejuvenate) the retained particle moves too. Leaving it unmoved while the
others move puts the sweep's stationary law slightly off the posterior (an
exact enumeration of a three-step chain: 1.7·10⁻⁴ in total variation, far
below what a statistical test sees). Instead its trajectory is taken as the
state after each barrier's move, and its state before the move is drawn
backwards through the move: the moves are reversible Metropolis–Hastings
kernels, so that is the move itself, run from the last barrier back to the
first before the sweep starts. The retained particle enters each barrier
with its pre-move values and leaves it with the post-move ones. The same
enumeration puts this construction at rounding error (clojure -M:test:enumerate). Not with ancestor sampling (PGAS).
SMCP3 (Lew et al. 2023) generalizes resample-move: at each observation a particle takes a step of a forward kernel K and is reweighted by a backward kernel L instead of accepting or rejecting. K may change any latent site the particle has passed — revise the past, propose from the data just seen — and the weight keeps the particle properly weighted:
w ← w · p(x')·L(u' | x') / (p(x)·K(u | x)) · |J|
(def random-walk-kernel
{:forward (fn [{mu :mu} trace] ; the particle's choices and trace
(let [u (dist/draw (dist/normal 0.0 0.05))]
{:updates {:mu (+ mu u)} ; new values
:log-q (dist/logpdf (dist/normal 0.0 0.05) u)
:reverse (- u) ; what L would draw to go back
:log-jacobian 0.0})) ; a random walk preserves volume
:backward (fn [choices' u'] (dist/logpdf (dist/normal 0.0 0.05) u'))})
(infer/smc-infer (model) 1000 {:anchors #{:mu} :smcp3 random-walk-kernel})
K draws from foerster.random (it runs in the particle's stream), must
return :log-jacobian (the log |det| of its deterministic map; 0.0 for a
random walk or a discrete move), and needs an
anchor at the earliest site it updates. Sites it does not update are kept,
and sites the replay reaches afresh are drawn from their prior. The weight's
variance depends on how well L reverses K: a symmetric random walk as L
multiplies the weights by π(x')/π(x), which is fine for steps small against
the posterior's spread and heavy-tailed for large ones; an L close to the
reversal of K under the target keeps them flat. Fresh worlds, not with
:retained.
The {:forward :backward} map is a pair (foerster.moves.pair), the one
shape for moves with auxiliary randomness: SMCP3 takes it, a population
transport takes it (program/transport-kernel source destination pair),
and an involutive MH move becomes one with pair/from-involution, so the
same move serves as an MCMC move (k/involutive-kernel) and as an SMCP3
step (a reversible jump, which names :removed sites, stays an MCMC move).
A program whose proposals come from an outside process — a language model's
turn, a simulator — with no density to report, and whose data is a score
(a verifier, a judge, a process reward model), is steered by factors that
are barriers: (factor w :barrier true) makes SMC park and resample there,
as at an observation. foerster.inference.steer/model builds the step loop:
(infer/smc-infer
(steer/model {:init s0
:step (fn [state] (spin …)) ; the next state, from the process
:value (fn [state] (spin …)) ; log ψ: the reward to come
:reward (fn [state] (spin …)) ; the final log potential
:done? (fn [state] …)})
16 {:resampling :stratified})
The target is p(trajectory)·exp(reward): the process's own law tilted by the
reward (Korbak et al. 2022). The process's randomness has no sample site, so
it acts as the proposal and its density cancels; the weights are the
factors. A value estimate twists the intermediate targets — each step adds
log ψ_t − log ψ_{t−1} and the end reward − log ψ_T — which resamples early
without changing the final target (twisted SMC); without one, SMC is
best-of-N weighted by the reward. The model starts at smc/start-site, so
every particle takes every step itself.
A step drawn from another proposal q — a smaller model whose tokens the
large one only scores, a learned proposal — returns
(steer/weighted state log-w) with log-w = log p − log q of the step; the
weight joins the step's factor and the target is still p·exp(reward).
Each particle's trace holds its states ([:steer/state t]) and reward
(:steer/reward); :record ((fn [state]), default the state) chooses
what is kept of a state — a language model's tokens rather than its KV
cache. foerster.amortized.learn/trajectories reads them back with the
particles' normalized weights, and learn/draws resamples them into
unweighted draws from the target: training data for the value estimates and
proposals that make the next search cheaper.
When a static parameter is conjugate to its data, the program can carry its
posterior instead of a sampled value (org.replikativ.foerster.model.conjugate):
score each datum under the posterior predictive, update the posterior in
closed form, and draw the parameter only if it is needed:
(spin (loop [c (conjugate/beta-bernoulli 1 1) remaining (seq ys)]
(if remaining
(let [y (first remaining)]
(observe (conjugate/predictive c) y)
(recur (conjugate/update c y) (next remaining)))
(sample (conjugate/posterior c) :id :p))))
Here ys contains numeric 0/1 Bernoulli observations. The conjugate
helper uses that convention; dist/flip instead produces booleans. The loop
tests sequence exhaustion separately from the observation value.
That is exact marginalization (Rao-Blackwellization; the conjugate case of
delayed sampling, Murray et al. 2018): under SMC the parameter never
degenerates, because it is never resampled, and streaming pays the same per
observation however long the stream. Families: normal-mean (known
observation sd), beta-bernoulli, gamma-poisson (predictive: the negative
binomial) and dirichlet-discrete. A parameter that is not conjugate to its
data needs PMMH or SMC² (below).
A state-space program with parameters θ — named sample sites — is filtered
best by SMC once θ is known; org.replikativ.foerster.inference.smc2 infers θ by
running SMC conditioned on it. With the parameter sites constrained, SMC's
evidence estimate is p(θ)·p̂(y | θ), so two parameter values compare by the
difference of their SMC log-evidences.
;; a Markov chain over θ, each proposal scored by an SMC of 100 particles
(smc2/pmmh (model) 100 2000 {:params #{:sigma} :scale 0.2 :burn 200})
;; online, for a program whose data arrive at stream sites
(smc2/smc2 (model) {:params #{:sigma} :n-theta 50 :n-x 100})
;; → a step like smc/stream's: :push the next observation, :measure over θ
:thetas is
the chain after burn-in, :acceptance the rate over every iteration.:ess-target, they are resampled — a duplicate forks its inner
population's worlds (smc/stream's :fork), so copies evolve
independently — and moved by PMMH steps that run a fresh inner filter on
the data so far. m/log-marginal of its measure is the evidence of the
data so far, observations the program makes before its first stream site
included.θ's prior is the program's own: a θ-particle starts from a simulation of the
program. A proposal that leaves the support — a distribution refusing its
parameters — has density zero and is rejected; one outside the prior's
support is rejected too, after its SMC ran. The default random walk moves
continuous parameters; give discrete ones a :propose of their own: a
foerster.proposal asked about θ (a map {address value}) as the old value
of its site — (proposal/random-walk s) steps every parameter — or a
function (fn [θ]) -> {:theta θ' :log-ratio …}. An asymmetric proposal needs
its log-q, the density of the reverse move. The
caller's :policy (proposals, other constraints) carries into every inner
SMC. Correlated PMMH, which needs SMC driven by explicit noise, is not
provided.
A proposal (org.replikativ.foerster.proposal) says where a latent
site's value comes from when it is not drawn from the site's own law. One
protocol serves every place a value is proposed: the particle methods'
:proposal, MH moves' and resample-move's :propose, and a policy's
:draw.
| Proposal | Draws |
|---|---|
proposal/prior | from the site's own law (the default) |
(proposal/distribution q) | from the distribution q, whatever the old value (an independence proposal) |
(proposal/random-walk s) | a Gaussian step of s from the old value (symmetric), the prior where there is none |
(proposal/only sites p) | p at the given addresses (a set or a predicate), the site's law elsewhere |
(proposal/from-fn propose log-q) | your own: (propose site) → {:value v :log-q lq} or {:value v :symmetric? true} |
;; importance sampling whose :mu is drawn near the data
(infer/importance-sampling (model) 1000
{:proposal (proposal/only #{:mu} (proposal/distribution (dist/normal 7.0 1.0)))})
The weights correct for the proposal (log p − log q), so the target is
unchanged, and an MH move scores its reverse with the proposal's own
log-q: asymmetric proposals are exact. The built-in proposals are right by
construction, and SMC's evidence stays :unbiased; a from-fn proposal is
your obligation, and the evidence is :unverified (foerster.contract).
(kernel-infer model kernel n opts) runs a kernel from
org.replikativ.foerster.moves.kernel (k):
| Kernel | Runs |
|---|---|
(k/single-site-mh-kernel iterations opts) | n independent Markov chains of single-site MH, proposals from the prior; {:only #{address …}} (or a predicate) restricts it to those sites |
(k/random-walk-mh-kernel iterations {:step-size 0.1}) | random-walk MH on continuous sites |
(k/block-gibbs-kernel iterations selector block-kernels classifier opts) | block Gibbs: classifier (fn [address entry]) assigns sites to block ids, selector (k/round-robin-selector, k/random-selector) picks a block per iteration, block-kernels {block-id kernel} moves it (k/prior-block-kernel, k/random-walk-block-kernel on continuous sites, or a custom kernel marked k/PSymmetricBlockKernel); a move that changes the block's membership is rejected |
(k/hmc-kernel iterations {:step-size 0.1 :steps 10}) | HMC on block sites (blocks), and nothing else; :only restricts it to some block sites |
(k/within-gibbs block-kernel) | hmc-kernel or nuts-kernel on the block sites, then one single-site MH move of another latent: a whole model with blocks and other latents (what infer :nuts runs) |
(k/nuts-kernel iterations {:burn … :target-accept 0.8 :max-depth 10 :metric :variance}) | the No-U-Turn sampler on block sites (and nothing else; :only restricts it), Stan's algorithm: multinomial trajectories, the generalized U-turn criterion, divergences, and during :burn dual-averaging step size and windowed diagonal-metric adaptation. :metric :fisher estimates the metric from both the warm-up positions and their gradients, √(Var q / Var ∇log p), the diagonal metric minimizing the Fisher divergence (nutpie; Seyboldt, Carlson & Carpenter 2024). On posteriordb, :fisher gives 2.8× the ESS/s on wells (logistic, 3020 rows). The default :variance gives 1.5× on earnings, a regression with strongly correlated coefficients, and 1.2× on mesquite. Try both on your model; the block's target must be its complete conditional (see below), and its gradient hand-written, compiled by foerster-raster (defdensity), or — an explicit opt-in for small blocks — numeric (block/with-numeric-gradient). Warm-up starts, and restarts after every metric window, from Stan's step size heuristic (init_stepsize), and dual averaging restarts at μ = log 10ε from the step size it found — at the start of warm-up too, a deliberate deviation from Stan, whose first μ comes from its configured step size (default 1). :step-size and :inverse-metric (diagonal, in θ) set where warm-up starts, or with :burn 0 the fixed values NUTS runs at; every retained draw's :sampler-stats report the :step-size and :inverse-metric it ran with, so a warmed-up run's tuning starts a frozen one. Each chain starts at its block sites' :init plus U(−2, 2) per coordinate of θ (Stan's initialization radius, from the chain's seeded stream); :jitter false starts exactly at :init. Jitter stays on with :burn 0, where every draw is kept, so resuming a frozen chain — :init a posterior draw, :burn 0, the exported :step-size and :inverse-metric — needs :jitter false (nuts) |
(k/delayed-acceptance-kernel iterations {address coarse} {:scales 0.1 :steps …}) | two-stage delayed acceptance on block sites (Christen & Fox 2005): a random walk on a cheap coarse log density of θ (a surrogate, a coarser model, a Laplace fit) proposes, one replay of the program decides on its own ratio corrected by the coarse one. The chain targets the program's exact posterior whatever coarse is (a poor one costs acceptance); the block need not be its complete conditional, since the replay scores every observation |
(k/involutive-kernel iterations {:propose :log-q :involution}) | involutive MH moves (foerster.moves.involutive), reversible jump included; chains only (population programs refuse it) |
(k/proposal-mh-kernel iterations proposal {:args […]}) | MH with a proposal program (gfi/mh-proposal), scored forward and backward; chains only |
(k/cycle iterations [k₁ k₂ …]) | each iteration runs k₁'s moves, then k₂'s, … — the composition K₂∘K₁ |
(k/mixture iterations [[w₁ k₁] [w₂ k₂] …]) | each iteration runs one kᵢ, picked with probability ∝ wᵢ |
(k/kernel iterations step opts) | a kernel of your own: step is (fn [trace context]), typically trace/mh-step with your own :select and :propose (extending) |
A chain refuses a kernel that leaves some latent site unmoved
(::transition/uncovered-sites, naming them): a block kernel alone on a
model with discrete latents would keep them at their initial values.
Compose it with a kernel for the other sites, or use within-gibbs.
Kernels that leave the posterior invariant form a monoid under composition,
and their mixtures with state-independent weights are invariant too: cycle
and mixture build a sampler from parts — HMC for the continuous block,
single-site moves for the discrete sites — and nest. Each part runs for its
own iterations.
The Markov-chain kernels run in fresh worlds or, under :world-policy :fork, each chain in a fork of one frozen fork of the caller's world
(incremental reuse only for runs that wrote no world state), and take
output options: :samples :final (default) emits each chain's last
state; :samples :all emits every state after the first :burn moves,
equally weighted. kernel-infer takes :seed, from which each chain's
seed is derived.
Random-walk MH moves a block site on its θ vector (foerster.moves.block-mh):
several Gaussian steps against the block's own density, written back to the
trace by one replay, as NUTS does. The per-coordinate scales adapt during
:burn, steering the acceptance towards 0.234. The block's target must be
its complete conditional; a write-back that changes the log probability of
anything outside the block is refused (::block-mh/incomplete-target).
Other latents take single-site random-walk moves as before.
NUTS and the random-walk block moves follow the block's own density and
write the result back without a correction on the rest of the trace. That
is exact only when the block's density is the complete conditional of its
latents: every factor the latents touch is inside the block, and no other
site depends on them. Latents upstream of the block (other blocks, or
discrete sites moved by MH) may change its inputs; nothing may sit
downstream. A block B whose inputs come from block A's latents makes A's
move change B's score, and A's move is refused
(::nuts/incomplete-target, ::block-mh/incomplete-target; the error's
:changed names the sites whose score moved). Adding p(B | A) to A's
density would not help: the trace scores B's site too, so the factor would
count twice.
In practice NUTS within Gibbs is NUTS on blocks that are leaves of the dependency graph. For a hierarchical model:
k/hmc-kernel) or delayed
acceptance (k/delayed-acceptance-kernel): both accept on the whole
trace's joint density, so their blocks need not be complete conditionals.Lifting the restriction needs the gradient of the downstream factors with respect to the block (composing block densities, or differentiating through the replay); it is not implemented.
org.replikativ.foerster.inference.optimize works on a block site's target directly
(a block/block-dist at its inputs), no sampling:
(optimize/map-estimate (block/block-dist b inputs) {:init θ0})
;; => {:latents {:mu 1.02 :sigma 0.48} :theta […] :log-density … :converged? true}
(def approx (optimize/laplace (block/block-dist b inputs) {:init θ0}))
(:covariance approx) ; of θ, the unconstrained coordinates
(:log-evidence approx) ; log p(θ*) + d/2 log 2π − ½ log det(−H)
((:draw approx)) ; {:mu … :sigma …}, mapped back through the transforms
map-estimate runs L-BFGS on the block's gradient. With constrained
latents it finds the mode in their own coordinates (σ itself), as Stan's
optimize; :jacobian? true gives the mode of the unconstrained θ.
laplace fits a Gaussian in θ at that mode, the Hessian by central
differences of the gradient (as Stan's laplace), and refuses a mode that
is not a strict maximum. It is exact for a Gaussian target and fast; check
it against NUTS when the posterior may be skewed or multimodal in θ.
pathfinder (Zhang, Carpenter, Gelman & Vehtari 2022) needs no mode. It
follows L-BFGS from a few jittered starts. At every iterate, the optimizer's
own inverse-Hessian estimate gives a Gaussian. A Monte Carlo ELBO, computed
with the same standard normals for every iterate, picks the best one. The
paths' Gaussians then propose draws, which are importance-resampled against
the target with Pareto-smoothed weights:
(def pf (optimize/pathfinder (block/block-dist b inputs) {:paths 4 :draws 1000}))
(:draws pf) ; θ vectors, resampled towards the target
(:pareto-k pf) ; above 0.7: the Gaussians miss part of the target
[(:mean pf) (:covariance pf)] ; the best Gaussian, e.g. to start NUTS
; or as a tempering reference
Its gradients are L-BFGS's own, one per optimizer step and line-search
trial; on top, every iterate costs :elbo-draws (30) log densities without
gradient, and the final draws one each. The inverse Hessian is kept as a
diagonal plus a low-rank term and drawn from by Zhang et al.'s Alg. 4, so
an iterate costs O(d·m²) for m ≤ :history pairs rather than a dense d × d
matrix and its Cholesky factor (at d = 100 about 4.7× faster than the dense
version it replaced, with the same approximations). The diagonal is
Zhang et al.'s Alg. 3 estimate recovered along the path by default, as Stan
and BlackJAX do; :diagonal :lbfgs uses L-BFGS's own γ·I. On posteriordb's
regressions its means and standard deviations come within a few percent of
the reference. Where k̂ is high, use it to initialize a sampler, not as the
posterior.
(smc/stream model n opts) (org.replikativ.foerster.inference.smc) is a CPS
operation resolving a step {:measure :push :fork :done? :close} once every
particle waits at its next :stream site; ((:push step) value) scores the
value and runs on (streaming notebook).
((:fork step)) resolves an independent copy of the population — every
waiting particle's world forked — to be pushed and closed on its own.
org.replikativ.foerster.runtime.gfi: simulate, generate, assess, update,
regenerate, mh, close! — CPS operations on traces of a model spin.org.replikativ.foerster.moves.involutive/step: an involutive MCMC move.(counterfactual/counterfactual model {:evidence :interventions :particles}): twin-world counterfactuals of the model spin.See the programmable and counterfactuals notebooks.
org.replikativ.foerster.sandbox is what code running in an interpreter
sandbox — dvergr's SCI sandbox for agents — gets of foerster: namespace maps
the host installs ((sandbox/namespaces host-opts)) and their help text
(sandbox/docs).
| Namespace | Functions |
|---|---|
infer/ | run (a spin of the posterior as plain data), start/status/result/cancel/forget for background jobs, and summary, query, values, weights, log-weights, ess, log-marginal |
model/ | block, site, constrain, numeric-gradient |
opt/ | map-estimate, laplace (with :draws), pathfinder |
diag/ | rhat, ess-bulk, ess-tail, mcse, chains, summary |
lp/ | normal, normal-glm, bernoulli-logit-glm, poisson-log-glm: log densities summed over data vectors on the host, about 5× faster than the same sum interpreted |
The methods available there are :importance, :smc, :mh, :rmh and
:nuts; the options are a data-only allowlist (no :policy, :kernel,
:executor or :world-opts). Results are plain data with metadata
stripped, and errors cross as their message and :type only. The host
decides the default :world-policy, caps the work of a run (:max-work:
particles, or iterations × chains) and of an optimizer
(:max-optimizer-work), and gives a session :seed from which every run
without its own seed derives one. A posterior whose particles all have zero
weight is refused (::sandbox/zero-likelihood).
Particle MCMC is not offered there yet: particle Gibbs resets the weights of each sweep, so a model with zero likelihood everywhere comes back as a uniform posterior rather than an error, and a sandboxed caller could not tell. Tempering, parallel tempering and enumeration are left out because their sizes do not bound their cost.
Complete keyed regions use incremental re-execution through the same public MH
and kernel paths. (k/finite-gibbs-kernel iterations address opts) composes
fixed-address finite-support Gibbs with other Markov descriptors. It requires
fixed site structure and unchanged auxiliary choices. See the
incremental inference guide for prefix fallback,
unsupported dynamic structures, resource ownership and complexity limits.
Can you improve this documentation?Edit on GitHub
cljdoc builds & hosts documentation for Clojure/Script libraries
| Ctrl+k | Jump to recent docs |
| ← | Move to previous article |
| → | Move to next article |
| Ctrl+/ | Jump to the search field |