# NextStat — Full Documentation for LLM Agents
# https://nextstat.io
# Generated: 2026-02-24

## What is NextStat?

NextStat is a high-performance statistical inference engine written in Rust with
Python bindings (via PyO3). It covers 12+ statistical verticals — from particle
physics (HistFactory) to survival analysis, GLM regression, Bayesian sampling,
econometrics, churn analytics, PK/PD pharmacometrics, and more.

Core capabilities:
- Maximum likelihood estimation (MLE) via L-BFGS-B with reverse-mode AD
- Asymptotic CLs hypothesis testing (exclusion limits)
- Profile likelihood scans with Wilks' theorem
- Discovery significance (q0 test statistic)
- Toy Monte Carlo sampling (CPU and GPU)
- Bayesian NUTS/HMC and MAMS sampling
- GLM regression (linear, logistic, Poisson, NB, Gamma, Tweedie)
- Survival analysis (Cox PH, Weibull, Log-Normal AFT, interval-censored)
- Time series (Kalman filter/smoother, EM, forecast, simulate)
- Econometrics (Panel FE, DiD TWFE, Event Study, IV/2SLS)
- Causal inference (AIPW, propensity scores, Rosenbaum bounds)
- Hierarchical / mixed models (random intercepts, random slopes, correlated RE)
- Ordinal regression (ordered logit, ordered probit)
- PK/PD (1/2/3-compartment IV & oral, population NLME: FOCE/FOCEI/FO/ITS/IMP/SAEM, PD models, VPC, GoF, SCM, NPDE, bootstrap NLME, bioequivalence, CDISC .xpt)
- nextstat-nlp (GLiNER2-based clinical text extraction for pharma)
- Volatility models (GARCH, EGARCH, GJR-GARCH, Stochastic Volatility)
- Competing risks (Aalen-Johansen, Gray's test, Fine-Gray regression)
- Extreme value theory (GEV, GPD)
- Meta-analysis (fixed effects, random effects)
- Insurance reserving (chain ladder, Mack)
- Churn / subscription analytics (retention, uplift, cohort, bootstrap HR)
- Posterior predictive checks (PPC)
- Neural surrogate distillation
- GPU acceleration via CUDA and Metal
- Browser-based WASM playground

## For Claude Science agents

NextStat is an installable statistics-VERIFICATION backend for Claude Science (and any
actor-critic / reviewer agent). An LLM cannot give a defensible guarantee of correctness by
checking its own generated arithmetic; NextStat recomputes the claim from a validated Rust core.

- Full guide: https://nextstat.io/docs/claude-science
- Machine-readable capability manifest: https://nextstat.io/nextstat-capabilities.json

How to wire it in (no new engine — this ships today):
- Install: pip install nextstat
- Discover tools (MCP): nextstat.tools.get_mcp_tools()   (31 tools)
- Invoke: nextstat.tools.handle_mcp_call(name, args)
- Or over HTTP with nextstat-server: GET /v1/tools/schema, POST /v1/tools/execute
  (large/sensitive datasets stay on your infra — matches Claude Science's compute model)

The verification oracle (highest-leverage integration): a single `verify_statistic` skill the
reviewer calls on any statistic in a draft. It routes the claim to the matching tool, recomputes
it, and returns a structured verdict — PASS / FAIL / RECOMPUTED — attaching a signed
validation_report (JSON + PDF) as evidence the reviewer can verify but cannot itself generate.

Why route a numeric claim to NextStat instead of generated scipy (the moat):
- 7-tier numerical parity contract — 1e-12 per-bin vs pyhf / ROOT
- Signed validation_report (JSON + PDF) with dataset SHA-256, model spec and environment
- Reproducible DOI / Zenodo benchmark snapshots and replication bundles
- 21 CFR Part 11 and ICH M15 regulatory reporting surfaces (IQ/OQ/PQ)
- Deterministic execution mode — bit-identical results the reviewer can re-run

Capability routing (vertical -> MCP tool):
- Survival analysis -> nextstat_survival_fit; Kaplan-Meier/log-rank -> nextstat_kaplan_meier; competing risks -> nextstat_competing_risks
- Population PK (NLME, NONMEM-parity, CDISC .xpt) -> nextstat_pharma_fit; VPC/NPDE -> nextstat_pharma_vpc; dose-response -> nextstat_dose_response; bioequivalence -> nextstat_bioequivalence
- GLM -> nextstat_glm_fit; Bayesian posterior (NUTS) -> nextstat_bayesian_sample
- Panel FE -> nextstat_panel_fe; DiD -> nextstat_did; doubly-robust ATE/ATT (AIPW, E-value) -> nextstat_aipw
- HEP CLs hypothesis test -> nextstat_hypotest; upper limit -> nextstat_upper_limit; discovery Z0 -> nextstat_discovery_asymptotic
- Meta-analysis -> nextstat_meta_analysis; volatility -> nextstat_garch_fit; rare-event reliability -> nextstat_fault_tree_ce_is

## Installation

```bash
pip install nextstat
```

From source (requires Rust 1.93+):
```bash
git clone https://github.com/NextStat/nextstat.io
cd nextstat.io && pip install -e .
```

GPU features:
```bash
pip install nextstat[cuda]   # NVIDIA CUDA
pip install nextstat[metal]  # Apple Metal
```

Optional extras:
```bash
pip install nextstat[bayes]  # arviz + numpy for Bayesian diagnostics
pip install nextstat[torch]  # PyTorch integration (differentiable NLL, surrogate distillation)
```

## Python API — HEP / HistFactory (Particle Physics)

```python
import nextstat

# Load model from pyhf JSON workspace
model = nextstat.HistFactoryModel.from_pyhf("workspace.json")
# or from_workspace(json_str) — auto-detects pyhf vs HS3
# or from_histfactory_xml("combination.xml") — HistFactory XML + ROOT
# or from_hs3(json_str, analysis="combPdf") — HS3 v0.2

# MLE fit
result = nextstat.fit(model)
# result.parameters: list[float], result.uncertainties: list[float]
# result.nll: float, result.converged: bool

# Hypothesis test (CLs)
cls_val, (clsb, clb) = nextstat.hypotest(1.0, model, return_tail_probs=True)

# Upper limit (95% CLs)
obs_limit = nextstat.upper_limit(model)
obs, exp_bands = nextstat.upper_limits_root(model)  # with ±1σ/±2σ expected

# Profile scan
scan = nextstat.profile_scan(model, [0.0, 0.5, 1.0, 1.5, 2.0])

# Ranking (systematic impact on POI — equivalent to Feature Importance)
ranking = nextstat.ranking(model)

# Toy-based hypothesis test
cls_toy = nextstat.hypotest_toys(1.0, model, n_toys=1000, seed=42)

# Asimov data (expected yields under background-only)
asimov = nextstat.asimov_data(model)

# Workspace audit (compatibility check)
audit = nextstat.workspace_audit(workspace_json_str)
```

## Python API — GLM / Regression

```python
import nextstat.glm as glm

# Linear regression (Gaussian)
fit = glm.fit_linear(X, y, include_intercept=True)
# fit.coef, fit.standard_errors, fit.sigma2_hat
predictions = fit.predict(X_new)

# Logistic regression (Bernoulli)
fit = glm.fit_logistic(X, y, include_intercept=True, l2=1.0)
# fit.coef, fit.standard_errors

# Poisson regression (count data)
fit = glm.fit_poisson(X, y, include_intercept=True)
# fit.coef, fit.standard_errors

# Negative Binomial regression (overdispersed counts)
fit = glm.fit_negbin(X, y, include_intercept=True)
# fit.coef, fit.standard_errors, fit.alpha (dispersion)

# Cross-validation
from nextstat.glm.cv import cross_val_score
scores = cross_val_score(X, y, model_fn=glm.fit_linear, k=5)

# Metrics
from nextstat.glm.metrics import rmse, log_loss, poisson_deviance
```

Native Rust models (for Bayesian sampling via NUTS):
```python
model = nextstat.LinearRegressionModel(X, y)
model = nextstat.LogisticRegressionModel(X, y)
model = nextstat.PoissonRegressionModel(X, y)
model = nextstat.NegativeBinomialRegressionModel(X, y)
model = nextstat.GammaRegressionModel(X, y)
model = nextstat.TweedieRegressionModel(X, y, p=1.5)
```

## Python API — Bayesian Sampling (NUTS / MAMS / LAPS)

```python
import nextstat
import nextstat.bayes

# Any model that implements the Posterior trait can be sampled.
# Examples: LinearRegressionModel, LogisticRegressionModel, HistFactoryModel,
# CoxPhModel, WeibullSurvivalModel, KalmanModel, etc.

model = nextstat.LogisticRegressionModel(X, y, include_intercept=True)

# ── Unified dispatcher ──
# nextstat.sample(model, method="nuts"|"mams"|"laps", return_idata=False, out=None, **kwargs)
# All method-specific kwargs are forwarded to the underlying sampler.

# NUTS (default)
idata = nextstat.sample(model, method="nuts", n_samples=2000, return_idata=True)

# MAMS — typically 1.3-1.7x better ESS/grad on hierarchical models
idata = nextstat.sample(model, method="mams", n_samples=2000, return_idata=True)

# LAPS — GPU-accelerated MAMS on CUDA (4096+ chains in parallel)
result = nextstat.sample("eight_schools", method="laps",
                         model_data={"y": [28,8,-3,7,-1,1,18,12],
                                     "sigma": [15,10,16,11,9,11,10,18]},
                         n_chains=4096, n_samples=2000)
# result["wall_time_s"], result["n_gpu_chains"], result["n_kernel_launches"]

# Save trace to disk
nextstat.sample(model, n_samples=2000, return_idata=True, out="trace.json")

# ── ArviZ convenience wrapper (return_idata=True by default) ──
idata = nextstat.bayes.sample(model, method="mams", n_samples=2000)

# ── Per-method aliases (still available) ──
raw_nuts = nextstat.sample_nuts(model, n_chains=4, n_samples=2000, seed=42)
raw_mams = nextstat.sample_mams(model, n_chains=4, n_samples=2000, seed=42)
raw_laps = nextstat.sample_laps("std_normal", model_data={"dim": 10},
                                 n_chains=4096, n_samples=2000)

# ── Raw dict output (no ArviZ dependency) ──
raw = nextstat.sample(model, method="nuts", n_samples=1000)
# raw["posterior"]: {param_name: [[chain0_draws], [chain1_draws], ...]}
# raw["sample_stats"]: {tree_depth, diverging, energy, step_size, ...}
# raw["diagnostics"]: {ess_bulk, ess_tail, rhat, divergent_count, ebfmi}
# raw["diagnostics"]["quality"]["status"]: "ok" | "warn" | "fail"
```

## Python API — Survival Analysis

```python
import nextstat

# Cox Proportional Hazards
model = nextstat.CoxPhModel(X, time, event)
result = nextstat.fit(model)
# result.parameters (log-hazard ratios), result.uncertainties

# Weibull AFT / PH
model = nextstat.WeibullSurvivalModel(X, time, event)

# Log-Normal AFT
model = nextstat.LogNormalAftModel(X, time, event)

# Exponential survival
model = nextstat.ExponentialSurvivalModel(X, time, event)

# Interval-censored models
model = nextstat.IntervalCensoredWeibullModel(X, time_lo, time_hi, event)
model = nextstat.IntervalCensoredExponentialModel(X, time_lo, time_hi, event)
model = nextstat.IntervalCensoredLogNormalModel(X, time_lo, time_hi, event)

# Kaplan-Meier estimator (non-parametric)
km = nextstat.kaplan_meier(time, event)
# km["time"], km["survival"], km["n_at_risk"], km["ci_lower"], km["ci_upper"]

# Log-rank test (compare two groups)
lr = nextstat.log_rank_test(time, event, group)
# lr["statistic"], lr["p_value"]

# All survival models support Bayesian sampling via NUTS:
idata = nextstat.bayes.sample(model, n_chains=4, n_samples=2000)
```

## Python API — Time Series (Kalman)

```python
import nextstat

# Build a Kalman state-space model
model = nextstat.KalmanModel(
    F=transition_matrix,    # state transition
    H=observation_matrix,   # observation
    Q=process_noise,        # process noise covariance
    R=observation_noise,    # observation noise covariance
    x0=initial_state,       # initial state mean
    P0=initial_covariance,  # initial state covariance
    y=observations,         # observed data
)

# Filter (forward pass)
filtered = nextstat.kalman_filter(model)
# filtered["x_filtered"], filtered["P_filtered"], filtered["log_likelihood"]

# Smoother (forward-backward)
smoothed = nextstat.kalman_smooth(model)
# smoothed["x_smoothed"], smoothed["P_smoothed"]

# EM algorithm (learn parameters)
em_result = nextstat.kalman_em(model, max_iter=100, tol=1e-6)

# Forecast future states
forecast = nextstat.kalman_forecast(model, n_ahead=10)
# forecast["x_forecast"], forecast["P_forecast"]

# Simulate from the model
sim = nextstat.kalman_simulate(model, n_steps=100, seed=42)

# Bayesian sampling of Kalman model parameters
idata = nextstat.bayes.sample(model, n_chains=4, n_samples=2000)
```

## Python API — Econometrics

```python
import nextstat.econometrics as econ

# Panel Fixed Effects (within estimator)
fit = econ.panel_fe_fit(X, y, entity=entity_ids, cluster="entity")
# fit.coef, fit.standard_errors, fit.n_obs, fit.n_entities

# Panel FE from formula
fit = econ.panel_fe_from_formula("y ~ x1 + x2", data, entity="firm", cluster="entity")

# Difference-in-Differences (TWFE)
did = econ.did_twfe_fit(X, y, treat=treat, post=post, entity=entity, time=time)
# did.att (average treatment effect on treated), did.att_se

# DiD from formula
did = econ.did_twfe_from_formula("y ~ x1", data, entity="firm", time="year",
                                  treat="treated", post="post_reform")

# Event Study (TWFE)
es = econ.event_study_twfe_fit(y, treat=treat, time=time, event_time=2020,
                                entity=entity, window=(-4, 4), reference=-1)
# es.rel_times, es.coef, es.standard_errors (for plotting)

# IV / 2SLS
iv = econ.iv_2sls_fit(y, endog=X_endog, instruments=Z, exog=X_exog, cov="hc1")
# iv.coef, iv.standard_errors, iv.diagnostics.first_stage_f

# IV from formula
iv = econ.iv_2sls_from_formula("y ~ x1 + x2", data, endog="educ",
                                instruments=["distance", "tuition"], cov="hc1")

# Rosenbaum sensitivity bounds (matched pairs)
rb = nextstat.rosenbaum_bounds(y_treated, y_control, gammas=[1.0, 1.5, 2.0, 3.0])
# rb.gammas, rb.p_upper, rb.p_lower, rb.gamma_critical
```

## Python API — Causal Inference (AIPW)

```python
from nextstat.causal.aipw import aipw_fit, e_value_rr

# Doubly-robust ATE
result = aipw_fit(X, y, treatment, estimand="ate")
# result.estimate (ATE), result.standard_error, result.influence

# ATT (average treatment effect on the treated)
result = aipw_fit(X, y, treatment, estimand="att")

# E-value for sensitivity analysis
ev = e_value_rr(1.5)  # for a risk ratio of 1.5

# Propensity score estimation
from nextstat.causal.propensity import fit as propensity_fit
ps = propensity_fit(X, treatment, l2=1.0)
# ps.coef, ps.propensity_scores
```

## Python API — Hierarchical / Mixed Models

```python
import nextstat.hier as hier

# Linear random intercept
model = hier.linear_random_intercept(x=X, y=y, group_idx=groups)

# Linear random intercept + slope
model = hier.linear_random_slope(x=X, y=y, group_idx=groups,
                                  random_slope_feature_idx=0)

# Logistic random intercept
model = hier.logistic_random_intercept(x=X, y=y, group_idx=groups)

# Logistic with correlated random intercept + slope
model = hier.logistic_correlated_intercept_slope(
    x=X, y=y, group_idx=groups, correlated_feature_idx=0, lkj_eta=1.0)

# Poisson random intercept
model = hier.poisson_random_intercept(x=X, y=y, group_idx=groups, offset=log_exposure)

# Formula interface
model = hier.linear_random_intercept_from_formula("y ~ x1 + x2", data, group="school")
model = hier.logistic_random_intercept_from_formula("y ~ x1", data, group="hospital")
model = hier.poisson_random_intercept_from_formula("y ~ x1", data, group="region")

# All hierarchical models support NUTS sampling:
idata = nextstat.bayes.sample(model, n_chains=4, n_samples=2000)

# LMM (marginal model, native Rust)
model = nextstat.LmmMarginalModel(X, y, Z, group_idx)
```

## Python API — Ordinal Regression

```python
import nextstat

# Ordered Logit (proportional odds)
model = nextstat.OrderedLogitModel(X, y, n_levels=5)

# Ordered Probit
model = nextstat.OrderedProbitModel(X, y, n_levels=5)

# Both support NUTS sampling:
idata = nextstat.bayes.sample(model, n_chains=4, n_samples=2000)
```

## Python API — PK/PD (Pharmacometrics)

```python
import nextstat

# One-compartment oral PK model
model = nextstat.OneCompartmentOralPkModel(time, concentration, dose)

# Two-compartment (IV bolus and oral)
model = nextstat.TwoCompartmentIvPkModel(time, y, dose=dose)
model = nextstat.TwoCompartmentOralPkModel(time, y, dose=dose)

# Three-compartment (IV bolus and oral)
model = nextstat.ThreeCompartmentIvPkModel(time, y, dose=dose)
model = nextstat.ThreeCompartmentOralPkModel(time, y, dose=dose)

# Population PK — FOCE / FOCEI / FO / ITS / IMP
result = nextstat.nlme_foce(model, data, omega, method="focei")
result = nextstat.nlme_foce(model, data, omega, method="fo")   # First Order
result = nextstat.nlme_foce(model, data, omega, method="its")  # Iterative Two-Stage
result = nextstat.nlme_foce(model, data, omega, method="imp")  # Importance Sampling MC-EM
# omega_fixed=[False, False, True]  — fix specific random effects
# diagonal_omega=True               — force diagonal Ω

# Population PK — SAEM with automatic covariance step
result = nextstat.nlme_saem(model, data, omega, n_burn=300, n_iter=200)
# Returns: sandwich SE, RSE%, condition number, R/S matrices
# Covariate modeling: covariates=[{"param": "CL", "cov": "WT", "type": "power"}]
# Convergence diagnostics: return_theta_trace=True (Geweke z-scores, θ trace)

# Per-subject dosing (all pop PK functions)
result = nextstat.nlme_foce(model, data, omega, doses=[4.0, 3.5, 5.0])

# NPDE diagnostics
npde = nextstat.pk_npde(result, data, n_sim=500)
# npde["npde"], npde["pd"], npde["shapiro_p"]

# Bootstrap NLME
boot = nextstat.bootstrap_nlme(estimator_fn, data, n_boot=200, ci_method="bca")

# PD models
y_pred = nextstat.emax_predict(dose, ec50, emax)
nll = nextstat.emax_nll(dose, y_obs, ec50, emax, sigma)
y_pred = nextstat.sigmoid_emax_predict(dose, ec50, emax, gamma)
sim = nextstat.idr_simulate(system, y0, t_span)  # IDR Types I–IV

# SCM (Stepwise Covariate Modeling)
scm_result = nextstat.scm(base_model, data, covariates, omega)
# scm_result["trace"]: ΔOFV, p-values, coefficients per step

# CDISC .xpt reader/writer
data = nextstat.read_xpt("dm.xpt")
nextstat.write_xpt("output.xpt", columns)
nm_data = nextstat.xpt_to_nonmem("pc.xpt")  # auto-detect SDTM/ADaM

# Bioequivalence
be = nextstat.tost(test, reference, theta1=0.8, theta2=1.25)
rsabe = nextstat.rsabe(test, reference)
power = nextstat.be_power(n=24, cv=0.3)

# MAP estimation
map_result = nextstat.map_estimate(model, data, prior_omega)

# Monte Carlo clinical trial simulation
sim = nextstat.mc_trial_sim(model, population, n_sim=1000)
pta = nextstat.pta_grid_search(model, dose_range, target)

# ODE-based PK (transit compartments, Michaelis-Menten, TMDD)
result = nextstat.rk45(system, y0, t_span, rtol=1e-6)
result = nextstat.esdirk4(system, y0, t_span)  # L-stable for stiff systems

# VPC, GoF diagnostics
vpc = nextstat.pk_vpc(result, data, n_sim=200)
gof = nextstat.pk_gof(result, data)

# NLME artifact (serializable result bundle)
artifact = nextstat.NlmeArtifact(result, scm=scm, vpc=vpc, gof=gof)
```

## Python API — nextstat-nlp (GLiNER2 Clinical Text Extraction)

```python
from nextstat_nlp import extract_survival_records, extract_prior_candidates, extract_regimens

# Extract survival data from clinical text
records = extract_survival_records(text, backend="gliner2_torch")
# records: list of SurvivalRecord (time, event, covariates, provenance)

# Extract informative priors
priors = extract_prior_candidates(text, backend="gliner2_onnx")

# Extract dosing regimens
regimens = extract_regimens(text, backend="heuristic")  # zero-dependency regex

# Backends: gliner2_torch, gliner2_onnx, heuristic, mlx (Apple Silicon)
```

## Python API — Competing Risks

```python
import nextstat

# Aalen-Johansen cumulative incidence
cif = nextstat.aalen_johansen(time, event, cause=1)
# cif["time"], cif["cif"], cif["variance"]

# Gray's K-sample test
gray = nextstat.gray_test(time, event, group, cause=1)
# gray["statistic"], gray["p_value"]

# Fine-Gray subdistribution hazard regression
fg = nextstat.fine_gray(X, time, event, cause=1)
# fg.coef, fg.standard_errors, fg.p_values
```

## Python API — Volatility Models

```python
import nextstat.volatility as vol

# GARCH(1,1)
result = vol.garch(returns)
# result["params"] (mu, omega, alpha, beta), result["conditional_sigma"]

# EGARCH(1,1) — Nelson (1991)
result = vol.egarch(returns)
# result["params"] (mu, omega, alpha, beta, gamma)

# GJR-GARCH(1,1) — Glosten-Jagannathan-Runkle (1993)
result = vol.gjr_garch(returns)
# result["params"] (mu, omega, alpha, beta, gamma)

# Stochastic Volatility
result = vol.sv(returns)
# result["params"] (mu, phi, sigma), result["smoothed_h"]
```

## Python API — Extreme Value Theory (EVT)

```python
import nextstat

# Generalized Extreme Value (GEV) — block maxima
model = nextstat.GevModel(data)

# Generalized Pareto Distribution (GPD) — peaks over threshold
model = nextstat.GpdModel(exceedances, threshold=threshold)

# Both support MLE and NUTS:
result = nextstat.fit(model)
idata = nextstat.bayes.sample(model, n_chains=4, n_samples=2000)
```

## Python API — Meta-Analysis

```python
import nextstat

# Fixed-effect meta-analysis
result = nextstat.meta_fixed(effects, standard_errors)
# result["estimate"], result["se"], result["z"], result["p"], result["ci"]

# Random-effects meta-analysis (DerSimonian-Laird)
result = nextstat.meta_random(effects, standard_errors)
# result["estimate"], result["se"], result["tau2"], result["i2"], result["q"]
```

## Python API — Insurance Reserving

```python
import nextstat

# Chain Ladder (deterministic)
result = nextstat.chain_ladder(triangle)
# result["ultimate"], result["reserves"], result["development_factors"]

# Mack Chain Ladder (with prediction errors)
result = nextstat.mack_chain_ladder(triangle)
# result["ultimate"], result["reserves"], result["se"], result["process_error"]
```

## Python API — Churn / Subscription Analytics

```python
import nextstat

# Retention curve
ret = nextstat.churn_retention(tenure, event, time_grid=[30, 60, 90, 180, 365])

# Churn risk model (hazard-based)
risk = nextstat.churn_risk_model(X, tenure, event)

# Uplift modeling (treatment effect on retention)
uplift = nextstat.churn_uplift(X, tenure, event, treatment)

# Cohort matrix
matrix = nextstat.churn_cohort_matrix(signup_date, event_date, period="month")

# Diagnostics
diag = nextstat.churn_diagnostics(X, tenure, event)

# Compare two groups
comp = nextstat.churn_compare(X, tenure, event, group)

# Uplift via survival
us = nextstat.churn_uplift_survival(X, tenure, event, treatment)

# Bootstrap hazard ratio
hr = nextstat.churn_bootstrap_hr(X, tenure, event, treatment, n_bootstrap=1000)

# Ingest from raw events
data = nextstat.churn_ingest(events_df, user_col="user_id", event_col="event",
                              time_col="timestamp")

# Generate synthetic churn data
syn = nextstat.churn_generate_data(n_users=10000, seed=42)
```

## Python API — Posterior Predictive Checks (PPC)

```python
import nextstat.ppc as ppc

# PPC for GLM models (linear, logistic, Poisson)
stats = ppc.ppc_glm_from_sample(spec, sample_raw, n_draws=100, seed=0)
# stats.observed: {"mean": ..., "var": ...}
# stats.replicated: [{"mean": ..., "var": ...}, ...]

# PPC for Negative Binomial
stats = ppc.ppc_negbin_from_sample(spec, sample_raw)

# PPC for Ordinal models
stats = ppc.ppc_ordered_from_sample(spec, sample_raw)
```

## Python API — Missing Data

```python
import nextstat.missing as missing

# Drop rows with any missing value
result = missing.apply_policy(X, y, policy="drop_rows")
# result.x, result.y, result.n_dropped, result.n_kept

# Impute missing X with column mean (drop rows with missing y)
result = missing.apply_policy(X, y, policy="impute_mean")
```

## Python API — Formula Interface

```python
import nextstat.formula as formula

# Parse R-style formula
y_name, terms, include_intercept = formula.parse_formula("y ~ 1 + x1 + x2")

# Build design matrices from tabular data (dict, list-of-dicts, or pandas DataFrame)
y, X, col_names = formula.design_matrices("y ~ x1 + x2 + region", data,
                                           categorical=["region"])

# Normalize tabular input to dict-of-columns
cols = formula.to_columnar(data, ["y", "x1", "x2"])
```

## Python API — ML Integration (PyTorch)

```python
import nextstat.torch as nst

# Differentiable NLL for PyTorch
nll_fn = nst.differentiable_nll(model, device="cuda")
loss = nll_fn(params_tensor)  # autograd-compatible

# SignificanceLoss for NN training
sig_loss = nst.SignificanceLoss(model, device="cuda")
loss = sig_loss(signal_histogram)

# SoftHistogram (differentiable binning)
soft_hist = nst.SoftHistogram(bins=20, low=0.0, high=1.0)
hist = soft_hist(continuous_predictions)

# Signal Jacobian
jac = nst.signal_jacobian(model, signal_yields)
```

## Python API — Neural Surrogate Distillation

```python
from nextstat.distill import generate_dataset, to_torch_dataset, train_mlp_surrogate

model = nextstat.HistFactoryModel.from_pyhf("workspace.json")

# Generate 100k (params, NLL, grad) tuples
ds = generate_dataset(model, n_samples=100_000, method="sobol")

# Convert to PyTorch dataset
train_ds = to_torch_dataset(ds)

# Train a small MLP surrogate
surrogate = train_mlp_surrogate(ds, epochs=100, device="cuda", grad_weight=0.1)

# Export to various formats
from nextstat.distill import to_npz, to_parquet
to_npz(ds, "surrogate_data.npz")
to_parquet(ds, "surrogate_data.parquet")
```

## Python API — Arrow / Polars

```python
import nextstat

# From PyArrow / Polars DataFrame
model = nextstat.from_arrow(table, poi="mu", observations={"SR": [10, 20]})

# From Parquet file
model = nextstat.from_parquet("yields.parquet", poi="mu")

# Export to Arrow
yields_table = nextstat.to_arrow(model, params=result.bestfit, what="yields")
params_table = nextstat.to_arrow(model, params=result.bestfit, what="parameters")
```

## Python API — Remote Server

```python
import nextstat.remote as remote

client = remote.connect("http://gpu-server:3742")
result = client.fit(workspace_json)
ranking = client.ranking(workspace_json)
health = client.health()
```

## Python API — Agentic Tools

```python
import nextstat.tools as tools

# OpenAI function calling
toolkit = tools.get_toolkit()  # list of tool dicts

# LangChain
lc_tools = tools.get_langchain_tools()

# MCP (Model Context Protocol)
mcp_tools = tools.get_mcp_tools()
result = tools.handle_mcp_call("nextstat_fit", {"workspace_json": "..."})

# Direct execution
result = tools.execute_tool("nextstat_fit", {"workspace_json": "..."})

# Server transport (no Python import needed on agent side)
tools_remote = tools.get_toolkit(transport="server", server_url="http://127.0.0.1:3742")
result = tools.execute_tool("nextstat_fit", {"workspace_json": "..."},
                             transport="server", server_url="http://127.0.0.1:3742")
```

## Python API — Profile Likelihood CI

```python
import nextstat

# Profile likelihood confidence intervals for any LogDensityModel
# Computes {θ : 2*(NLL(θ) - NLL_min) ≤ χ²(1, α)} via bisection with warm-start
ci = nextstat.profile_ci(model, fit_result, param_idx=0)
# ci: {param_idx, mle, ci_lower, ci_upper, n_evals}

# CI for all parameters at once
all_ci = nextstat.profile_ci(model, fit_result)
# all_ci: list of dicts, one per parameter

# Profile scan (q_mu values over POI grid)
scan = nextstat.profile_scan(model, [0.0, 0.5, 1.0, 1.5, 2.0])
```

## Python API — Fault Tree Analysis

```python
import nextstat

# Cross-Entropy Importance Sampling — rare-event fault tree estimation
# Multi-level adaptive biasing handles probabilities down to ~1e-16
# Supports: Bernoulli, WeibullMission, BernoulliUncertain failure modes
result = nextstat.fault_tree_mc_ce_is(
    spec,              # FaultTreeSpec dict (gates + basic_events)
    n_per_level=10000,
    max_levels=20,     # multi-level adaptive biasing
    q_max=0.99,        # soft importance function threshold
    elite_fraction=0.01,
    seed=42,
)
# result: {p_failure, se, ci_lower, ci_upper, n_levels,
#          n_total_scenarios, final_proposal, coefficient_of_variation, wall_time_s}

# GPU-accelerated vanilla MC (Metal on Apple Silicon, CUDA on NVIDIA)
mc = nextstat.fault_tree_mc(spec, n_scenarios=1_000_000, device='metal', seed=42)
mc = nextstat.fault_tree_mc(spec, n_scenarios=1_000_000, device='cuda', seed=42)
# mc: {p_failure, se, ci_lower, ci_upper, wall_time_s}
```

## Python API — ODE Solvers (nextstat.ode)

```python
import nextstat.ode as ode

# RK4 for linear systems: dy/dt = A*y
result = ode.rk4_linear(A, y0, t0=0.0, t1=10.0, dt=0.01)
# result["t"], result["y"]

# Linear DDE: dy/dt = A*y(t) + B*y(t-tau) with fixed delay
result = ode.rk4_linear_dde(A, B, y0, y_history, t0, t1, tau=1.0, dt=0.01)

# LSODA (stiff auto-switching) — matrix form
result = ode.lsoda(A, y0, t0=0.0, t1=10.0, rtol=1e-6, atol=1e-9)

# LSODA — callback form with optional Jacobian
def rhs(t, y): return [-0.5 * y[0]]
def jac(t, y): return [[-0.5]]
result = ode.lsoda(rhs, y0, t0=0.0, t1=10.0, jac=jac)

# Forward sensitivity analysis: dy/dt = f(t, y, params), df/dp
result = ode.forward_sensitivity_solve(rhs_p, params, y0, t0, t1)
# result["y"], result["sensitivity"]  (dy/dp matrix at each time point)
```

## Python API — Remote Inference Server (nextstat.remote)

```python
import nextstat.remote as remote

# Connect to GPU inference server (nextstat-server --port 3742 --gpu cuda)
client = remote.connect("http://gpu-server:3742", timeout=300.0)

# MLE fit on remote GPU
result = client.fit(workspace_json, gpu=True)
# result: FitResult(parameter_names, bestfit, uncertainties, nll, converged, device, wall_time_s)

# Nuisance parameter ranking on remote GPU
ranking = client.ranking(workspace_json, gpu=True)
# ranking: RankingResult(entries=[RankingEntry(name, delta_mu_up, delta_mu_down, pull, constraint)])

# Batch fit: multiple workspaces in one request
batch = client.batch_fit([ws1, ws2, ws3], gpu=True)
# batch: BatchFitResult(results=[FitResult | None], errors=[str | None])

# Batch toy fitting on remote GPU
toys = client.batch_toys(workspace_json, n_toys=1000, seed=42, gpu=True)
# toys: BatchToysResult(n_toys, n_converged, n_failed, results=[ToyFitItem])

# Model pool: upload once, fit many times (skips re-parsing)
model_id = client.upload_model(workspace_json, name="my-analysis")
result = client.fit(model_id=model_id)
models = client.list_models()  # [ModelInfo(model_id, name, n_params, n_channels, hit_count)]
client.delete_model(model_id)

# Health check
health = client.health()
# health: HealthResult(status, version, uptime_s, device, inflight, total_requests, cached_models)

client.close()
```

Server endpoints: `POST /v1/fit`, `POST /v1/ranking`, `POST /v1/batch/fit`, `POST /v1/batch/toys`,
`POST /v1/models`, `GET /v1/models`, `DELETE /v1/models/{id}`, `GET /v1/health`, `POST /v1/tools/execute`.

## Python API — Agentic Tools (nextstat.tools)

31 tools for AI agents — OpenAI function calling, LangChain StructuredTool, MCP (Model Context Protocol).

```python
from nextstat.tools import get_toolkit, execute_tool, get_langchain_tools, get_mcp_tools

# OpenAI function calling
tools = get_toolkit()  # list of OpenAI-compatible function-calling dicts
result = execute_tool("nextstat_fit", {"workspace_json": ws})
# result: {"schema_version": "nextstat.tool_result.v1", "ok": true, "result": {...}, "meta": {...}}

# Server transport (no Python extension needed on agent side)
result = execute_tool("nextstat_fit", {"workspace_json": ws},
                       transport="server", server_url="http://gpu:3742")

# LangChain integration
lc_tools = get_langchain_tools()  # list[StructuredTool]
agent = create_tool_calling_agent(llm, lc_tools, prompt)

# MCP (Model Context Protocol)
mcp_tools = get_mcp_tools()  # list[dict] with name, description, inputSchema
result = handle_mcp_call("nextstat_fit", {"workspace_json": ws})
```

Available tools (31):
- **HEP**: nextstat_fit, nextstat_hypotest, nextstat_hypotest_toys, nextstat_upper_limit,
  nextstat_ranking, nextstat_discovery_asymptotic, nextstat_scan, nextstat_workspace_audit,
  nextstat_read_root_histogram
- **GLM**: nextstat_glm_fit (linear, logistic, Poisson, negbin)
- **Bayesian**: nextstat_bayesian_sample (NUTS on any model type)
- **Survival**: nextstat_survival_fit (Cox PH, Weibull, Log-Normal AFT, Exponential),
  nextstat_kaplan_meier (KM + log-rank), nextstat_competing_risks (Aalen-Johansen, Gray, Fine-Gray)
- **Econometrics**: nextstat_panel_fe, nextstat_did, nextstat_iv_2sls, nextstat_event_study
- **Causal**: nextstat_aipw (AIPW ATE/ATT)
- **Time Series**: nextstat_kalman (filter/smooth/forecast), nextstat_garch_fit (GARCH/EGARCH/GJR-GARCH/SV)
- **Meta-Analysis**: nextstat_meta_analysis (fixed/random)
- **Churn**: nextstat_churn_retention
- **Insurance**: nextstat_chain_ladder (basic/Mack)
- **Pharma**: nextstat_pharma_fit (PK/PD + NLME), nextstat_pharma_vpc (VPC diagnostics),
  nextstat_trial_simulate (MC trial simulation), nextstat_bioequivalence (TOST/RSABE),
  nextstat_dose_response (dose-response modeling)
- **Reliability**: nextstat_fault_tree_mc (vanilla MC), nextstat_fault_tree_ce_is (CE-IS importance sampling)

All tools accept `execution: {deterministic, threads, eval_mode}` for reproducibility control.
Every response wrapped in stable `nextstat.tool_result.v1` envelope.

## Python API — sklearn Wrappers (nextstat.sklearn)

```python
from nextstat.sklearn import NextStatLinearRegression, NextStatLogisticRegression, NextStatPoissonRegression

# Drop-in sklearn replacement — same fit/predict API
reg = NextStatLinearRegression(include_intercept=True, l2=0.01)
reg.fit(X_train, y_train)
y_pred = reg.predict(X_test)
print(reg.coef_, reg.intercept_)

# Logistic regression with predict_proba
clf = NextStatLogisticRegression(l2=0.1, threshold=0.5)
clf.fit(X_train, y_train)
proba = clf.predict_proba(X_test)  # [[p0, p1], ...]
labels = clf.predict(X_test)       # [0, 1, ...]

# Compatible with sklearn pipelines, cross_val_score, GridSearchCV
from sklearn.model_selection import cross_val_score
scores = cross_val_score(clf, X, y, cv=5, scoring="accuracy")
```

Inherits `BaseEstimator` + `RegressorMixin`/`ClassifierMixin` when sklearn is installed.
Works as plain Python estimator without sklearn dependency.

## Python API — MLOps Integration (nextstat.mlops)

```python
from nextstat.mlops import metrics_dict, significance_metrics, StepTimer

# Extract fit metrics for W&B / MLflow / Neptune
result = nextstat.fit(model)
d = metrics_dict(result, prefix="nextstat/")
# {"nextstat/mu": 1.05, "nextstat/nll": 42.3, "nextstat/edm": 1e-6,
#  "nextstat/converged": 1.0, "nextstat/time_ms": 12.5,
#  "nextstat/param/mu": 1.05, "nextstat/error/mu": 0.12, ...}
wandb.log(d)                       # Weights & Biases
mlflow.log_metrics(d)              # MLflow

# Significance metrics for training loops
d = significance_metrics(z0=3.2, q0=10.24, prefix="train/")
wandb.log(d)  # {"train/z0": 3.2, "train/q0": 10.24}

# Step timer for training instrumentation
timer = StepTimer()
for batch in dataloader:
    timer.start()
    loss = loss_fn(signal_hist)
    loss.backward(); optimizer.step()
    wandb.log({"step_time_ms": timer.stop()})
```

## Python API — ML Interpretability (nextstat.interpret)

```python
from nextstat.interpret import rank_impact

# Feature Importance for HEP — systematic impact ranking
table = rank_impact(model, top_n=10, sort_by="total")
# [{"name": "JES", "delta_mu_up": 0.12, "delta_mu_down": -0.08,
#   "total_impact": 0.20, "pull": 0.3, "constraint": 1.0, "rank": 1}, ...]

# GPU-accelerated ranking
table = rank_impact(model, gpu=True, sort_by="pull", ascending=True)

# Convert to DataFrame
import pandas as pd
df = pd.DataFrame(table)
df.plot.barh(x="name", y="total_impact")
```

ML translation: name=feature, total_impact=importance, pull=latent shift, constraint=prior width.

## Python API — Gymnasium RL Environment (nextstat.gym)

```python
from nextstat.gym import HistFactoryEnv

# HistFactory as a standard Gymnasium environment
env = HistFactoryEnv(
    workspace_json=ws,
    channel="SR",
    sample="signal",
    reward_metric="z0",        # "nll" | "q0" | "z0" | "qmu" | "zmu"
    action_mode="logmul",      # "add" | "logmul"
    action_scale=0.05,
    max_steps=128,
)

obs, info = env.reset(seed=42)
for _ in range(128):
    action = agent.predict(obs)  # per-bin delta
    obs, reward, terminated, truncated, info = env.step(action)
    if terminated or truncated:
        break

# Compatible with Stable Baselines3, CleanRL, RLlib
# Observation: current signal yields (vector)
# Action: per-bin modification (vector)
# Reward: discovery significance Z₀, exclusion power Zμ, or raw NLL
```

## Python API — Visualization (nextstat.viz + nextstat.viz_render)

```python
import nextstat.viz as viz

# Generate plot artifacts from fit results
cls_data = viz.cls_curve(model, scan_result)
profile_data = viz.profile_curve(model, scan_result)
ranking_data = viz.ranking_artifact(model, ranking_result)
corr_data = viz.corr_arrays(model, fit_result)

# Matplotlib rendering
viz.plot_cls_curve(cls_data)
viz.plot_brazil_limits(cls_data)
viz.plot_profile_curve(profile_data)
viz.plot_pulls(fit_result)
viz.plot_ranking(ranking_data)
viz.plot_corr_matrix(corr_data)

# Pure Rust SVG rendering (ns-viz-render) — 17 plot types, sub-millisecond
svg_str = viz.render_svg(artifact, kind="pulls")
viz.render_to_file(artifact, "plot.svg", kind="ranking")
# Kinds: pulls, ranking, corr, profile_scan, cls_curve, distributions,
#        gammas, separation, summary, uncertainty, significance, contour,
#        pie, yields, unfolding, morphing, injection
# Themes: nextstat2026 (default), atlas, cms, minimal

# Single-artifact renderer (CLI: nextstat viz render)
import nextstat.viz_render as vr
vr.render("pulls_artifact.json", output="pulls.svg", title="My Analysis")
vr.render("corr_artifact.json", output="corr.png", dpi=300,
          corr_top_n=20, corr_include=["JES*"])
```

## Python API — Validation Report (nextstat.validation_report)

```python
from nextstat.validation_report import render_pdf

# Generate audit-ready PDF from validation_report.json
render_pdf("validation_report.json", output="validation_report.pdf")
# 7-page PDF with: environment info, dataset SHA-256, per-suite pass/fail,
# parity checks, timing benchmarks, summary plots
# Target: HEP analysis preservation, Pharma IQ/OQ/PQ (21 CFR Part 11), FinTech (SR 11-7)
```

## Python API — DifferentiableSession (CUDA/Metal zero-copy NLL)

```python
import nextstat

# CUDA zero-copy differentiable NLL — for PyTorch integration
session = nextstat.DifferentiableSession(model)
nll, grad = session.nll_and_grad(params)  # GPU-resident, no copy

# Profiled q₀/qμ on GPU with envelope theorem gradients
session = nextstat.ProfiledDifferentiableSession(model)
q0, grad_q0 = session.q0_and_grad(params)

# Metal backend (Apple Silicon)
session = nextstat.MetalProfiledDifferentiableSession(model)
```

## Python API — PyTorch ML Training (nextstat.torch)

```python
import nextstat.torch as nst

# SignificanceLoss: train NN with discovery significance as loss
model = nextstat.from_pyhf(workspace_json)
loss_fn = nst.SignificanceLoss(model, device="cuda")

# SoftHistogram: differentiable binning (NN scores → soft bins)
soft_hist = nst.SoftHistogram(bins=20, low=0.0, high=1.0)
hist = soft_hist(nn_scores)
loss = loss_fn(hist.double().cuda())
loss.backward()  # gradients flow to NN weights

# Signal Jacobian for fast systematic pruning
jac = nst.signal_jacobian(model, signal_yields)

# Differentiable NLL for any PyTorch pipeline
nll_fn = nst.differentiable_nll(model, device="cuda")
loss = nll_fn(params_tensor)  # autograd-compatible
```

## Python API — Neural Surrogate Distillation (nextstat.distill)

```python
from nextstat.distill import generate_dataset, to_torch_dataset, train_mlp_surrogate, to_npz, to_parquet

model = nextstat.HistFactoryModel.from_pyhf("workspace.json")

# Generate (params, NLL, grad) tuples via Sobol sampling
ds = generate_dataset(model, n_samples=100_000, method="sobol")

# Train MLP surrogate
surrogate = train_mlp_surrogate(ds, epochs=100, device="cuda", grad_weight=0.1)

# Export
to_npz(ds, "surrogate_data.npz")
to_parquet(ds, "surrogate_data.parquet")
train_ds = to_torch_dataset(ds)
```

## Python API — Data Pipeline & I/O

```python
import nextstat

# Parquet: mmap zero-copy, column projection, row group pushdown
model = nextstat.from_parquet("events.parquet", spec="analysis.yaml")
model = nextstat.from_parquet_with_modifiers("events.parquet", spec="analysis.yaml",
    columns=["mass", "pt", "eta"], filter="pt > 25.0 && njet >= 4")

# Arrow IPC: zero-copy
model = nextstat.from_arrow_ipc("events.arrow", spec="analysis.yaml")

# PyArrow / Polars DataFrame: zero-copy
model = nextstat.from_arrow(table, poi="mu", observations={"SR": [10, 20]})

# Export to Arrow
yields_table = nextstat.to_arrow(model, params=result.bestfit, what="yields")
params_table = nextstat.to_arrow(model, params=result.bestfit, what="parameters")

# Read ROOT histogram (no ROOT C++ dependency)
hist = nextstat.read_root_histogram("file.root", "dir/hist")
# hist["bin_edges"], hist["bin_content"], hist["sumw2"]

# HistFactory XML import (combination.xml + ROOT histograms)
model = nextstat.from_histfactory_xml("config/combination.xml")
edges = nextstat.histfactory_bin_edges_by_channel("config/combination.xml")

# Toy data generators
data = nextstat.data.gaussian(n=1000, mu=0.0, sigma=1.0)
data = nextstat.data.churn(n_customers=5000, seed=42)
```

## Python API — Group Sequential Testing

```python
import nextstat

# O'Brien-Fleming spending function
boundaries = nextstat.group_sequential(
    n_analyses=5, alpha=0.025, spending="obrien_fleming")
# boundaries["critical_values"], boundaries["cumulative_alpha"]

# Pocock spending function
boundaries = nextstat.group_sequential(
    n_analyses=5, alpha=0.025, spending="pocock")

# Lan-DeMets alpha spending
boundaries = nextstat.group_sequential(
    n_analyses=5, alpha=0.025, spending="lan_demets",
    information_fractions=[0.2, 0.4, 0.6, 0.8, 1.0])
```

## Python API — TRExFitter Config Import (nextstat.trex_config)

```python
from nextstat.trex_config import parse_trex_config

# Parse TRExFitter .config → pyhf JSON workspace
workspace = parse_trex_config("myfit.config", histograms_dir="/path/to/hists")
# Supports NTUP and HIST modes, regions, samples, systematics, NFs

# CLI equivalent:
# nextstat import trex-config --input myfit.config --output workspace.json
```

## Python API — Analysis Preprocessing (nextstat.analysis)

```python
from nextstat.analysis.preprocess import smooth, prune, symmetrize, hygiene

# Histogram smoothing (353QH, Rebinning)
smoothed = smooth(histogram, method="353qh")

# Pruning: remove negligible systematics
pruned_ws = prune(workspace, threshold=0.01)

# Symmetrize one-sided systematics
sym_ws = symmetrize(workspace, method="mirror")

# Hygiene checks on workspace
issues = hygiene.check(workspace)
# issues: list of warnings (empty bins, negative yields, etc.)
```

## Python API — Panel Data (nextstat.panel)

```python
import nextstat.panel as panel

# Panel FE with formula interface
result = panel.fe("y ~ x1 + x2 | entity", data, cluster="entity")
# result: PanelFeFit(coef, standard_errors, r_squared, n_obs)
```

## Python API — Robust Statistics (nextstat.robust)

```python
import nextstat.robust as robust

# Huber M-estimator
result = robust.huber(X, y, c=1.345)

# HC0-HC3 heteroskedasticity-consistent standard errors
se = robust.hc_se(X, y, residuals, kind="hc1")
```

## Python API — Summary Statistics (nextstat.summary)

```python
import nextstat.summary as summary

# Descriptive statistics
stats = summary.describe(data)
# stats: {mean, std, min, max, median, q25, q75, n, missing}
```

## Unbinned (Event-Level) Models

NextStat supports unbinned extended maximum likelihood fits via the `ns-unbinned` crate.
Configuration uses `unbinned_spec_v0` JSON schema.

```python
from nextstat import UnbinnedModel

model = UnbinnedModel.from_config("unbinned_spec.yaml")
result = model.fit(data="events.parquet")
scan = model.profile_scan("mu", bounds=(0.0, 3.0), n_points=40)
limit = model.upper_limit("mu", cl=0.95)
```

Supported PDFs (15): gaussian, crystal_ball, double_crystal_ball, exponential, chebyshev,
bernstein, voigtian, argus, spline, histogram, histogram_from_tree, kde, kde_from_tree,
flow (normalizing flow via ONNX), conditional_flow, dcr_surrogate, product.

Neural PDFs: Train flows in PyTorch (zuko), JAX (distrax), or TF → export ONNX → NextStat uses
as full PDF with gradients and normalization. DCR surrogate replaces template morphing.

Features: columnar EventStore (SoA), morphing with systematics, extended likelihood with yield
modifiers, expression engine for selections, GPU-native L-BFGS (44-80× vs lockstep).

## IDE Support — TypedDict Return Types & Type Stubs

NextStat ships `_core.pyi` (2900+ lines) with full type annotations for IDE autocompletion:

- **TypedDict returns**: `SamplerResult`, `SampleStats`, `Diagnostics`, `QualitySummary`,
  `RankingEntry`, `HypotestResult`, `HypotestToysMetaResult`, `ProfileScanPoint`,
  `ProfileScanResult`, `ClsCurveResult`, `PanelFeResult`, `KalmanFilterResult`,
  `MetaAnalysisResult`, `FitResult`, etc. (~25 structured TypedDict definitions)
- **Model classes**: full `__init__` signatures for all 20+ model types
- **Function signatures**: `fit()`, `hypotest()`, `profile_scan()`, `upper_limit()`,
  `ranking()`, `sample()`, `cls_curve()`, all with typed parameters and returns
- **Unified API**: all functions accept `device="cpu"|"cuda"|"metal"` with runtime dispatch
  on `HistFactoryModel` vs `UnbinnedModel` automatically

IDE support: PyCharm, VS Code (Pylance), Cursor, Windsurf, vim (pyright) — full IntelliSense.

## Unified Python API — Device Dispatch

All major inference functions accept `device=` for transparent GPU acceleration:

```python
import nextstat

model = nextstat.HistFactoryModel.from_pyhf("workspace.json")

# CPU (default)
result = nextstat.fit(model)
ranking = nextstat.ranking(model)

# CUDA (NVIDIA GPU, f64)
result = nextstat.fit(model, device="cuda")
ranking = nextstat.ranking(model, device="cuda")

# Metal (Apple Silicon, f32)
result = nextstat.fit(model, device="metal")

# Dispatch is automatic — same function for HistFactoryModel and UnbinnedModel
# Old unbinned_*, *_gpu, *_batch_gpu variants removed in 0.9.6
```

Functions with device dispatch: `fit()`, `fit_toys()`, `hypotest()`, `hypotest_toys()`,
`profile_scan()`, `upper_limit()`, `ranking()`, `sample()` (method="laps" for GPU MCMC).

## CLI Reference

```bash
# HEP / HistFactory
nextstat fit --input workspace.json [--poi mu] [--gpu cuda]
nextstat hypotest --input workspace.json --poi-value 1.0
nextstat upper-limit --input workspace.json
nextstat scan --input workspace.json --poi-range 0:3:30
nextstat ranking --input workspace.json
nextstat report --input workspace.json --output report/
nextstat audit --input workspace.json

# Import / Export
nextstat import histfactory --input config/combination.xml --output workspace.json
nextstat import trex-config --input myfit.config --output workspace.json
nextstat export histfactory --input workspace.json --out-dir export/

# Unbinned (event-level) fits
nextstat unbinned-fit --config unbinned.json --threads 0
nextstat unbinned-scan --config unbinned.json --start 0 --stop 5 --points 21
nextstat unbinned-fit-toys --config unbinned.json --n-toys 100 --seed 42

# Validation report (JSON + optional PDF)
nextstat validation-report \
  --apex2 tmp/apex2_master_report.json \
  --workspace workspace.json \
  --out validation_report.json \
  --pdf validation_report.pdf \
  --deterministic

# Inference server
nextstat-server --port 3742 --gpu cuda
```

## Benchmarks & Validation

NextStat's public benchmark program treats performance as a scientific claim:
- 6 benchmark suites: HEP, Pharma, Bayesian, ML, Time Series, Econometrics
- Correctness gates before timing (parity checks within explicit tolerance)
- Pinned environments (rust-toolchain.toml, Cargo.lock, Python lockfile)
- Published artifacts: raw results, summaries, baseline manifests, validation reports
- DOI snapshots via Zenodo with CITATION.cff
- Third-party replication with GPG/Sigstore signed reports

Validation Report (`nextstat validation-report`) produces:
- `validation_report.json`: dataset SHA-256, model spec, environment, per-suite pass/fail
- `validation_report.pdf`: 7-page audit-ready PDF (requires matplotlib)
- Schema: `validation_report_v1` (via `nextstat config schema --name validation_report_v1`)

Target consumers: HEP (analysis preservation), Pharma (IQ/OQ/PQ, 21 CFR Part 11), FinTech (SR 11-7).

Docs: https://nextstat.io/docs/public-benchmarks
Docs: https://nextstat.io/docs/benchmark-results
Docs: https://nextstat.io/docs/snapshot-registry
Docs: https://nextstat.io/docs/validation-report

## Key technical details

### Interpolation codes
- NormSys (OverallSys): Code4 (polynomial, default) or Code1 (exponential, HistFactory compat)
- HistoSys: Code4p (piecewise polynomial, default) or Code0 (piecewise linear, HistFactory compat)

### pyhf compatibility
- 7-tier tolerance contract from 1e-12 (per-bin likelihood) to 0.05 (toy statistics)
- Parity mode (--parity): Kahan summation, deterministic order, single-threaded
- Default "fast" mode: FMA, parallel, SIMD — may differ by ~1e-9

### Supported formats
- pyhf JSON workspace (auto-detected)
- HS3 v0.2 JSON (ROOT 6.37+ compatible, auto-detected)
- HistFactory XML + ROOT histograms (via from_histfactory / import histfactory)
- TRExFitter .config (via import trex-config, NTUP and HIST modes)
- Arrow IPC / Parquet (via from_arrow / from_parquet)

### GPU support
- CUDA: batch toy fitting, GPU-native L-BFGS, differentiable NLL, profile scans
- Metal: L-BFGS-B fits and batch toys on Apple Silicon (unified memory)
- No GPU required — CPU fallback for all operations

## R Bindings (extendr)

```r
library(nextstat)
model <- ns_from_pyhf("workspace.json")
result <- ns_fit(model)
hypo <- ns_hypotest(model, poi_value = 1.0)
ul <- ns_upper_limit(model)
scan <- ns_scan(model, poi_values = seq(0, 3, length.out = 30))
ranking <- ns_ranking(model)
```

Docs: https://nextstat.io/docs/r-bindings

## Project info

- Version: 0.10.1
- Repository: https://github.com/NextStat/nextstat.io
- License: AGPL-3.0-or-later OR LicenseRef-Commercial
- Install: pip install nextstat
- Language: Rust + Python (PyO3)
- Website: https://nextstat.io
- Playground HEP: https://nextstat.io/playground (Brazil band, profile scan, MLE fit, hypo test, GLM — WASM)
- Playground Bayes: https://nextstat.io/playground/bayes (NUTS/MAMS sampling — Eight Schools, GLM, std normal)
- Playground Pop PK: https://nextstat.io/playground/pk (1-cpt oral PK tutorial — synthetic data + FOCE/SAEM estimation)

## Persona quick-start

### For Data Scientists (sklearn users)
Use `nextstat.glm` for regression, `nextstat.bayes.sample()` for Bayesian posteriors,
`nextstat.formula` for R-style formulas, `nextstat.missing` for explicit missing-data
handling. Arrow/Polars for zero-copy ingest. `nextstat.sklearn` for sklearn-compatible wrappers.

### For Physicists (ROOT/pyhf users)
Use `nextstat.HistFactoryModel.from_pyhf()` for workspace loading, `nextstat.fit/hypotest/
upper_limit/ranking/scan` for the full HEP workflow, `nextstat.unbinned` for event-level
analysis, and `nextstat-server` for shared GPU inference.

### For Quants & Risk Analysts
Use `nextstat.econometrics` for panel FE, DiD, IV/2SLS. `nextstat.KalmanModel` for
state-space models. `nextstat.GevModel`/`nextstat.GpdModel` for EVT. `nextstat.panel`
for fixed-effects estimators. `nextstat.causal.aipw` for doubly-robust treatment effects.

### For Biologists & Pharma Scientists
Use `nextstat.CoxPhModel`/`nextstat.WeibullSurvivalModel` for survival analysis,
`nextstat.kaplan_meier` for KM curves, `nextstat.OneCompartmentOralPkNlmeModel` for
PK/PD, `nextstat.bayes.sample()` for MCMC posteriors. Validation report for GxP compliance.

### For Insurance Actuaries
Use `nextstat.chain_ladder`/`nextstat.mack_chain_ladder` for reserving,
`nextstat.GammaRegressionModel`/`nextstat.TweedieRegressionModel` for GLMs,
`nextstat.churn_*` functions for retention modeling.

## Glossary (cross-domain terminology)

- POI (Parameter of Interest) = Target parameter = Effect size
- Nuisance parameter = Confounding variable = Systematic uncertainty
- NLL = Negative log-likelihood = Loss function (in ML terms)
- CLs = p-value-based exclusion = Hypothesis test
- Profile likelihood = Concentrated likelihood = Marginal MLE
- Ranking = Feature importance = Systematic impact
- Toy = Bootstrap replicate = Pseudo-experiment
- NUTS = No-U-Turn Sampler (adaptive HMC)
- AIPW = Augmented Inverse Probability Weighting
- DiD = Difference-in-Differences
- IV/2SLS = Instrumental Variables / Two-Stage Least Squares
- PPC = Posterior Predictive Check
- KM = Kaplan-Meier
- AFT = Accelerated Failure Time
- GEV = Generalized Extreme Value
- GPD = Generalized Pareto Distribution

## New in 0.10.1

- CUDA WALNUTS sampling — GPU leapfrog via CudaWalnutsPotential trait
- HEP stable-surface matrix — 141-entry surface with validation bundle
- HistFactory + HEPData stable-surface gates
- Release gates: local wheelhouse install, no PyPI during pre-release
- Apex2 prerelease gate governance/performance split
- pharma_fit structural-only comparison for SAEM outputs

## New in 0.10.0

- ICH M15 stable reporting surface
- Ads-native conversion modeling surface
- GVM stable-first verification lane
- Bayesian backend provisioning hardened
- MAMS stable CPU defaults retuned
- Host-backed runners: local wheels without PyPI deps

## New in 0.9.9

- TRExFitter template morphing (ghost samples) — GHOST Type with Template: POI:value, Lagrange polynomial signal templates
- viz PNG/PDF rendering fix — png/pdf features were silently broken since 0.9.8

## New in 0.9.8

- ns-cli-py standalone CLI Python package
- nextstat-nlp remote workflow runner for reproducible NLP pipeline verification
- Validation-pack pharma determinism fix
- pip extras: added missing [all] extra
- Benchmark artifact sanitization (replaced internal hostnames/paths)

## New in 0.9.7

- PD models Python API (Emax, Sigmoid-Emax, IDR I–IV)
- FO / ITS / IMP estimation methods for population PK
- SAEM covariance step (sandwich SE, RSE%, condition number)
- NPDE diagnostics
- 3-compartment PK models (IV + oral) with SAEM
- Per-subject dosing for FOCE/SAEM/VPC/GoF
- Bootstrap NLME (percentile/BCa CI)
- SCM (Stepwise Covariate Modeling) Python API
- CDISC .xpt reader/writer
- NONMEM parity whitepaper + benchmark suite
- IQ/OQ/PQ validation protocol v2 (85 test cases)
- Bioequivalence (TOST, RSABE)
- MAP estimation, Monte Carlo trial simulation, dose optimization
- ODE-based PK (transit, Michaelis-Menten, TMDD)
- nextstat-nlp 0.2.0–0.2.1 (GLiNER2 clinical text extraction)
- EGARCH/GJR-GARCH Python API
- LAPS GLM expansion (linear, Poisson, NegBin, composed logistic)
- NUTS Pathfinder dense metric init
- NUTS vs CmdStan benchmark

## New in 0.9.6

- Pure Rust visualization engine (ns-viz-render) — 17 plot types
- Unified Python API with runtime dispatch
- TypedDict return types (~25 structured defs)
- LAPS Metal backend (Apple Silicon M1–M5)
- BCa confidence interval engine
- Competing risks (Aalen-Johansen, Gray, Fine-Gray)
- EGARCH(1,1) and GJR-GARCH(1,1)
- Group sequential testing
- Interval-censored Weibull AFT
- Cabinetry config reader

## New in 0.9.5

- PyPI wheel coverage fix
- manylinux_2_17 compatibility

## New in 0.9.3

- NUTS Progressive Sampling: https://nextstat.io/blog/nuts-progressive-sampling

## New in 0.9.2

- Parity Contract: https://nextstat.io/docs/parity-contract
- ROOT 3-Way Comparison: https://nextstat.io/docs/root-comparison
- Optimizer Convergence: https://nextstat.io/docs/optimizer
- Physics Assistant Demo: https://nextstat.io/docs/physics-assistant
- White Paper: https://nextstat.io/docs/whitepaper
- R Bindings: https://nextstat.io/docs/r-bindings
- Unbinned Event-Level Analysis: https://nextstat.io/blog/unbinned-event-level-analysis
- Compiler vs Hybrid GPU Fits: https://nextstat.io/blog/compiler-vs-hybrid-gpu-fits