Statistical inference,
the next step.
Statistical inference at Rust speed. From A/B tests to particle physics, one engine for frequentist and Bayesian methods. SIMD, CUDA, Metal — from data to results.
No C++ dependencies.
No install headaches.
Native readers for ROOT TTree, HS3, Arrow IPC, and Parquet — all with zero external dependencies. mmap I/O, rayon-parallel decompression. No conda, no shared-library pain.
One command.
Full pipeline.
ROOT ntuples → histograms → HistFactory workspace → fit → CLs limits. Inspired by TRExFitter's workflow, reimplemented in Rust for native speed. No intermediate files.
nextstat fit --input workspace.jsonpyhf-compatible.
Validated to match.
pyhf set the standard for pure-Python HistFactory. NextStat reads the same JSON workspaces and validates against pyhf with a 7-tier tolerance contract (Δ < 1e-8). Barlow-Beeston lite for bin-by-bin stat errors. No rewriting required.
50+ statistical methods.
One engine.
Reproducibility
built in.
Deterministic mode with --threads 1 and --seed for bit-exact reproducibility. Immutable run bundles. Every toy, every fit, every scan — traceable.
Built for physics.
Ready for everything.
NextStat started in high-energy physics. The same engine powers survival analysis, econometrics, and time series.
High-energy physics
- HistFactory workspaces (pyhf JSON)
- Profile likelihoods & CLs limits
- Ntuple → workspace pipeline
- TRExFitter config compatibility
Differentiable & RL
- PyTorch: zero-copy CUDA gradients through NLL
- Gymnasium RL environment for signal optimization
- Profiled CLs as differentiable loss
- Complementary to neos / MODE approaches
Applied statistics
- GLM / regression / hierarchical Bayes
- Survival (KM, log-rank, Cox PH, parametric)
- Time series (Kalman / EM)
- Econometrics (DiD, IV/2SLS, AIPW)
Compute backends
- CUDA kernels for NVIDIA GPUs (f64)
- Metal shaders for Apple Silicon (f32)
- SIMD + Apple Accelerate (vDSP/vForce)
- WASM playground in the browser
One engine. Five verticals.
Kaplan-Meier, log-rank, Cox PH, GLM, and bootstrap share the same Rust core. A single pip install nextstat unlocks survival analysis across pharma, insurance, epidemiology, reliability, and subscription analytics.
Pharma & Biotech
Clinical trial survival endpoints, KM curves, treatment arm comparison, Cox PH hazard ratios.
Insurance & Actuarial
Loss reserving, Gamma/Tweedie GLM for pricing, extreme value (GEV/GPD) for reinsurance.
Epidemiology
Population survival, cohort studies, meta-analysis with heterogeneity diagnostics.
Reliability Engineering
Component failure analysis, Weibull fits, time-to-event with right censoring.
Subscription & Churn
Cohort retention curves, churn risk models, causal uplift estimation (AIPW).
import nextstat
# Kaplan-Meier survival estimate
km = nextstat.kaplan_meier(times, events, conf_level=0.95)
print(f"Median survival: {km['median']:.1f}")
# Log-rank test: compare two treatment arms
lr = nextstat.log_rank_test(times, events, groups)
print(f"p = {lr['p_value']:.4f} (chi² = {lr['chi_squared']:.2f})")
# Cox PH hazard ratios
fit = nextstat.survival.cox_ph.fit(times, events, x, robust=True)
print(fit.hazard_ratio_confint())Your analysis is an RL environment.
Train an RL agent to optimize signal histograms for maximum discovery significance. NextStat wraps the full HistFactory likelihood as a standard Gymnasium environment — bringing reinforcement learning to statistical analysis workflows.
from nextstat.gym import make_histfactory_env
env = make_histfactory_env(
workspace_json,
channel="SR",
sample="signal",
reward_metric="z0", # discovery significance
action_mode="logmul",
max_steps=64,
)
obs, info = env.reset(seed=42)
for _ in range(64):
action = agent.predict(obs)
obs, reward, done, _, info = env.step(action)
if done:
breakGymnasium-compatible HistFactory environment
Standard Gymnasium API — works with Stable Baselines3, CleanRL, RLlib, or any Gym-compatible agent out of the box.
Real physics as reward
NLL, discovery significance (q₀/Z₀), exclusion power (qμ/Zμ) — not toy proxies.
Gradient-free + gradient-based
Use PPO/SAC via Gymnasium for black-box search, or switch to nextstat.torch for end-to-end differentiable pipelines. Complementary to neos and MODE's gradient-based approaches.
One line to switch reward
reward_metric="nll" for fast loops, "z0" for profiled discovery, "zmu" for exclusion — same env.
Train your NN with physics as the loss.
Instead of training a classifier with cross-entropy and then running a statistical test, NextStat lets you differentiate through the test itself. Your neural network directly optimises discovery significance — with all systematics profiled out automatically.
from nextstat.torch import SignificanceLoss, SoftHistogram
model = nextstat.from_pyhf(workspace_json)
loss_fn = SignificanceLoss(model, "signal")
soft_hist = SoftHistogram(bin_edges=torch.linspace(0, 1, 11))
for batch_x, batch_w in dataloader:
scores = classifier(batch_x) # NN → scores
histogram = soft_hist(scores, batch_w) # → soft bins
loss = loss_fn(histogram.double().cuda()) # → -Z₀
loss.backward() # gradients → NN
optimizer.step()Differentiate through the test statistic
Your loss is the profiled discovery significance Z₀ — not a proxy. Systematic uncertainties are marginalised in the forward pass.
SoftHistogram bridges NNs and binned analysis
KDE or sigmoid differentiable binning. NN scores become soft histograms with full gradient flow.
DLPack / Array API interop
as_tensor() bridges JAX, CuPy, Arrow, and NumPy arrays to PyTorch — zero-copy on GPU via DLPack.
MLOps built in
metrics_dict() for W&B / MLflow. rank_impact() for systematic feature importance. signal_jacobian() for fast pruning.
Every event counts.
Unbinned analysis preserves full information from every event — no loss from binning. 13 parametric PDFs, extended unbinned likelihood, declarative YAML specification, and the same optimizer stack as HistFactory.
YAML Specification
# unbinned_spec.yaml — parametric unbinned analysis
observables:
- name: mass
range: [5.0, 6.0]
unit: GeV
processes:
- name: signal
pdf:
type: crystal_ball
params: { mu: mu_sig, sigma: 0.02, alpha: 1.5, n: 2.0 }
yield: { expr: "mu * 500.0" }
- name: background
pdf:
type: chebyshev
params: { c1: c1_bkg, c2: c2_bkg }
yield: { expr: "5000.0" }Python
from nextstat import UnbinnedModel
model = UnbinnedModel.from_config("unbinned_spec.yaml")
result = model.fit(data="events.parquet")
# Profile scan over μ
scan = model.profile_scan("mu", bounds=(0.0, 3.0), n_points=40)
# CLs upper limit
limit = model.upper_limit("mu", cl=0.95)CLI
# Fit
nextstat unbinned-fit spec.yaml --data events.parquet
# Profile scan
nextstat unbinned-scan spec.yaml --poi mu --range 0:3 --points 40
# Toy studies on GPU
nextstat unbinned-fit-toys spec.yaml --ntoys 2000 --gpu cuda13 parametric PDFs out of the box
Gaussian, Crystal Ball, Double Crystal Ball, Exponential, Chebyshev, Bernstein, Voigtian, ARGUS, Spline, KDE (1D–3D), Product PDF. Analytical normalization where possible, Gauss-Legendre quadrature for the rest.
EventStore: SoA columnar layout
Observables stored in Structure-of-Arrays format for cache-friendly SIMD. Automatic ingestion from Parquet, ROOT TTree, or NumPy.
Declarative YAML specification
unbinned_spec_v0 — same approach as HistFactory JSON. Processes, PDFs, yields, rate modifiers (NormSys, WeightSys), per-event systematics — all in one file.
Unified LogDensityModel trait
Same optimizer (L-BFGS), same profile scan, same CLs — binned and unbinned stacks are fully compatible.
GPU acceleration (CUDA + Metal)
NLL reduction on GPU for thousands of events. Batch toy fitting entirely on device. Unbinned --gpu reads Parquet directly.
Templates out. Flows in.
Normalizing flows and DCR surrogates replace binned templates in unbinned analysis. Train a model in any framework, export ONNX — NextStat uses it as a full PDF with gradients, normalization, and sampling.
# unbinned_spec.yaml — neural PDF instead of template
processes:
- name: signal
pdf:
type: conditional_flow # neural flow
manifest: models/flow_manifest.json
context_params: [alpha_jes] # NP as context
yield: { expr: "mu * 1200.0" }
- name: background
pdf:
type: dcr_surrogate # neural surrogate
manifest: models/dcr_manifest.json
yield: { expr: "3400.0" }from nextstat import UnbinnedModel
model = UnbinnedModel.from_config("unbinned_spec.yaml")
result = model.fit()
# Profile scan over μ
scan = model.profile_scan("mu", bounds=(0.0, 3.0), n_points=40)Any framework → ONNX → NextStat
Train flows in PyTorch (zuko), JAX (distrax), or TensorFlow — export to ONNX, NextStat picks it up automatically.
DCR Surrogate: bin-free morphing
Neural Direct Classifier Ratio replaces template morphing. Smooth, continuous systematics — no bins, no interpolation.
FAIR-HUC: validated training
Distillation protocol from HistFactory templates with normalization checks (PIT/KS) and closure tests. The surrogate is verifiable.
Feature-gated: zero overhead
--features neural activates neural PDFs. Without the flag — the binary contains no ONNX Runtime, size and speed are unchanged.
Beyond NUTS. MAMS.
Metropolis-Adjusted Microcanonical Sampler (Robnik et al., 2025) — a microcanonical MCMC method that reports 2–7× ESS per gradient evaluation over NUTS on published benchmarks. NextStat ships a production-ready implementation with auto-tuning and multi-chain support.
import nextstat
# Any model that implements dim/nll/grad_nll
model = nextstat.EightSchoolsModel(y, sigma)
# MAMS: arXiv:2503.01707 — next-gen MCMC
result = nextstat.sample_mams(
model,
n_chains=4,
n_warmup=1000,
n_samples=2000,
seed=42,
)
# Same diagnostics as NUTS: R-hat, ESS, divergences
print(result.summary())Isokinetic dynamics on the unit sphere
Velocity lives on S^{d-1} with projected gradient updates. Fixed trajectory length — no tree building, predictable cost per transition.
Auto-tuned warmup
3-phase warmup: binary search for step size ε, Welford diagonal preconditioner, then decoherence length L tuning via ESS/gradient optimisation.
Overflow-safe MH criterion
Log-sum-exp formulation via ζ = exp(−2δ) avoids cosh/sinh overflow for large gradients — stable where NUTS diverges.
Drop-in NUTS replacement
Same Chain output, same diagnostics (R-hat, ESS, divergences), same ArviZ export. Switch one function call.
Churn analysis. Causal answers.
A complete subscription analytics pipeline for SaaS, telco, and insurance — from data ingestion to causal uplift estimation. Survival models, hazard ratios, AIPW treatment effects, and retention matrices, all compiled to native speed.
import nextstat
# Generate or ingest real customer data
data = nextstat.churn_generate_data(
n_customers=5000, seed=42
)
# Cohort retention: stratified KM + log-rank
ret = nextstat.churn_retention(
data.times, data.events, data.groups
)
# Cox PH risk model → hazard ratios + CIs
risk = nextstat.churn_risk_model(
data.times, data.events,
data.covariates, data.names
)
# Causal uplift: did the intervention reduce churn?
uplift = nextstat.churn_uplift(
data.times, data.events,
data.treated, data.covariates,
horizon=12.0,
)
print(f"ATE = {uplift.ate:.4f}, Γ = {uplift.gamma}")End-to-end pipeline
From Parquet/CSV ingestion with YAML column mapping through KM retention, Cox PH risk model, to AIPW causal uplift — one library, zero glue code.
Causal uplift with Rosenbaum sensitivity
AIPW-based treatment effect estimation with built-in sensitivity analysis. Know exactly how much hidden confounding would overturn your result.
Cohort retention matrix
Life-table style retention matrix by cohort — the exact artifact growth teams need. JSON-serialisable, ready for dashboards.
Bootstrap CIs via Rayon
Parallel bootstrap hazard ratios across all CPU cores. 1000 resamples in seconds, not minutes.
Population PK. Rust speed.
Pharmacokinetic and pharmacodynamic modeling with NLME estimation, NONMEM-compatible dataset reader, stepwise covariate modeling, and VPC diagnostics. The same L-BFGS + analytical gradients engine that powers HEP fits — applied to clinical pharmacology.
import nextstat
# 1-compartment oral PK model
model = nextstat.OneCompartmentOralPkModel(
times=times, observations=conc,
dose=320.0, tau=24.0,
)
# MLE fit → CL, V, KA estimates
fit = nextstat.fit(model)
print(fit.parameters) # [CL, V, KA]
# Population PK via NLME (Laplace approx.)
pop = nextstat.OneCompartmentOralPkNlmeModel(
subjects=subjects,
dose=320.0, tau=24.0,
)
pop_fit = nextstat.fit(pop)
# Visual Predictive Check
vpc = nextstat.vpc(pop, pop_fit, n_sim=500)1-compartment oral PK + population NLME
Individual and population PK models with Laplace approximation. Fixed and random effects on CL, V, KA.
FOCE / FOCEI estimation
First-order conditional estimation — the FDA gold standard for nonlinear mixed-effects models. Production-grade implementation.
NONMEM dataset reader
Read .csv with EVID/MDV/AMT/DV columns directly. Multi-dose regimens, infusions, and dose adjustments. No format conversion needed.
Stepwise covariate modeling
Automated forward addition + backward elimination with LRT. Identify which patient covariates explain PK variability.
VPC + diagnostics
Visual predictive checks, goodness-of-fit plots. The standard pharma toolkit for model qualification.
Zero-copy. Disk to GPU.
Parquet and Arrow IPC readers with mmap zero-copy, column projection, row group predicate pushdown, and parallel decode. Data flows from disk to fit — or straight to GPU — without intermediate allocations.
import nextstat
# Zero-copy mmap Parquet reader
model = nextstat.from_parquet(
"events.parquet",
spec="analysis.yaml",
)
# Or with column projection + row group pushdown
model = nextstat.from_parquet_with_modifiers(
"events.parquet",
spec="analysis.yaml",
columns=["mass", "pt", "eta"],
filter="pt > 25.0 && njet >= 4",
)
# Same API for Arrow IPC
model = nextstat.from_arrow_ipc(
"events.arrow",
spec="analysis.yaml",
)
# Fit — data never leaves mmap until needed
result = nextstat.fit(model)mmap zero-copy reads
memmap2-backed reader — no heap allocation for file I/O. 100K events with 4 columns decoded in 557 µs.
Column projection
Read only the columns your analysis needs. 35–50% speedup on wide tables by skipping irrelevant columns entirely.
Row group predicate pushdown
Prune row groups using min/max statistics before decoding. Up to 5× speedup on selective filters.
Parallel row group decode
Rayon-based parallel decompression. Benefits datasets with >1M events and many row groups.
GPU-direct upload
SoA f64 arrays go straight to CudaSlice<f64> (CUDA) or MTLBuffer (Metal, f32). No intermediate copies between disk and GPU.
R users. Welcome.
NextStat ships an R package built with extendr. HistFactory fits, GLMs, Kalman filters, GARCH, and stochastic volatility — all calling the same Rust core. CRAN-ready packaging with vendored crates.
library(nextstat)
# HistFactory fit from pyhf JSON
fit <- nextstat_fit("workspace.json")
print(fit$parameters)
# GLM: logistic regression
glm <- nextstat_glm_logistic(X, y)
print(glm$coefficients)
# Time series: Kalman filter + GARCH
kf <- nextstat_kalman(y, F, H, Q, R, x0, P0)
garch <- nextstat_garch(returns)
# Stochastic volatility
sv <- nextstat_sv(returns)
print(sv$volatility)11 functions, CRAN-ready
roxygen2 documentation, testthat test suite with 48 assertions, configure script, getting-started vignette, and NEWS.md. Ready for CRAN submission.
HEP + applied stats in one package
nextstat_fit() / nextstat_hypotest() / nextstat_upper_limit() for HistFactory alongside nextstat_glm_logistic() / nextstat_kalman() / nextstat_garch() for general statistics.
Vendored Rust crates
Self-contained tarball — no internet required during R CMD INSTALL. make nsr-vendor-sync keeps the bundled crate sources in sync.
Same engine, different language
Identical Rust core as the Python and WASM bindings. R users get the same speed and numerical accuracy — not a reimplementation.
The agent calls. NextStat answers.
GPT, Claude, Llama, a local model via Ollama — any agent can invoke NextStat as a tool. OpenAI function calling, LangChain, and MCP out of the box. Locally or through a GPU server.
OpenAI
from nextstat.tools import get_toolkit, execute_tool
import openai, json
tools = get_toolkit() # OpenAI-compatible definitions
resp = openai.chat.completions.create(
model="gpt",
messages=[{"role": "user", "content": "Fit the workspace and show POI"}],
tools=tools,
)
for call in resp.choices[0].message.tool_calls:
result = execute_tool(
call.function.name,
json.loads(call.function.arguments),
)
# result: { ok, result, error, meta }MCP
from nextstat.tools import get_mcp_tools, handle_mcp_call
mcp_tools = get_mcp_tools() # MCP-compatible definitions
# ... MCP server receives a call ...
result = handle_mcp_call("nextstat_fit", {"workspace_json": ws})LangChain
from nextstat.tools import get_langchain_tools
tools = get_langchain_tools() # list[StructuredTool]
agent = create_tool_calling_agent(llm, tools, prompt)8 tools — full analysis cycle
fit, hypotest, upper_limit, ranking, discovery, scan, workspace_audit, read_root_histogram. An agent can run a complete analysis from ROOT file to publication.
OpenAI · LangChain · MCP
Single tool registry with adapters for three major standards. get_toolkit() — OpenAI, get_langchain_tools() — LangChain StructuredTool, get_mcp_tools() — Model Context Protocol.
local + server transport
Tools run in-process via Python or remotely via nextstat-server (HTTP). Automatic fallback: if the server is unavailable, the call executes locally.
Deterministic envelope
Every response is a stable JSON nextstat.tool_result.v1: { ok, result, error, meta }. Parity mode + threads=1 guarantees reproducibility across calls.
Physics Assistant — demo
Full example: ROOT histograms → workspace → anomaly scan → discovery → upper limits. Runs locally, via server, or in Docker.
Performance
at scale.
Toy-based CLs on a real-scale HEP model. Same workspace, same physics. Wall-clock time.
Install in seconds.
Run in minutes.
pip install nextstat. No conda. No root. No compiler.
NextStat builds on the pioneering work of pyhf, TRExFitter, ROOT, and the HistFactory specification.
We are grateful to their authors and communities.