Contributing¶

Build from source¶

The package is a compiled Rust extension built with maturin, with uv managing the environment and lockfile.

make install           # uvx maturin develop
make install-release   # uvx maturin develop --release (use for benchmarking)

Day-to-day commands¶

Prefer the Makefile targets so flags stay consistent with CI:

make test       # POLARS_MAX_THREADS=4 uv run --group testing pytest tests
make typing     # pyrefly + pyright + mypy (all three, as in CI)
make lint       # prek hooks (ruff, rumdl, ryl) + cargo fmt (nightly) + clippy
make benchmark  # polars_stats vs scipy comparison report (tools/benchmarks/)

make test caps POLARS_MAX_THREADS=4 on purpose: it forces multi-thread, multi-chunk execution so the chunk- and thread-invariance of sample actually gets exercised.

It does not pin a query engine, so it runs whichever one your environment resolves (POLARS_ENGINE_AFFINITY, else the polars default). CI runs both, because the two chunk a plugin's inputs differently: in-memory calls it once over the whole column, streaming once per morsel. A chunk-boundary or input-length bug can therefore pass under one and fail under the other. Pin an engine with:

POLARS_ENGINE_AFFINITY=in-memory uv run --group testing pytest tests
POLARS_ENGINE_AFFINITY=streaming uv run --group testing pytest tests

To preview these docs locally:

uv run --group docs zensical serve

Repository layout¶

polars-stats/
├── Cargo.toml
├── pyproject.toml
├── rust-toolchain.toml
├── src/
│   ├── lib.rs                # pymodule entry + global allocator
│   ├── rng.rs                # per-row RNG + the sampler drivers every plugin shells over
│   └── distributions/        # one Rust file per distribution
├── polars_stats/
│   ├── __init__.py           # public exports
│   └── distributions/
│       ├── _base.py          # ABCs + coercion helpers
│       └── _<name>.py        # one Python class per distribution
├── tests/
│   ├── _registry.py          # one DistSpec row per distribution, plus the shared contract frame
│   ├── conftest.py           # the fixtures that deliver it
│   ├── distributions/        # flat, one file per behavioural contract, parametrised over the rows
│   ├── property/             # hypothesis-based invariant tests
│   └── scipy_parity/         # scipy reference-oracle tests
├── tools/                    # hand-run tooling, outside the package and the test suite
│   ├── accuracy/audit.py     # the mpmath tail-accuracy audit (`make audit`)
│   └── benchmarks/           # the polars_stats vs scipy harness (`make benchmark`)
└── docs/                     # this documentation, one directory per Diataxis quadrant

Stack¶

Rust runtime¶

Deliberately excluded: rand_chacha (replaced by rand_pcg; per-row ChaCha20 construction made sampling markedly slower, see Design notes), ndarray (Polars is Arrow-native), scirs2-stats (pre-1.0). There is no direct rayon dependency either: the multi-draw fill in rng.rs parallelises through the re-export in polars-core, on the same POOL Polars uses.

Python runtime¶

Only polars>=1.15. No other runtime dependencies.

Dev and CI tooling¶

• Dev dependencies are grouped in pyproject.toml (testing, tools, typing, docs) and installed with uv sync --group ....
• Tests run under pytest with scipy + numpy as the parity oracle (tests/scipy_parity/) and hypothesis for property tests; the tools/benchmarks/ comparison report measures wall-clock time and peak memory against scipy.stats.
• Python is checked by ruff (lint + format) and three type checkers in CI (mypy, pyright, pyrefly); Rust by cargo fmt (nightly) and cargo clippy --all-features --all-targets -- -D warnings. Prose and config are linted by rumdl (Markdown), ryl (YAML), codespell, typos, and blacken-docs, wired through .pre-commit-config.yaml and run by prek both locally and in CI, so the hook versions pinned there are the single source of truth. Do not bypass the hooks with --no-verify unless asked.

Adding a distribution¶

Copy an existing pair, modify, repeat. Use Uniform (_uniform.py / uniform.rs) as the canonical continuous example and Bernoulli as the discrete one.

Each distribution opened up for contribution gets a GitHub issue carrying its exact spec and checklist: the parameterisation, the scipy equivalent, every closed form, the contract it settles, and the traps found the last time someone tried. That issue is canonical if it conflicts with this section, so please say so on the issue when it does. Ones labelled good first issue need no new Rust math; the rest need a special-function port. Resolve the choices the issue flags (a validator arity, a scipy reparameterisation, an out-of-regime convention) before writing code rather than halfway through, write the scipy-parity test first, and keep a pull request to one distribution.

1. Rust. Add src/distributions/<name>.rs, register it in src/distributions/mod.rs, and implement a #[polars_expr] only for methods that need it. Which ones those are, method by method, is tabulated in Design notes. Whatever you add, every plugin function owes three things:

   • align_inputs(inputs)? first, before any cast, for any plugin with more than one input. Polars broadcasts nothing into a plugin and try_*_elementwise truncates to its shortest input, so skipping it silently drops every row past the first. tests/distributions/broadcast_test.py catches a missed one. The value-keyed drivers in mod.rs skip it when every parameter is length 1: they check the constant once and never expand it. Route through them rather than writing that branch yourself.
   • is_elementwise=True, which the shared register_plugin shim fixes for you. An aggregating plugin would break over and group_by, so it is a guard rather than a default.
   • Null in, null out. The try_*_elementwise drivers give you that; a raw .into_iter() over chunks does not, which is the one way to lose the contract while writing ordinary-looking Rust.
   • Always the samplers. Every one is a one-line shell over a driver in src/rng.rs; none resolves a seed, coerces an input or writes a row loop itself, which is what keeps seeding, null propagation (a null in any input nulls the row) and the invalid-parameter error contract in one place. The per-row <name>_sample / <name>_samples (multi-draw, backing samples) take the parameter columns plus a row index as the last input and call sample_per_row_binary / samples_per_row_binary (one parameter, passing its ParamDomain) or sample_per_row_ternary / samples_per_row_ternary (two, passing each parameter's coercer and the check_params pass), plus a build that constructs the row's draw state. The constant-parameter fast paths <name>_sample_scalar / <name>_samples_scalar build once from SampleScalarKwargs<<Name>Params> / SamplesScalarKwargs<<Name>Params> and hand sample_by_index / samples_by_index the row index and a draw closure. Every driver takes its output dtype from what draw returns, so no call site names a polars type. Give the distribution one named fn draw and call it from all four shells; that shared call, not a test, is what keeps them byte-identical.
   • When a method needs a special function (erf, log-gamma, regularized incomplete beta/gamma, ...) or has no elementary closed form: bind it in statrs (pdf / pmf, cdf, ppf, ln_pdf / ln_pmf, native sf, native median). Write one named *_value body per method, fn cdf_value(dist: &Dist, v: f64) -> Option<f64>, and call it from two one-line shells: the per-row <name>_<method> over the distribution's own value_keyed helper (itself one line over value_keyed_ternary in src/distributions/mod.rs, or value_keyed_binary for one parameter, passing the coercers, check_params and build_dist), and the constant-parameter <name>_<method>_scalar over kwargs.value_keyed(&inputs[0], <method>_value), an inherent method on the distribution's <Name>Params struct that builds once through Self::build and maps the body through value_keyed_scalar. The driver's constant branch (every parameter length 1) calls the same value_keyed_scalar with the same body, so a Python scalar and a length-1 literal or aggregate run one instantiation and agree bit for bit. Putting the builder on the kwargs struct is what makes the hoist structural: a shell cannot build the distribution itself, so it cannot rebuild per row.

   • When a method is an elementary closed form (no special function: a Normal's 0.5*log(2*pi*e*sigma^2), a Bernoulli's 1 - p at 0, mean = n*p, ...), the test is whether every row reaches the validator, not whether Polars can express the arithmetic:

     • The validated parameter is read unconditionally (Bernoulli.mean is the validator's own output), or it is named in a when(...) condition (_moment's when(validator.is_not_null()) gate): leave it in Python as a Polars expression. A Rust binding for arithmetic Polars does natively is dead FFI surface. In practice this is the moments: mean, variance, std, median, entropy.
     • The validated parameter is read only inside a when / then / otherwise arm, which is what every method that branches on the evaluation value looks like: compute it in Rust. From polars 1.44 an arm is masked to null on the rows it does not select, so the validator never sees an invalid parameter on those rows and the method silently returns a value instead of raising (pola-rs/polars#29005), and no Polars expression can spell a guard that survives the optimiser. In practice this is the value-keyed set: pdf / pmf, log_pdf / log_pmf, cdf, log_cdf, sf, log_sf, ppf, isf. DiscreteUniform, Bernoulli, Exponential, Geometric and Uniform are the worked examples: the base routes every value-keyed method to its <name>_<method> Rust body, and only the moments stay in _<name>.py.
     • Split such a body into a derive that turns the parameters into their branch answers and a select that picks one by the evaluation value. Where the branch answers are fixed once the parameters are, the table is non-generic and derive is an associated constructor (bernoulli.rs's Mass::pmf / Mass::at, uniform.rs's Density::pdf); where a branch still depends on the evaluation point, the table carries an Arm type parameter and derive is a free function (exponential.rs's, geometric.rs's, pareto.rs's and weibull.rs's derive_cdf over the shared Sides<Arm> in mod.rs, whose floor the four set to 0, 1, the scale and 0; uniform.rs's derive_cdf / Regions::at). An inverse needs no table at all: derive returns the arm and select is the shared on_unit_interval. derive runs once per call on the constant path and once per row when a parameter is a column; a body that recomputed the parameter-only terms per row regressed the constant-parameter path by up to 195% at 10M rows.
     • The per-row shell is one line over value_keyed_derived_binary(inputs, &DOMAIN, derive, select) (one parameter) or value_keyed_derived_ternary(inputs, check_params, derive, select) (two); these wrap value_keyed_binary / value_keyed_ternary with derive as an infallible build, so there is one row loop per arity in the crate. The _scalar twin is kwargs.value_keyed(&inputs[0], derive, select), an inherent method that checks the constant and calls value_keyed_scalar(value, &derive(p), select), the same call the driver's constant branch makes.
       • State each parameter's domain once as a ParamDomain constant (ParamDomain::finite("mu"), ParamDomain::positive("sigma"), ParamDomain::probability("p"), or a literal for any other rule) and a joint constraint as a PairDomain. Every column-parameter plugin runs check_column / check_columns over its parameter columns before its row loop, and every <Name>Params::build runs check once, so an invalid parameter raises the same ComputeError in both regimes whatever else its row holds. Keep a build_dist(...) -> PolarsResult<Dist> around the statrs constructor for the row loops; behind the pass it cannot fail.
     • Python. Add polars_stats/distributions/_<name>.py, subclassing ContinuousDistribution or DiscreteDistribution. In __init__, coerce each parameter with coerce_param / coerce_n (types only, never validate values at construction, let invalid values raise in Rust) and store the fast-path bundle self._scalar_kwargs = scalar_kwargs(...). The base owns all routing (sample, samples, _value_plugin, _param_plugin, _validated, _moment); a subclass declares only what is distribution-specific:

   • _distribution_name: ClassVar[DistributionName] (e.g. "normal"), the prefix of every plugin the base calls: the samplers (<name>_sample / _sample_scalar / _samples / _samples_scalar), the value-keyed hooks (<name>_pdf, <name>_ln_cdf, ...) and the validator. Constant parameters route to the _scalar twins off _scalar_kwargs. Add the name to DistributionName in polars_stats/_typing.py.

   • _param_exprs, the coerced parameters as a tuple in plugin-input order; the first sets the output's root name.
   • Nothing for the value-keyed methods. The base hooks (_pdf / _pmf, _log_*, _cdf, _sf, _ppf, _isf) each call self._value_plugin("<method>", value), which resolves <name>_<method> (ln_ for the log forms). Override a hook only where the Rust body does not exist (Beta and Binomial compose log_cdf / log_sf as log(cdf) / log(sf)).
   • _validated_params, the validating plugin call the moments gate on. In Rust, a plugin returning a reused quantity that raises on invalid parameters and nulls on a null one (uniform_range returns max - min, normal_sigma the validated sigma); its suffix must be in ParamFunction. Cast each input, run the domain pass, then return the quantity: validated_param(inputs, &DOMAIN) returns the checked column itself for one parameter (bernoulli_proba); validated_pair(inputs, coerce_a, check_params) returns the second parameter where both are present (normal_sigma, beta_params); param_keyed(inputs, coerce_a, coerce_b, check_params, |a, b| ...) returns any other per-row quantity, null where either parameter is null (uniform_range, beta_entropy). It takes no output name: polars resolves an expression's output name from its first input. In Python, _validated_params is self._validated("<suffix>", <the same quantity as a polars expr>), which runs the plugin per row for column parameters and once, on length-1 literals, for constants; or self._param_plugin("<suffix>") when the plugin's own output is the answer on both routings (DiscreteUniform.support_size, whose Float64 count must not be recomputed in polars).
   • The closed-form moments, as self._moment(<formula over the raw parameters>): one validator mention however many times the formula names a parameter. A moment whose formula is the validated quantity reads _validated_params directly (Bernoulli.mean, Exponential's five, Uniform.range). A moment with no closed form is its own parameter-keyed plugin, self._param_plugin("entropy"), validating inside.

   Never override the public pdf / cdf / ... methods, nor the base-owned sample / samples / _value_plugin / _param_plugin / _validated / _moment. Export the class from polars_stats/__init__.py. 3. Tests. A new distribution touches its own files and the shared registry.

   The registry row is the one thing that cannot be forgotten quietly

   Missing it used to be the only failure mode here with no signal at all: the suite did not know your distribution existed, so it could not report that it had skipped it. tests/distributions/registry_test.py now cross-checks polars_stats.__all__ against the rows and fails when they disagree.

   All of:

   • tests/_registry.py: one DistSpec row, and it is the only registry there is. The row carries name, cls, continuous, parameters (each with its constructor keyword, its column in the shared contract frame, and its type), a param_strategy, an example, eval_range, bounds, on_support_point, sample_dtype, the finite out-of-domain invalid table, and integration_bounds (continuous) or support (discrete). Add a column per parameter to CONTRACT_FRAME in the same file; no two parameters of one distribution may share one. Add the name to DistributionName in polars_stats/_typing.py, which is what types the row and the four sparse tables below. Every shared contract reads the row, and tests/distributions/registry_test.py fails if a distribution is exported without one, so this is the one place the old silent-skip cannot happen.
   • Possibly a row in one of the four sparse tables, which a DistSpec field cannot hold because most distributions have no entry. support_test.py::_DENSITY_AT_ENDPOINT (the density at a finite support endpoint, whose value is the distribution's own rather than a saturated constant) is required when the support has a finite endpoint, and fails red when missing. _registry.py::DEGENERATE (parameterisations where the mass collapses onto one point), _registry.py::ULP_TOLERANT_MOMENTS (moments that are not bit-exact across the two parameter routings) and _registry.py::UNDEFINED_MOMENTS (moments the distribution does not have, pinned to null on every valid row) are all opt-in; add to ULP_TOLERANT_MOMENTS only from a failing bit-exactness assertion, never to quiet one.
   • Any bespoke facts, appended to the file that owns their subject in tests/distributions/: an algebraic reduction to another distribution goes in identities_test.py, a numerical-regime fact with no scipy oracle in precision_test.py. Everything a shared contract already states needs no new test at all; see tests/README.md for the file-by-file map.
   • tests/scipy_parity/<name>_test.py: one Case per method against scipy.stats.<name> through the shared _harness.py (default absolute tolerance 1e-12, relaxed per Case to 1e-9 / 1e-6 for erf-based or binary-search-ppf methods).
   • Update the umbrella issue when the change merges.

Numerical stability¶

Every method must be accurate in the regime it exists to serve. log_sf exists for the deep tail, so a log_sf that returns -inf there does not work, even though every test passes. Composing log_cdf as cdf().log() (as Beta and Binomial still do) is a stopgap, not an implementation: shipping one is a decision to justify. Never write isf as ppf(1 - q), because the loss would happen before your code runs: 1 - q resolves to 1.1e-16 absolute, so the tail mass is already quantised to 1.1e-16 / q relative and no inverse can recover it. Solve against q itself, via a symmetry, a closed form, or entering a two-sided solve from the other end.

A composed method inherits the weakest part's range, and the composition is often wider than the part. std() defaults to variance().sqrt(), and a variance that legitimately overflows can hide a standard deviation that does not. Ask what the domain of the composed quantity is, not what the domain of its parts is.

Reassociate before reaching for log space or a branch. Computing a density as exp(log_pdf) fixes the subnormal range but costs one to two orders of magnitude everywhere else, because exp(log(rate)) does not round-trip to rate, so it needs a threshold and when/then/otherwise evaluates both sides. (rate * exp(t / 2)) * exp(t / 2) is the same product with no branch and no threshold constant. Rearrange so the single unavoidable rounding happens last.

The test suite cannot catch this for you. Three structural blind spots, each of which has produced a real defect:

• parity grids are finite and curated, so they never probe the extreme regime;
• scipy is not a valid oracle in the tails: its own logcdf / logsf are naive for the incomplete-beta and incomplete-gamma families, so parity passes while both libraries return -inf;
• property tests assert shape, not accuracy: monotonicity and mass-integrates-to-one are satisfied by a value that is relatively wrong by 1e-3.

Choose the oracle deliberately, and say which one you used and why. scipy where it is finite and known accurate; mpmath at high precision, or an exact closed-form identity, beyond that. Never assert against a scipy value that is itself saturated. An oracle for a discrete inverse needs exact rational arithmetic, since a rounding there becomes a jump between support points rather than a small error.

Two acceptable outcomes, and no third. Fix the algorithm, or document the caveat in Numerical accuracy, quantified with a regime and a magnitude ("relative error ~1.1e-16 / q", not "may be inaccurate in the tails"). That page is the single home for accuracy caveats: no per-class docstring blocks, no restating the same limit in several places. A runtime warning is never the fix: it cannot fire per-row from inside the engine, and a Python-side scalar-only warning would break the scalar/column symmetry the architecture guarantees.

Claim a tolerance and justify it. 1e-12 for elementary closed forms, 1e-10 for special-function methods (2e-10 through erfc, which is what statrs holds), 1e-9 for log-scale, 1e-8 for discrete log-mass and support-sum entropy, 1e-6 or integer-valued for a discrete binary-search ppf. A relaxed tolerance needs a one-line reason in the Case, not a shrug.

statrs' inverse_cdf is not a safe default. For the discrete families it is a binary search, so its parity tolerance is the search's convergence rather than a formula's: state the 1e-6 loosening in the Case. Where it is not a binary search it has been wrong in several ways, including relatively wrong by 6e-3 for Gamma in the low-quantile tail, and panicking, hanging and saturating for Beta. A bounded solve with every Newton proposal clamped into a bisection bracket is the pattern to copy.

Run make audit for a new distribution, and add its oracle to the registry in tools/accuracy/audit.py. A distribution absent from the audit is unaudited, exactly as one absent from tests/_registry.py is untested; make audit now fails outright on a distribution with no oracle. Never bound the sweep by what the implementation is known to be bad at: that is the defect's own shape used as a bound on the instrument. Sweep extreme parameters too, not only extreme inputs.

Ask what class a finding belongs to before closing it, and treat "this method is exempt" as a hypothesis to probe rather than an argument to accept. Both rules were bought the same way: a pdf-in-log-space defect fixed once and found again elsewhere, and three distributions cleared from the isf fix by reasoning that each fell to the first probe aimed at it.

Where the recipes live. exponential.rs's derive_ln_sf (an exact closed form), its derive_cdf and LogNormal.variance (the sinh identity standing in for the expm1 Polars does not expose), its derive_ln_cdf and uniform.rs's derive_ln_cdf (log1p on the near-certain side), geometric.rs's smallest_support_point (the one-ulp tie rule of a discrete inverse, decided in the log domain both sides entered through), normal.rs's ln_erfc (a special function ported to log space, the pattern for the hard cases), the isf_value bodies in normal.rs (a symmetry) and lognormal.rs (composing one), Sides::read_at in mod.rs with pareto.rs's and weibull.rs's log_ratio (one distribution read as another on a transformed variate, the transform formed from an exact difference so the support edge keeps its digits), and weibull.rs's second_moment_terms (a difference of two special-function values replaced by its own series where the two cancel, with the crossover chosen where the direct form still has the digits the series keeps).

Conventions¶

Prose, in PRs and docs: no em dashes or double hyphens (use commas, colons, or parentheses); lead with the answer; state uncertainty explicitly; end non-trivial answers with a short "Blind spots" section.

Code: KISS and YAGNI; production-grade type hints and explicit error handling; default to no comments, add one only when the why is non-obvious. Do not refactor speculatively, do not introduce backwards-compatibility shims unless asked, and do not add the DistKind dispatch macro (rejected for v1, see Design notes).

Tests assert on pl.Series / pl.DataFrame via polars.testing, not Python lists. For random output, assert the null mask, not values. Genuinely scalar results are fine to read with .item(...) and compare via pytest.approx. scipy parity stays on numpy assert_allclose.

Git: one change per distribution; commit messages mirror the existing history (feat:, fix:, refactor:); add a co-authored-by trailer for AI-generated commits. Do not use destructive git operations without explicit approval.