stanli

The Stan Language Interpreter. Compile and sample Stan models with no C++ toolchain on the machine.

pip install stanli

That is the whole install. No compiler, no make, no CmdStan checkout, no multi-minute first-run build. One wheel, one shared library, under eight megabytes. In the current Eight Schools benchmark, a first complete 1,000-warmup, 1,000-draw run takes 0.03 s in stanli versus 3.4 s to build and run the model with CmdStan - roughly 100x faster from source to CSV.

import stanli

model = stanli.Model(stan_file="eight_schools.stan", data="data.json")
fit = model.sample(seed=1, chains=4, warmup=1000, samples=1000)

fit["mu"].mean()        # every draw of a column, chains concatenated
fit.draws("mu")         # (chains, draws), for a trace plot

Sampling reports CmdStan-shaped progress every 100 transitions by default, followed by per-chain warm-up, sampling, and total times:

Chain [1] Iteration:    1 / 2000 [  0%]  (Warmup)
Chain [1] Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain [1] Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain [1] Iteration: 2000 / 2000 [100%]  (Sampling)

Chain [1] Elapsed Time: 0.821 seconds (Warm-up)
                        0.169 seconds (Sampling)
                        0.990 seconds (Total)

Set refresh=0 for a completely quiet run, or another positive integer to change the update interval. Progress is written through Python's sys.stdout, so notebook output and contextlib.redirect_stdout() work normally. Reporting only observes completed transitions: changing refresh does not change draws, sampler statistics, generated-quantity RNG streams, or reproducibility. If any post-warmup transition diverges or saturates the maximum treedepth, the final output also reports the aggregate count. Those counts cover all transitions, including transitions omitted by thinning.

Chains and convergence

Four chains by default, run in parallel, because R-hat needs more than one and a single-chain run cannot be checked for convergence at all. Eight schools does all four in about 70 ms. Threading changes nothing about the answer: each chain owns its executor and its RNG stream, so the draws come out byte-identical to a sequential run.

print(fit.summary())

name                Mean       MCSE     StdDev         5%        50%        95%   ESS_bulk   ESS_tail      R_hat
mu                4.4600     0.0532     3.1705    -0.7414     4.5519     9.5384       3586       2847      1.000
tau               3.4752     0.0635     3.1612     0.2192     2.6680     9.6313       2160       1874      1.001

R-hat is rank-normalized split-R-hat and ESS is the bulk/tail pair (Vehtari et al. 2021), computed by stan's own estimators, so the numbers agree with stansummary rather than approximating it.

print(fit.diagnose())

No divergent transitions.
No transitions saturated the maximum treedepth of 10.
E-BFMI is above 0.3 in every chain.
R-hat is below 1.01 for every parameter (worst 1.002, theta.6).
Bulk ESS is at least 100 per chain for every parameter (worst 2160, tau).
Tail ESS is at least 100 per chain for every parameter (worst 1874, tau).
No problems detected.

Those are the checks a Bayesian workflow actually turns on, including E-BFMI, the one that catches a badly explored heavy tail, which R-hat and ESS are both blind to. The pieces are reachable individually too: fit.divergences, fit.max_treedepth_hits, fit.stepsize and fit.ebfmi() are per-chain arrays, and fit.to_arviz() hands off an InferenceData with the sampler stats attached.

The mode, and where to start

r = model.optimize(seed=1)
r["mu"], r.lp          # every CSV column at the mode, and the lp there
r.unconstrained        # the point on the sampler's scale

fit = model.sample(inits=r.unconstrained)   # start the chains there

L-BFGS, stan's own, the one behind CmdStan's optimize. It returns the posterior mode. CmdStan's optimize defaults to jacobian=0, the penalized maximum likelihood, and stanli cannot offer that: the change-of-variables Jacobian is folded into the graph when the model is lowered. jacobian=False raises rather than quietly handing back the other quantity.

How it works

Every Stan model is a composition of a fixed vocabulary of operations: densities, constraint transforms, linear algebra, elementwise math. stanli ships those precompiled and turns each model into data, a static graph of ops over flat preallocated buffers, instead of generating and compiling C++ per model. The graph doubles as the autodiff tape, so a reverse sweep is a backwards loop over an array, and steady-state gradient evaluation allocates nothing.

Two things are not reimplemented, which is what makes the results trustworthy: the compiler is the real stanc3 plus the stanli OCaml pipeline, embedded in the runtime on macOS and Linux and packaged as an executable on Windows, and the math is unmodified stan-math, the same code CmdStan runs.

Correctness

Nothing here ships on "looks close".

118 of 120 posteriordb models are differentially verified against CmdStan: same model, same data, same evaluation point, comparing the log density and every single gradient component. 41 agree bitwise. The worst deviation across the entire corpus is 2.6e-12 relative.

The two exceptions are documented rather than hidden. sir's ODE solution dips about 1e-9 below a declared lower bound at the shared evaluation point, where CmdStan rejects it too; kronecker_gp matches on the log density and 436 of 438 gradients, differing on the two that flow through eigenvectors of a nearly degenerate covariance matrix.

Full per-model accuracy table: docs/corpus-status.md

Performance

Per-gradient latency against CmdStan, same models, same evaluation point, both sides -O3 with FP contraction pinned off:

The wins come from op granularity. CmdStan's var tape allocates, walks, and frees one node per scalar operation per leapfrog step; stanli pays a fixed cost per op, and a vectorized statement over N elements amortizes that to nothing. Across the whole posteriordb corpus the median is 2.91x and

Repeated independent work produces the largest wins. Dense kernels and serial recurrences land closer to parity because there is less work to fuse, though the measured HMM, ARMA, and GARCH gradients are now all at least as fast as CmdStan. The only gradient losses are the three ODE models at 0.87-0.90x; their first complete runs still win once CmdStan's model build is included.

Method and full table: docs/benchmarks.md

API

The surface is small on purpose.

import stanli

# A path to a .stan file, or the model source directly.
model = stanli.Model(stan_file="model.stan", data="data.json")
model = stanli.Model(stan_code=src, data={"J": 8, "y": y, "sigma": sigma})

model.n_unconstrained               # length of the unconstrained vector
model.constrained_names             # ['mu', 'tau', 'theta.1', ...]

lp, grad = model.log_prob_grad(q)   # sampling log density and its gradient

fit = model.sample(seed=1, warmup=1000, samples=1000, delta=0.8,
                   refresh=100)
fit["theta.1"]                      # ndarray, chains concatenated

data accepts a path to a JSON file or a dict of Python scalars, lists, and numpy arrays. sample returns every column CmdStan's CSV would carry (constrained parameters, transformed parameters, generated quantities, with RNG draws streamed per chain), named the way CmdStan names them, so theta declared as vector[8] arrives as theta.1 through theta.8. Sampler columns (lp__, divergent__, ...) are reachable by name too.

Platforms

Wheels for macOS (arm64 and x86_64), Linux (x86_64 and aarch64, manylinux_2_28) and Windows (x86_64). The Windows wheel is built under mingw-w64, because stan-math does not build under MSVC (the same reason RStan ships through RTools). Windows wheels built from this revision run stanli-compile.exe as a short-lived subprocess and keep pristine stanc.exe beside it for one rollback cycle; there is no OCaml compiler DLL. Python prefers stanli-compile.exe whenever it is present and selects stanc.exe only when it is absent. Launch errors, compiler errors, and empty output are reported; an invalid portable result is rejected when decoded. None causes an automatic retry through stock stanc. The public API works the same way as the embedded compiler path.

The installed library is 29.7 MB: over half of it is the density kernels, about a quarter the embedded stanc3, and the interpreter and NUTS together are about 410 KB. That is the trade this design makes: ship the compiler and every kernel once, so nothing is ever built on the user's machine.

Limits

Stated plainly:

• The sampler is Stan's own NUTS with diagonal-metric adaptation, and optimize() is Stan's L-BFGS. No variational inference or Pathfinder yet.
• inits are on the unconstrained scale. Constrained inits would need the inverse parameter transforms, which do not exist here yet.
• optimize(jacobian=False) (CmdStan's default penalized maximum likelihood) raises; see above.

What is here is verified against CmdStan model by model, and every number on this page is reproducible from the repository.

• Source, issues, and roadmap: github.com/seantalts/stanli
• License: BSD-3-Clause, matching Stan's own.