Кому будет полезно

• Если вы живёте в 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):

Протокол: 3 seeds, медиана; для JAX — block_until_ready. Конфиг BlackJAX: NUTS + dual averaging, target_accept=0.9, L=π√d, 500 warmup + 1000 samples.

Но скорость — не главное. Главное — алгоритмическая эффективность. Сколько полезных независимых сэмплов вы получаете на один вызов градиента?

Только модели, где оба сэмплера сходятся (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):

Вертикаль 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‑контракта; конфигурация важна):

Дисклеймер: результаты зависят от оптимизатора и настроек. Мы отдельно проверяем совпадение 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?

1. Zero‑cost abstractions: generics без виртуальных вызовов, инлайнинг через монomorphization

2. Ownership model: zero‑alloc hot path без GC‑пауз. Наш оптимизатор делает 0 аллокаций в основном цикле (NllScratch + nll_reuse())

3. SIMD: автовекторизация через wide crate, Accelerate на macOS

4. PyO3: прямой FFI в Python без swig/cython/ctypes

5. 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}")

---

Бенчмарки: сводная таблица

Мы старались сделать бенчмарки воспроизводимыми (протокол, версии, 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 проект.

GitHub

Документация

PyPI