Кому будет полезно
Если вы живёте в Python и одновременно используете
statsmodels,lifelines,pyhf,PyMC/BlackJAX,linearmodels(или что-то похожее).Если вам важны воспроизводимость и понятная валидация численных оптимизаций (особенно в HEP).
Если вам интересна архитектура «одно вычислительное ядро → много задач» и практические hot paths (AOT, SIMD, zero-copy).
TL;DR: Если под капотом у многих статистических моделей — log‑likelihood и градиент, то почему вокруг этого столько разрозненных API и рантаймов? NextStat — AOT‑компилированное ядро на Rust с Python‑обёрткой: один вычислительный слой для likelihood‑based задач (в разных доменах). Мы разделили inference‑логику и численный backend, добавили детерминированный режим (parity) и контракт толерансов для валидации. На наших бенчмарках ускорения получаются от ~4× до сотен раз в конкретных сценариях (ниже — протоколы и оговорки).
Немного честности про науку
В физике высоких энергий работают очень очень умные люди. Математика — на уровне. Статистика — на уровне. Теория — космос. Но если говорить откровенно — системная инженерия там чаще бывает не в приоритете. Не потому что никто не умеет. А потому что у людей другие приоритеты.
Когда ты PhD, ты думаешь:
о сигналах и фонах,
о систематиках,
о дедлайнах,
о том, как бы защититься.
Ты не думаешь:
А насколько оптимально реализован мой likelihood?
Если код работает — отлично. Если сходится — прекрасно. Если публикация вышла — это победа.
Так принято.
Будучи PhD‑студентом, я никогда не задавалась вопросом, удобен ли инструмент. Есть RooFit? Есть RooStats? Есть TRExFitter? Есть pyhf? Есть PyMC? Значит, так принято. Cross‑check в нескольких фреймворках? Да, нормально.
Фит здесь — в RooFit.
Перепроверка — в pyhf.
Байесовская версия — в PyMC / BlackJAX.
Профиль — отдельным скриптом.
Bootstrap — ещё где‑то.
И это, кстати, правильный научный подход. Независимые реализации — это хорошо.
Но потом ты становишься постдоком. И внезапно понимаешь, что делаешь одну и ту же работу трижды.
Маленькая боль, о которой редко пишут
Вы запускаете фит. Он считается. Вы меняете один nuisance parameter. Запускаете снова. Он считается.
И вот вы стоите с чашкой кофе и думаете: «Ну ладно, всё равно считается…».
Это, конечно, шутка. (Но если честно — не совсем.)
Важно отметить: скорость — не единственный критерий
В научной статистике скорость — не главное. Есть вещи важнее:
воспроизводимость,
проверяемость,
прозрачность,
десятилетия валидации,
открытый код,
отсутствие чёрного ящика.
И я не жалуюсь на существующую систему. Она сделала огромную работу. Но сейчас 2026 год. И стало очевидно: часть этих вещей можно сделать архитектурно чище — не ломая научный подход к проверкам. Так родилась идея NextStat.
Проблема: зоопарк пакетов
Типичный requirements.txt для проекта, где нужна статистика:
statsmodels>=0.14 # GLM, OLS lifelines>=0.28 # Cox PH, Kaplan‑Meier pyhf>=0.7 # HEP HistFactory pymc>=5.0 # Байесовский вывод blackjax>=1.0 # MCMC на JAX linearmodels>=5.0 # Panel FE, IV/2SLS
Шесть пакетов. Шесть разных API. Шесть разных подходов к обработке ошибок. Шесть импортов, которые тянут за собой большой стек зависимостей. В итоге окружение разрастается, а установка иногда превращается в отдельный мини‑проект.
И самое неприятное: каждый из этих пакетов решает свою задачу по‑своему. Хотите логистическую регрессию? statsmodels.api.Logit. Хотите её же в байесовской постановке? Перепишите модель с нуля в PyMC. Хотите Cox PH? Отдельная библиотека, отдельный API, отдельная система типов.
А ведь под капотом у всех — одно и то же: функция правдоподобия + оптимизатор.
Одно ядро - все задачи
NextStat — это:
Rust‑ядро: один оптимизатор (L‑BFGS‑B), один MCMC‑движок (NUTS / MAMS / LAPS), одна система типов
Python‑обёртка:
pip install nextstat, 2–5 строк кода на задачуGPU (где это уместно): один параметр
device="cuda"— и батчевые сценарии могут уйти на ускорительБез JIT: AOT‑компиляция, холодный старт = горячий старт
Идея простая: если у вас есть быстрая реализация nll(params) → float и grad_nll(params) → Vec<float>, то оптимизация, сэмплирование, профилирование, доверительные интервалы — всё это одинаково для любой модели. GLM, Cox PH, HistFactory, PK‑модель — всё это частные случаи одного и того же интерфейса:
trait LogDensityModel { fn nll(&self, params: &[f64]) -> Result<f64>; fn grad_nll(&self, params: &[f64]) -> Result<Vec<f64>>; fn dim(&self) -> usize; }
Вы реализуете этот трейт — и получаете бесплатно: MLE‑фит, NUTS/MAMS‑сэмплинг, профиль‑сканирование, BCa‑бутстреп, GPU‑батч, ranking, диагностику.
Что это дает на практике
Вертикаль 1: Байесовский сэмплинг
Конкурент: BlackJAX / NumPyro (JAX)
Проблема: JAX‑based сэмплеры при первом запуске могут тратить минуты на JIT‑компиляцию XLA‑графа. Это не баг — это архитектурное решение: JAX компилирует Python‑код в оптимизированные kernels на лету. Для research‑ноутбуков, где вы запускаете модель один раз и ждёте — часто терпимо. Для production‑пайплайнов, CI/CD или частых перезапусков — это может стать заметной стоимостью.
NextStat: ядро скомпилировано заранее (AOT). Первый запуск = сотый запуск.
import nextstat as ns model = ns.EightSchoolsModel( [28, 8, -3, 7, -1, 1, 18, 12], [15, 10, 16, 11, 9, 11, 10, 18] ) # CPU, 4 цепи, NUTS — как в Stan result = ns.sample(model, method="nuts", n_samples=2000) # GPU (если доступен), большие батчи цепей result = ns.sample(model, method="laps", device="cuda")
Числа (пример одного стенда: Tesla V100, 4096 цепей, включая компиляцию/инициализацию у JAX):
Модель | NextStat (cold) | BlackJAX (cold) | Ускорение |
StdNormal 10d | 0.48 c | 343 с | 715x |
Eight Schools | 0.49 с | 437 с | 889x |
Neal Funnel 10d | 0.45 c | 408 с | 918x |
GLM logistic | 26.4 с | 545 с | 21x |
Протокол: 3 seeds, медиана; для JAX — block_until_ready. Конфиг BlackJAX: NUTS + dual averaging, target_accept=0.9, L=π√d, 500 warmup + 1000 samples.
Но скорость — не главное. Главное — алгоритмическая эффективность. Сколько полезных независимых сэмплов вы получаете на один вызов градиента?
Модель | ESS/grad (NS) | ESS/grad (BJ) | Ratio |
StdNormal 10d | 0.285 | 0.072 | 4.0x |
Eight Schools | 0.267 | 0.002 | 117x |
GLM logistic | 0.285 | 0.0004 | 748x |
Только модели, где оба сэмплера сходятся (R‑hat < 1.05).
Вертикаль 2: Эконометрика
Конкурент: linearmodels (Python)
Проблема: Panel FE с кластерными стандартными ошибками на 16 000 наблюдениях — 4 секунды в linearmodels. Для одного прогона терпимо. Для бутстрепа на 999 повторений — 66 минут.
NextStat: Rust‑ядро + numpy hot paths (PyO3 Vec<Vec<f64>> десериализация оказалась в 21 раз медленнее, чем numpy — мы это измерили и переписали).
import nextstat as ns # Panel Fixed Effects с формулой — как в R fit = ns.econometrics.panel_fe_from_formula( "wage ~ experience + education", data, entity="firm_id", cluster="entity", ) print(f"β = {fit.coef}") print(f"SE = {fit.standard_errors}") # DiD с wild cluster bootstrap boot = ns.econometrics.did_twfe_wild_cluster_bootstrap( x=None, y=y, treat=treat, post=post, entity=entity, time=time, n_boot=999, weight_dist="webb6", ) print(f"ATT = {boot.att:.3f}, p = {boot.p_value:.4f}")
Числа (EPYC 7502P, 16k obs):
Задача | NextStat | linearmodels | Ускорение |
Panel FE, cluster SE | 0.39 с | 4.2 с | 10.8x |
IV/2SLS, HAC | 0.12 с | 0.51 с | 4.2x |
Вертикаль 3: Survival / Cox PH
Конкурент: lifelines (Python)
import nextstat as ns # Cox PH с Efron ties и робастной ковариацией — 3 строки fit = ns.survival.cox_ph.fit( times, events, x, ties="efron", robust=True, ) print(f"HR = {fit.hazard_ratios()}") print(f"95% CI = {fit.hazard_ratio_confint()}") # Кривые выживания для новых пациентов surv = fit.predict_survival(x_new, times=[12, 24, 36, 60]) # Диагностика: тест пропорциональных рисков ph = ns.survival.cox_ph_ph_test(times, events, x) for row in ph: print(f" feature {row['feature']}: p={row['p']:.4f}")
Те же данные, тот же результат, тот же API — но под капотом Rust L‑BFGS‑B вместо scipy.optimize, и аналитические градиенты вместо конечных разностей.
Вертикаль 4: HEP / HistFactory
Конкурент: pyhf (Python), ROOT (C++)
Это нишевая, но показательная история. HistFactory — стандарт ЦЕРН для статистических моделей в физике частиц. Модели с 50–250 параметрами, сотни бинов гистограмм, сложные систематические неопределённости.
import nextstat as ns # Из pyhf workspace — один вызов model = ns.from_pyhf(workspace_json) # Или напрямую из ROOT/XML (без libroot!) model = ns.from_histfactory_xml("combination.xml") # Fit result = ns.fit(model) print(f"μ = {result.bestfit[0]:.4f} ± {result.uncertainties[0]:.4f}") # Profile scan — сетка значений POI scan = ns.profile_scan( model, poi="mu", scan_values=[0.0, 0.5, 1.0, 1.5, 2.0], )
Ключевое: NextStat читает часть ROOT‑форматов нативно на Rust (без зависимости от ROOT C++ runtime) — в объёме, необходимом для HistFactory‑пайплайна. Полный список поддерживаемых типов и ограничений мы держим в документации. I/O при этом изолирован от inference‑ядра.
Числа (pyhf как reference, проверка в рамках tolerance‑контракта; конфигурация важна):
Модель | NextStat | pyhf | Ускорение |
1-канал, 4 param | 0.08 мс | 3 мс | 37x |
3-канала, 184 param | 1.1 мс | 970 мс | 880x |
Profile scan, 100 pts | 3 с | 120 с | 40x |
Дисклеймер: результаты зависят от оптимизатора и настроек. Мы отдельно проверяем совпадение NLL/expected data в parity‑режиме в пределах контракта и фиксируем протокол в артефактах.
Вертикаль 5: Фарма / PK
Конкурент: NONMEM, nlmixr2
Фармакокинетика — это ODE‑модели с mixed effects. Стандарт индустрии — NONMEM (Fortran, проприетарная лицензия, ~$15k/год). Open‑source альтернатива — nlmixr2 ®.
import nextstat as ns # 2-компартментная PK модель (пероральное введение) model = ns.TwoCompartmentOralPkModel(times, concentrations, doses) # SAEM (Stochastic Approximation EM) result = ns.nlme_saem(model, data, n_iter=300) # Диагностика: GoF + VPC ns.pk_gof(result) ns.pk_vpc(result, n_sim=500)
Аналитические градиенты для 2-компартментных моделей (IV и oral): eigenvalue chain rule вместо finite differences. 21 тест, все проходят.
Вертикаль 6: GLM
import nextstat as ns # Логистическая регрессия — MLE model = ns.LogisticRegressionModel(x, y) result = ns.fit(model) # Та же модель — байесовский вывод (NUTS) posterior = ns.sample(model, method="nuts", n_samples=4000) # Или на GPU (LAPS, 4096 цепей) posterior = ns.sample(model, method="laps", device="cuda")
Один и тот же объект модели — MLE‑фит, NUTS, LAPS. Без переписывания.
Архитектура: почему Rust, а не C++ / Zig / Go
Почему не C++?
ROOT написан на C++. Физики из HEP работают с ним каждый день. 30 лет кодовой базы, ручное управление памятью, header‑файлы, cmake, ABI‑несовместимость. Не хотим повторять.
Почему не Python + C extensions
NumPy/SciPy работают через BLAS/LAPACK — но каждый вызов из Python в C — это GIL + overhead на маршаллинг. Для горячего цикла оптимизатора (тысячи вызовов NLL) это убивает производительность.
Почему Rust?
Zero‑cost abstractions: generics без виртуальных вызовов, инлайнинг через монomorphization
Ownership model: zero‑alloc hot path без GC‑пауз. Наш оптимизатор делает 0 аллокаций в основном цикле (
NllScratch+nll_reuse())SIMD: автовекторизация через
widecrate, Accelerate на macOSPyO3: прямой FFI в Python без swig/cython/ctypes
Cargo: единая система сборки, один
cargo testдля всех 8 crate'ов
Схема
┌─────────────────────────────────────────────────────────┐ │ Python API │ │ ns.fit() · ns.sample() · ns.econometrics · ns.survival │ └───────────────────────┬─────────────────────────────────┘ │ PyO3 (zero-copy) ┌───────────────────────┴─────────────────────────────────┐ │ Rust ядро │ │ ┌──────────┐ ┌──────────┐ ┌────────────┐ │ │ │ ns-core │ │ ns-ad │ │ ns-prob │ │ │ │ модели │ │ автодифф │ │ распред. │ │ │ └──────────┘ └──────────┘ └────────────┘ │ │ ┌──────────┐ ┌──────────┐ ┌────────────┐ │ │ │ns-compute│ │ns-infer │ │ns-translate│ │ │ │CUDA/Metal│ │NUTS/MAMS │ │ROOT→Model │ │ │ └──────────┘ └──────────┘ └────────────┘ │ └─────────────────────────────────────────────────────────┘ │ ┌───────────────┼───────────────┐ ▼ ▼ ▼ ┌─────────┐ ┌──────────┐ ┌──────────┐ │ CPU │ │ CUDA │ │ Metal │ │ SIMD │ │ V100+ │ │ Apple │ │ Rayon │ │ 4096ch │ │ f32 PoC │ └─────────┘ └──────────┘ └──────────┘
Где мы НЕ лучше (честный раздел)
Было бы нечестно не сказать о том, где NextStat пока проигрывает или имеет ограничения.
GPU не всегда побеждает CPU
Для малых моделей (< 150 параметров) CPU на Rayon в 250–640 раз быстрее GPU. Причина: CUDA kernel launch overhead доминирует при малом объёме вычислений. Мы это измерили, приняли и рекомендуем GPU только для:
Больших HistFactory моделей (150+ параметров)
Массового MCMC (4096 цепей параллельно)
Дифференцируемых пайплайнов (PyTorch layer)
Centered Neal Funnel не сходится
MAMS‑сэмплер (Microcanonical Adaptive Monte Carlo) использует фиксированную метрику. На задачах с мультимасштабной геометрией (centered Neal Funnel: σ_v=3, σ_x варьируется от 0.007 до 150) — не сходится. Это фундаментальное ограничение алгоритма, не баг.
Решение: NCP‑параметризация (non‑centered parameterization), которая NextStat автоматически использует по умолчанию.
Apple Metal бесполезен для f64
Apple Silicon не имеет аппаратной поддержки double precision. Metal шейдеры работают в f32 — для HEP‑задач, где нужна точность до 1e-8, это не годится. Metal остаётся как proof‑of‑concept для задач, где f32 достаточно (ML, графика).
Экосистема
NextStat не заменяет ArviZ для визуализации, Stan для автоматического дифференцирования произвольных моделей, или формульный DSL statsmodels. Мы не пытаемся быть «PyMC на Rust». Мы — вычислительное ядро, которое делает одно и делает это быстро.
Молодой проект
Документация и тестовое покрытие ещё растут. Некоторые вертикали (фарма, эконометрика) — baseline‑уровня, не production. Мы открыты к обратной связи.
Как попробовать
pip install nextstat
Пример 1: Байесовский вывод (30 секунд до результата)
import nextstat as ns # Eight Schools — классическая иерархическая модель model = ns.EightSchoolsModel( y=[28, 8, -3, 7, -1, 1, 18, 12], sigma=[15, 10, 16, 11, 9, 11, 10, 18], ) result = ns.sample(model, method="nuts", n_samples=2000, seed=42) # Результат — словарь с posterior, diagnostics, sample_stats diag = result["diagnostics"] print(f"R-hat max: {max(diag['r_hat'].values()):.4f}") print(f"ESS min: {diag['quality']['min_ess_tail']:.0f}")
Пример 2: Cox PH (survival analysis)
import nextstat as ns import numpy as np # Синтетические данные rng = np.random.default_rng(42) n = 500 x = rng.standard_normal((n, 3)) true_beta = [0.5, -0.3, 0.1] hazard = np.exp(x @ true_beta) times = rng.exponential(1 / hazard) events = rng.random(n) < 0.8 # 20% цензурирование fit = ns.survival.cox_ph.fit(times.tolist(), events.tolist(), x.tolist(), robust=True) print(f"Estimated HR: {fit.hazard_ratios()}") print(f"True β: {true_beta}")
Пример 3: HistFactory (физика частиц)
import nextstat as ns # Из pyhf JSON workspace workspace = { "channels": [{"name": "SR", "samples": [...]}], "measurements": [...], "observations": [...], } model = ns.from_pyhf(workspace) result = ns.fit(model) print(f"μ = {result.bestfit[0]:.4f} ± {result.uncertainties[0]:.4f}") print(f"NLL = {result.nll:.6f}")
Бенчмарки: сводная таблица
Вертикаль | Конкурент | Задача | NS | Конкурент | Ускорение |
Bayes MCMC | BlackJAX | Eight Schools, V100, cold | 0.49 с | 437 с | 889x |
Bayes MCMC | BlackJAX | GLM logistic, V100, cold | 26.4 с | 545 с | 21x |
Эконометрика | linearmodels | Panel FE, 16k obs | 0.39 с | 4.2 с | 10.8x |
HEP | pyhf | 184-param fit | 1.1 мс | 970 мс | 880x |
HEP | pyhf | Profile scan, 100 pts | 3 с | 120 с | 40x |
Мы старались сделать бенчмарки воспроизводимыми (протокол, версии, seeds). Код и конфиги планируем опубликовать.
Личное
Когда ты PhD, ты принимаешь инструменты как данность. Так принято. Так делают все. Работает — значит хорошо. Когда ты постдок, ты начинаешь замечать системные издержки. Повторяющийся код. Дублирующиеся реализации. Непонятные тормоза. Когда ты инженер — тебе уже хочется это починить. NextStat — это попытка соединить научный опыт и инженерную дисциплину. Не переписать мир. А сделать вычислительное ядро аккуратным, быстрым и единым. И да — теперь, когда я запускаю фит, я не иду за кофе.
Честно говоря, PDF с графиками в Matplotlib иногда генерируется дольше, чем считается likelihood. Мы в команде и эту часть постепенно ускоряем — просто это уже другая история и другой слой пайплайна.
Хотя иногда за кофе я всё равно иду. Просто теперь — не по техническим причинам.
Что дальше?
NextStat — проект в активной разработке. Roadmap на ближайшие месяцы:
WASM Playground: интерактивные демо прямо в браузере (уже есть PoC)
Расширение survival: стратифицированный Cox, time‑varying coefficients
Фарма production: полный NONMEM‑паритет по FOCE/SAEM
Больше GLM: квантильная регрессия, zero‑inflated модели
Inference server: gRPC/REST API для production‑деплоя
Следующая статья
В следующей части — deep dive в архитектуру zero‑alloc Rust‑ядра: как мы добились 0 аллокаций в горячем цикле оптимизатора, как работает compiled modifier pipeline для HistFactory, и почему наш L‑BFGS‑B находит более глубокий минимум, чем scipy SLSQP на больших моделях.
NextStat — open‑source проект.
