From eb8178b794d57f8f2037222529b33012bf903354 Mon Sep 17 00:00:00 2001
From: J Xu <jxucoder@gmail.com>
Date: Tue, 18 Aug 2026 01:41:03 -0700
Subject: [PATCH] docs: reposition around distributional regression

The docs led with "Boost every parameter", which is not a statistical term
and read as a slogan rather than a claim. Retitle everything around
distributional regression and varying-coefficient models, the established
names for what NaturalBoost/WeibullAFT and FormulaBoost actually fit, and
state the differentiator plainly: the full K x K natural gradient rather
than a diagonal approximation.

Slim the README to a lean overview that links out for detail, and stop
leading with the NGBoost speedup ratio on every pitch surface, since the
suite is still incomplete: three UCI datasets are unmeasured, and
XGBoostLSS and LightGBMLSS are not compared at all. The measured ratio now
lives only on the benchmarks and GPU pages. Comparison sections describe
what each library targets instead of what it cannot do.

New pages: benchmarks, how-it-works, formulaboost, survival. Em dashes
removed throughout.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md                                  |  22 +-
 CLAUDE.md                                     |  57 +++--
 README.md                                     | 219 +++++++---------
 benchmarks/bench_formula.py                   |   2 +-
 benchmarks/bench_survival.py                  |   2 +-
 docs/api/models.md                            |  28 +++
 docs/api/openboost.md                         |   5 +-
 docs/benchmarks.md                            | 134 ++++++++++
 docs/cookbook/custom-distribution.md          |   6 +-
 docs/cookbook/custom-growth-strategy.md       |  14 +-
 docs/cookbook/custom-loss.md                  |   4 +-
 docs/cookbook/device-loss.md                  |   6 +-
 docs/cookbook/index.md                        |  37 +--
 docs/getting-started/gpu-setup.md             | 102 ++++----
 docs/getting-started/installation.md          |  68 +++--
 docs/getting-started/quickstart.md            | 214 +++++++---------
 docs/index.md                                 | 165 +++++++-----
 docs/migration/from-xgboost.md                | 136 +++++-----
 docs/tutorials/custom-loss.md                 |  12 +-
 docs/tutorials/uncertainty.md                 |   7 +-
 docs/user-guide/evaluation/metrics.md         |   3 +-
 docs/user-guide/formulaboost.md               | 138 ++++++++++
 docs/user-guide/how-it-works.md               | 109 ++++++++
 docs/user-guide/model-persistence.md          |  13 +-
 docs/user-guide/models/dart.md                |   5 +-
 docs/user-guide/models/gam.md                 |  18 +-
 docs/user-guide/models/gradient-boosting.md   |   5 +-
 docs/user-guide/models/linear-leaf.md         |   7 +-
 docs/user-guide/models/multiclass.md          |   4 +-
 .../naturalboost/custom-distributions.md      |   7 +-
 docs/user-guide/naturalboost/distributions.md |   5 +-
 docs/user-guide/naturalboost/overview.md      |  57 +++--
 docs/user-guide/sklearn-integration.md        |   4 +-
 docs/user-guide/survival.md                   | 107 ++++++++
 docs/user-guide/training/callbacks.md         |   5 +
 examples/README.md                            | 235 +++++-------------
 mkdocs.yml                                    |  16 +-
 pyproject.toml                                |  15 +-
 src/openboost/__init__.py                     |  38 +--
 src/openboost/_models/_formula.py             |   8 +-
 src/openboost/_models/_survival.py            |   4 +-
 41 files changed, 1265 insertions(+), 778 deletions(-)
 create mode 100644 docs/benchmarks.md
 create mode 100644 docs/user-guide/formulaboost.md
 create mode 100644 docs/user-guide/how-it-works.md
 create mode 100644 docs/user-guide/survival.md

diff --git a/CHANGELOG.md b/CHANGELOG.md
index b1482a4..1a449e0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,8 +7,28 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+- `FormulaBoost`: boost parameters of a user formula `y = f(θ, x)` with
+  finite-difference Jacobian and damped GGN (`precond='full'|'diag'|'plain'`).
+- `WeibullAFT`: right-censored Weibull AFT that boosts both scale `λ(z)`
+  and shape `k(z)` using expected-Fisher natural gradient.
+- Unified `fit_boosting` trainer and `Objective` protocol shared by
+  NaturalBoost, FormulaBoost, and WeibullAFT.
+- Capability benchmarks `benchmarks/bench_formula.py` and
+  `benchmarks/bench_survival.py`; probabilistic speed/quality suite
+  `benchmarks/bench_probabilistic.py`.
+
 ### Changed
 
+- Documentation and package description framed as distributional
+  regression (GAMLSS / NGBoost) and varying-coefficient models: README,
+  docs home, quickstart, GPU setup, NaturalBoost guide, XGBoost
+  migration, and a new benchmarks page. Headline numbers: 1229× vs
+  NGBoost on A100 at 90K (NLL tied); FormulaBoost ~21× better
+  extrapolation than black-box GBDT; WeibullAFT recovers `k(z)`.
+- Performance documentation now quotes only results backed by committed
+  benchmark runs (Modal A100 speed/capability, CPU UCI quality).
 - Made multi-round `fit_trees_batch` recompute gradients and hessians from targets.
 - Consolidated batch configuration and training state into one canonical module.
 - Added manually triggered CUDA verification on Modal GPUs.
@@ -16,8 +36,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   cannot trace a user NLL, instead of returning placeholder gradients.
 - Unsupported `sample_weight` inputs now raise on CUDA, distributed, and
   multi-GPU training paths instead of being ignored.
-- Performance documentation now quotes only results backed by committed raw
-  benchmark artifacts.
 
 ## [1.0.0rc1] - 2026-01-20
 
diff --git a/CLAUDE.md b/CLAUDE.md
index b8e57b9..83d6ef9 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -4,7 +4,19 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 
 ## Project Overview
 
-OpenBoost is a GPU-native, all-Python gradient boosting library (~20K lines). It uses Numba JIT for CPU kernels and CuPy/numba-cuda for GPU acceleration. Designed as a research-friendly alternative to XGBoost/LightGBM with full Python source.
+OpenBoost is GPU gradient boosting for distributional regression (~20K
+lines of Python). Each parameter of the statistical model is a
+histogram-tree ensemble, trained with the full K x K natural gradient /
+GGN rather than a diagonal approximation. Numba JIT on CPU,
+CuPy/numba-cuda on GPU.
+
+Flagship models: `NaturalBoost` (GAMLSS-style distributions),
+`FormulaBoost` (varying-coefficient formula), `WeibullAFT` (censored
+Weibull). Mean-regression GBDT, GAM, DART, and linear-leaf models share
+the same tree engine.
+
+Not a faster XGBoost clone. XGBoost/LightGBM remain the right tool for
+ordinary MSE/logloss GBDT.
 
 ## Commands
 
@@ -44,31 +56,40 @@ uv build                                        # Build wheel/sdist
 ### Layer Overview
 
 ```
-Models (_models/)  →  Core (_core/)  →  Backends (_backends/)
-    ↓                     ↓                    ↓
-GradientBoosting     fit_tree()          _cpu.py (Numba JIT)
-NaturalBoost         histograms          _cuda.py (CuPy kernels)
-OpenBoostGAM         split finding
-DART, LinearLeaf     growth strategies
+Models (_models/)  →  Trainer (_trainer.py)  →  Core (_core/)  →  Backends
+    ↓                      Objective protocol         ↓              ↓
+NaturalBoost           fit_boosting()            fit_tree()     _cpu.py
+FormulaBoost           _objectives.py            histograms     _cuda.py
+WeibullAFT             Distribution / Formula /  split finding
+GradientBoosting       AFT / Loss                growth
+OpenBoostGAM, DART, LinearLeaf
 ```
 
 ### Data Layer (`_array.py`)
-`BinnedArray` is the fundamental data structure — quantile-bins continuous features into uint8 (max 255 bins). Missing values encode as `MISSING_BIN = 255`. Native categorical feature support with auto-detection of string/object columns. All tree-building operates on binned data.
+`BinnedArray` is the fundamental data structure: it quantile-bins continuous features into uint8 (max 255 bins). Missing values encode as `MISSING_BIN = 255`. Native categorical feature support with auto-detection of string/object columns. All tree-building operates on binned data.
 
 ### Core (`_core/`)
-- **`_tree.py`** — `fit_tree()`, `fit_tree_gpu_native()`, `fit_tree_symmetric()`: the main tree-fitting entry points
-- **`_primitives.py`** — Low-level histogram building, split finding, sample partitioning
-- **`_growth.py`** — Three growth strategies: `LevelWiseGrowth` (XGBoost-style), `LeafWiseGrowth` (LightGBM-style), `SymmetricGrowth` (CatBoost-style)
+- **`_tree.py`**: the main tree-fitting entry points `fit_tree()`, `fit_tree_gpu_native()`, `fit_tree_symmetric()`
+- **`_primitives.py`**: Low-level histogram building, split finding, sample partitioning
+- **`_growth.py`**: Three growth strategies: `LevelWiseGrowth` (XGBoost-style), `LeafWiseGrowth` (LightGBM-style), `SymmetricGrowth` (CatBoost-style)
 
 ### Backend Dispatch (`_backends/`)
 `get_backend()` / `set_backend()` switch between CPU and CUDA implementations. Same interface, different kernels. Control via `OPENBOOST_BACKEND` env var or `set_backend('cuda')`. Use `backend_context('cpu')` context manager for temporary switches.
 
 ### Models (`_models/`)
-- **`_boosting.py`** — `GradientBoosting`, `MultiClassGradientBoosting`: main model classes with full callback/eval_set support
-- **`_sklearn.py`** — sklearn-compatible wrappers (`OpenBoostRegressor`, `OpenBoostClassifier`, `OpenBoostDARTRegressor`, `OpenBoostGAMRegressor`, `OpenBoostDistributionalRegressor`, `OpenBoostLinearLeafRegressor`)
-- **`_distributional.py`** — `NaturalBoost`: distributional GBDT with callbacks/eval_set support
-- **`_dart.py`** — `DART`: dropout boosting with callbacks/eval_set support
-- **`_linear_leaf.py`**, **`_gam.py`** — Specialized model variants (both support callbacks/eval_set/early stopping; GAM also supports pairwise interactions, smoothing, and monotone shape constraints)
+- **`_boosting.py`**: `GradientBoosting`, `MultiClassGradientBoosting`
+- **`_distributional.py`**: `NaturalBoost` / `DistributionalGBDT` (facade over the unified trainer)
+- **`_formula.py`**: `FormulaBoost` with user `f(θ, x)`, FD Jacobian, GGN preconditioner
+- **`_survival.py`**: `WeibullAFT` with censored NLL, expected-Fisher natural gradient
+- **`_sklearn.py`**: sklearn wrappers (`OpenBoostRegressor`, `OpenBoostClassifier`, `OpenBoostDARTRegressor`, `OpenBoostGAMRegressor`, `OpenBoostDistributionalRegressor`, `OpenBoostLinearLeafRegressor`)
+- **`_dart.py`**, **`_linear_leaf.py`**, **`_gam.py`**: DART, linear-leaf, GAM
+
+### Trainer (`_trainer.py`) + objectives (`_objectives.py`)
+Single `fit_boosting(objective, X, y, ...)` loop. Facades supply an
+`Objective` (`DistributionObjective`, `FormulaObjective`,
+`WeibullAFTObjective`, `LossObjective`). GPU path: device-resident raw
+scores when the objective is device-capable; `fit_tree_gpu_native` per
+channel.
 
 ### Persistence (`_persistence.py`)
 `PersistenceMixin` provides `save()`/`load()` on all models. Generic `ob.load(path)` auto-detects model class from saved state.
@@ -85,12 +106,12 @@ DART, LinearLeaf     growth strategies
 ## Key Conventions
 
 - **Python 3.10+** target. Ruff rules: E, F, I, UP, B, SIM (line length 100; E501, E402, F821 ignored).
-- **uv only** for package management — never `pip install` or `conda`.
+- **uv only** for package management: never `pip install` or `conda`.
 - All Numba-jitted functions use `@njit` or `@cuda.jit`. CPU kernels are in `_backends/_cpu.py`, CUDA in `_backends/_cuda.py`.
 - Test environment variable `OPENBOOST_BACKEND=cpu` forces CPU backend in CI.
 - Tests use `pytest-xdist` (`-n auto --dist loadfile`) for parallel execution. Shared fixtures are in `tests/conftest.py` (session-scoped datasets, function-scoped gradients).
 - **GPU-native builder** (`fit_tree_gpu_native`) does not support missing values or categorical features. The training loop in `_boosting.py` auto-falls back to `fit_tree()` with a warning when the data has NaN or categorical columns.
-- **Callbacks**: all models — `GradientBoosting`, `MultiClassGradientBoosting`, `DART`, `NaturalBoost`/`DistributionalGBDT`, `LinearLeafGBDT`, and `OpenBoostGAM` — support `callbacks` and `eval_set` in `fit()`.
+- **Callbacks**: every model (`GradientBoosting`, `MultiClassGradientBoosting`, `DART`, `NaturalBoost`/`DistributionalGBDT`, `FormulaBoost`, `WeibullAFT`, `LinearLeafGBDT`, `OpenBoostGAM`) supports `callbacks` and `eval_set` in `fit()`. FormulaBoost eval tuples are `(X, y, model_input)`; WeibullAFT tuples are `(X, y[, event])`.
 - **`random_state`**: `GradientBoosting` and `DART` accept `random_state` for reproducibility. Sklearn wrappers pass it through. DART also accepts `seed` (alias).
 - **`suggest_params()`**: Returns sklearn-style names by default (`n_estimators`). Pass `style='core'` to get core API names (`n_trees`).
 - **Profiling**: `ProfilingCallback` wraps core primitives with timers. Enable via callback or `OPENBOOST_PROFILE=1` env var. Reports go to `logs/` as JSON.
diff --git a/README.md b/README.md
index 61b6d72..01f7ddb 100644
--- a/README.md
+++ b/README.md
@@ -1,181 +1,128 @@
 # OpenBoost
 
-**The hackable gradient boosting platform — probabilistic predictions, interpretable GAMs, and custom algorithms in readable Python, with CPU and CUDA tree backends.**
+**GPU gradient boosting for distributional regression.**
 
-> **Note:** OpenBoost is in active development. APIs may change between releases. Use at your own risk.
+Every parameter of `F(y | x)` gets its own tree ensemble, updated with the full
+`K×K` natural gradient rather than a diagonal approximation. `FormulaBoost`
+extends the same engine to varying-coefficient formulas `y = f(θ(z), x)`.
 
-## Why OpenBoost?
+> 1.0.0rc1. APIs may still move, so install with `--pre` until 1.0.
 
-For standard GBDT, use XGBoost/LightGBM — they're highly optimized C++.
+## Install
 
-For GBDT **variants** (probabilistic predictions, interpretable GAMs, custom algorithms), OpenBoost provides reusable Python primitives and a CUDA tree-building path:
-
-- **NaturalBoost**: full-distribution prediction with a GPU tree path. On the committed CPU comparison, OpenBoost and NGBoost are comparable (0.8-1.3x wall-clock, quality within ~1%) — see [Benchmarks](#benchmarks)
-- **OpenBoostGAM**: interpretable main effects with an optional GPU training path; use the included harness to measure speed and accuracy on your workload
-- **Your own algorithms**: custom losses, distributions, and tree-growth strategies are registration APIs (`register_loss`, `register_distribution`, `register_growth_strategy`), not C++ forks — see the [cookbook](https://jxucoder.github.io/openboost/cookbook/custom-loss/)
-
-Plus: ~20K lines of readable Python. Modify, extend, and build on — no C++ required.
-
-| | XGBoost / LightGBM | OpenBoost |
-|---|---|---|
-| **Code** | 200K+ lines of C++ | ~20K lines of Python |
-| **GPU** | Added later | Native from day one |
-| **Customize** | Modify C++, recompile | Modify Python, reload |
-
-## What You Can Build
-
-OpenBoost provides primitives (histograms, binning, tree fitting) that you combine into algorithms:
-
-- **Standard GBDT** — drop-in gradient boosting with selectable growth strategies (`growth='levelwise' | 'leafwise' | 'symmetric'`), early stopping, and callbacks
-- **Distributional GBDT** — predict full probability distributions with [NGBoost](https://arxiv.org/abs/1910.03225)-style natural gradient boosting
-- **Interpretable GAMs** — explainable feature effects inspired by [EBM](https://arxiv.org/abs/1909.09223)
-- **DART** — [dropout regularization](https://arxiv.org/abs/1505.01866) for reduced overfitting
-- **Linear-leaf models** — linear models in tree leaves for better extrapolation
-- **Your own algorithms** — custom losses, distributions, or entirely new methods
+```bash
+pip install --pre openboost              # core
+pip install --pre "openboost[cuda]"      # GPU trees
+pip install --pre "openboost[sklearn]"   # sklearn wrappers
+```
 
-The core tree-building paths support CPU and CUDA backends. Some features and
-model stages remain CPU-only or deliberately fall back to CPU; see the model
-guides for those boundaries. All models support `save()`/`load()` persistence,
-and most support callbacks and early stopping.
+Python 3.10+. NVIDIA GPU optional (CUDA 11/12).
 
-## Quick Start
+## Distributional regression
 
-**High-level API:**
+**NaturalBoost** predicts a distribution instead of a point. Every parameter
+of `F(y | x)` gets its own tree ensemble, trained by natural gradient.
 
 ```python
 import openboost as ob
 
-model = ob.GradientBoosting(n_trees=100, max_depth=6, random_state=42)
-model.fit(X_train, y_train,
-          eval_set=[(X_val, y_val)],
-          callbacks=[ob.EarlyStopping(patience=10)])
-predictions = model.predict(X_test)
+model = ob.NaturalBoostNormal(n_trees=500, max_depth=3, learning_rate=0.03)
+model.fit(X_train, y_train)
+
+mean = model.predict(X_test)
+lo, hi = model.predict_interval(X_test, alpha=0.1)   # 90% interval
 ```
 
-**sklearn-compatible:**
+**WeibullAFT** takes the same idea to right-censored survival. Both the scale
+and the shape vary by covariate, so the hazard shape is per row rather than
+one global hyperparameter.
 
 ```python
-from openboost import OpenBoostRegressor
-from sklearn.model_selection import GridSearchCV
+import openboost as ob
 
-# Works with GridSearchCV, Pipeline, cross_val_score, etc.
-model = OpenBoostRegressor(n_estimators=100, random_state=42)
-search = GridSearchCV(model, {"max_depth": [4, 6, 8]}, cv=5)
-search.fit(X_train, y_train)
+model = ob.WeibullAFT(n_trees=300, max_depth=3)
+model.fit(Z_train, time_train, event=observed)   # 1 = event, 0 = censored
 
-# Also available: OpenBoostClassifier, OpenBoostDARTRegressor,
-# OpenBoostGAMRegressor, OpenBoostDistributionalRegressor
+params = model.predict_params(Z_test)            # {scale, shape}
+t_hat = model.predict(Z_test)                    # median time
+s = model.predict_survival(Z_test, t=5.0)        # S(5 | z)
 ```
 
-**Hyperparameter suggestions:**
-
-```python
-# Auto-suggest params based on dataset characteristics
-params = ob.suggest_params(X_train, y_train, task='regression', style='core')
-model = ob.GradientBoosting(**params)
-```
+## Varying-coefficient models
 
-**Low-level API** (full control over the training loop):
+**FormulaBoost** boosts the coefficients of a formula you write. Given
+`y = f(θ, x)`, the trees learn `θ(z)` while `x` enters only through `f`.
 
 ```python
+import numpy as np
 import openboost as ob
 
-X_binned = ob.array(X_train)
-pred = np.zeros(len(y_train), dtype=np.float32)
-
-for round in range(100):
-    grad = 2 * (pred - y_train)  # your gradients
-    hess = np.ones_like(grad) * 2
-    tree = ob.fit_tree(X_binned, grad, hess, max_depth=6)
-    pred += 0.1 * tree(X_binned)
-```
-
-## Installation
-
-```bash
-# Current release candidate (recommended while 1.0 is in prerelease)
-pip install --pre openboost
+def sales(theta, x):
+    a, b = theta
+    return a * x ** (1.0 / (1.0 + np.exp(-b * x)))
 
-# With GPU support
-pip install --pre "openboost[cuda]"
+model = ob.FormulaBoost(
+    formula=sales, n_params=2, links=("log", "identity"),
+    param_names=("a", "b"), precond="full",
+)
+model.fit(Z_train, y_train, model_input=x_train)
 
-# With sklearn integration
-pip install --pre "openboost[sklearn]"
+params = model.predict_params(Z_test)              # per-row a(z), b(z)
+yhat = model.predict(Z_test, model_input=x_new)
 ```
 
-`pip install openboost` without `--pre` installs the older stable release.
-
-## Documentation
+## Mean regression
 
-Full docs, tutorials, and API reference: **[jxucoder.github.io/openboost](https://jxucoder.github.io/openboost)**
+The single-parameter case of the same engine: `GradientBoosting`,
+`OpenBoostGAM`, DART, linear-leaf models, and sklearn wrappers. They are here
+because they share the trainer and the tree code, not because they beat
+XGBoost or LightGBM at plain MSE or logloss. Those are optimized C++ and
+should stay your default for point estimates.
 
-- [Getting Started](https://jxucoder.github.io/openboost/getting-started/installation/)
-- [User Guide](https://jxucoder.github.io/openboost/user-guide/models/gradient-boosting/)
-- [API Reference](https://jxucoder.github.io/openboost/api/openboost/)
-- [Examples](./examples/)
-
-## Benchmarks
+## How it compares
 
-### Committed comparison: NaturalBoost vs NGBoost on CPU
+Each library optimizes for a different target. NGBoost introduced
+natural-gradient distributional boosting and stays close to sklearn on CPU.
+XGBoost is the reference for fast mean regression, and its custom-objective
+API takes a diagonal Hessian, which is a sound trade for that goal but cannot
+represent the off-diagonal coupling between formula parameters; `survival:aft`
+likewise holds the Weibull shape fixed across rows. OpenBoost gives up C++
+speed on plain regression in exchange for the full metric and an open model
+class.
 
-The repository includes one current, auditable third-party comparison:
-`benchmarks/results/ngboost_comparison_20260720.json`. It uses fixed seeds,
-identical boosting budgets, and the same train/test splits.
+|                             | NGBoost                                   | XGBoost                                | OpenBoost                          |
+| --------------------------- | ----------------------------------------- | -------------------------------------- | ---------------------------------- |
+| What varies with covariates | Distribution parameters (fixed catalogue) | The mean, or a diagonal custom objective | Distribution or formula parameters |
+| Metric                      | Natural gradient                          | Diagonal Hessian                       | Fisher / full GGN                  |
+| GPU trees                   | No                                        | Yes                                    | Yes                                |
+| Weibull shape `k(z)`        | n/a                                       | Global hyperparameter                  | Per row                            |
 
-| Dataset | OpenBoost / NGBoost fit time | Result |
-|---|---|---|
-| Synthetic heteroscedastic, 10K | 16.4s / 18.8s (1.15x) | OpenBoost slightly better NLL/CRPS/RMSE |
-| Synthetic heteroscedastic, 50K | 74.1s / 95.3s (1.29x) | NGBoost slightly better NLL/CRPS/RMSE |
-| California Housing, 20.6K | 30.6s / 25.0s (0.82x) | OpenBoost slightly better NLL/CRPS/RMSE |
-
-The honest read is CPU parity: neither implementation wins every dataset, and
-quality is within roughly 1% in this run.
-
-Reproduce it with:
-
-```bash
-OPENBOOST_BACKEND=cpu uv run --with ngboost python benchmarks/bench_ngboost_comparison.py
-```
-
-### GPU benchmark harnesses
-
-GPU comparisons are available in `benchmarks/bench_gpu.py` and
-`benchmarks/compare_gpu.py`. Third-party GPU speedups are intentionally not
-quoted here until the exact raw result artifact and environment metadata are
-committed alongside the claim.
-
-```bash
-# Local CUDA GPU
-uv run python benchmarks/bench_gpu.py --task all --scale medium
+## Benchmarks
 
-# Modal A100
-uv run modal run benchmarks/bench_gpu.py --task all --scale medium
-```
+Early and incomplete, so read them as directional rather than settled.
+NaturalBoost matches NGBoost's NLL on the UCI datasets measured so far and
+trains in seconds on an A100 at sizes where NGBoost, which is CPU-only, takes
+most of an hour. FormulaBoost and WeibullAFT recover parameter surfaces that a
+diagonal-Hessian objective cannot. Three UCI datasets have not been measured,
+and XGBoostLSS and LightGBMLSS are not in the comparison yet.
 
-NaturalBoost's CUDA acceleration applies to histogram-based tree building;
-distribution gradients and Fisher/natural-gradient calculations still run on
-CPU. Benchmark end-to-end fit time, accuracy, and calibration on the workload
-you actually care about.
+Numbers, caveats, and reproduce commands are on the
+[benchmarks page](https://jxucoder.github.io/openboost/benchmarks/).
 
-## Roadmap
+## Documentation
 
-**Train-many optimization**: OpenBoost now has a correctness-first API that shares
-binned data across hyperparameter configurations. The next milestone is fusing
-histogram and split work across configurations on GPU, with the sequential path
-serving as the behavioral reference.
+**[jxucoder.github.io/openboost](https://jxucoder.github.io/openboost)**
 
-## References
+- [Quickstart](https://jxucoder.github.io/openboost/getting-started/quickstart/)
+- [How it works](https://jxucoder.github.io/openboost/user-guide/how-it-works/)
+- [NaturalBoost](https://jxucoder.github.io/openboost/user-guide/naturalboost/overview/)
+- [FormulaBoost](https://jxucoder.github.io/openboost/user-guide/formulaboost/)
+- [Weibull AFT](https://jxucoder.github.io/openboost/user-guide/survival/)
+- [Benchmarks](https://jxucoder.github.io/openboost/benchmarks/)
+- [API reference](https://jxucoder.github.io/openboost/api/openboost/)
 
-OpenBoost implements and builds on ideas from these papers:
 
-- **Gradient Boosting**: Friedman, J. H. (2001). [Greedy Function Approximation: A Gradient Boosting Machine](https://projecteuclid.org/euclid.aos/1013203451). *Annals of Statistics*.
-- **XGBoost**: Chen, T., & Guestrin, C. (2016). [XGBoost: A Scalable Tree Boosting System](https://arxiv.org/abs/1603.02754). *KDD*.
-- **LightGBM**: Ke, G., et al. (2017). [LightGBM: A Highly Efficient Gradient Boosting Decision Tree](https://papers.nips.cc/paper/6907-lightgbm-a-highly-efficient-gradient-boosting-decision-tree). *NeurIPS*.
-- **CatBoost**: Prokhorenkova, L., et al. (2018). [CatBoost: Unbiased Boosting with Categorical Features](https://arxiv.org/abs/1706.09516). *NeurIPS*.
-- **NGBoost**: Duan, T., et al. (2020). [NGBoost: Natural Gradient Boosting for Probabilistic Prediction](https://arxiv.org/abs/1910.03225). *ICML*.
-- **EBM**: Nori, H., et al. (2019). [InterpretML: A Unified Framework for Machine Learning Interpretability](https://arxiv.org/abs/1909.09223).
-- **DART**: Rashmi, K. V., & Gilad-Bachrach, R. (2015). [DART: Dropouts meet Multiple Additive Regression Trees](https://arxiv.org/abs/1505.01866). *AISTATS*.
 
 ## License
 
-Apache 2.0
+Apache 2.0
\ No newline at end of file
diff --git a/benchmarks/bench_formula.py b/benchmarks/bench_formula.py
index 22a5687..43e6147 100644
--- a/benchmarks/bench_formula.py
+++ b/benchmarks/bench_formula.py
@@ -1,4 +1,4 @@
-"""Formula capability benchmark: the yardstick for "Boost every parameter".
+"""Formula capability benchmark: varying-coefficient FormulaBoost vs baselines.
 
 FormulaBoost fits a user formula ``y ~ f(theta(z), x)`` where each parameter
 surface ``theta_k(z)`` is its own boosting ensemble and ``x`` enters only
diff --git a/benchmarks/bench_survival.py b/benchmarks/bench_survival.py
index fa47e2d..1bf8ba8 100644
--- a/benchmarks/bench_survival.py
+++ b/benchmarks/bench_survival.py
@@ -1,6 +1,6 @@
 """Survival capability benchmark: Weibull AFT with covariate-dependent shape.
 
-The second capability proof for "Boost every parameter". Data is Weibull AFT
+Capability proof beyond NGBoost's model class. Data is Weibull AFT
 where BOTH the scale lambda(z) AND the shape k(z) depend on covariates:
 
     T ~ Weibull(scale=lambda(z), shape=k(z)),  right-censored at ~35%.
diff --git a/docs/api/models.md b/docs/api/models.md
index 6a1958b..6800613 100644
--- a/docs/api/models.md
+++ b/docs/api/models.md
@@ -2,6 +2,34 @@
 
 All OpenBoost model classes.
 
+## Distributional and varying-coefficient models
+
+### FormulaBoost
+
+::: openboost.FormulaBoost
+    options:
+      show_root_heading: true
+      members:
+        - __init__
+        - fit
+        - predict
+        - predict_params
+
+### WeibullAFT
+
+::: openboost.WeibullAFT
+    options:
+      show_root_heading: true
+      members:
+        - __init__
+        - fit
+        - predict
+        - predict_params
+        - predict_quantile
+        - predict_median
+        - predict_survival
+        - nll
+
 ## Standard GBDT
 
 ### GradientBoosting
diff --git a/docs/api/openboost.md b/docs/api/openboost.md
index b9f6877..57af52a 100644
--- a/docs/api/openboost.md
+++ b/docs/api/openboost.md
@@ -15,9 +15,10 @@ print(ob.get_backend())  # "cuda" or "cpu"
 X_binned = ob.array(X, n_bins=256)
 
 # Models
-model = ob.GradientBoosting(n_trees=100)
 model = ob.NaturalBoostNormal(n_trees=100)
-model = ob.OpenBoostGAM(n_rounds=500)
+model = ob.FormulaBoost(formula=f, n_params=2, links=("log", "identity"))
+model = ob.WeibullAFT(n_trees=300)
+model = ob.GradientBoosting(n_trees=100)
 ```
 
 ## Data Layer
diff --git a/docs/benchmarks.md b/docs/benchmarks.md
new file mode 100644
index 0000000..35b9b44
--- /dev/null
+++ b/docs/benchmarks.md
@@ -0,0 +1,134 @@
+# Benchmarks
+
+Three gates. All numbers below are from committed Modal runs (A100 for
+speed and capability, CPU for UCI quality). This page is the source of
+truth that the README summarizes.
+
+## Speed: NaturalBoost vs NGBoost
+
+Heteroscedastic Normal, 80 features, 500 trees, learning rate 0.03,
+depth 3. OpenBoost on a Modal A100; NGBoost on CPU (it has no GPU).
+
+| n_train | OpenBoost (A100) | NGBoost (CPU) | speedup | OB NLL | NGB NLL | cov90 |
+|--------:|-----------------:|--------------:|--------:|-------:|--------:|------:|
+| 45,000 | 3.01s | 1413.93s | 470× | 2.107 | 2.105 | 0.884 |
+| 90,000 | 2.21s | 2715.58s | **1229×** | 2.108 | 2.102 | 0.899 |
+| 450,000 | 5.40s | n/a | n/a | 2.100 | n/a | 0.901 |
+| 900,000 | 6.94s | n/a | n/a | 2.104 | n/a | 0.902 |
+
+NGBoost was not run past 100K (45 minutes at 90K). Quality is tied at
+every size both finished. Linear-scale extrapolation of NGBoost to 900K
+is ~hours; we do not quote that as a measured speedup.
+
+**Honest reading.** The 1229× is GPU OpenBoost vs CPU NGBoost, because
+that is the product comparison. On CPU the two libraries are ~parity
+(0.8–1.3× wall-clock, NLL/CRPS/RMSE within ~1%) on the older committed
+CPU comparison (`benchmarks/results/ngboost_comparison_20260720.json`).
+
+```bash
+# A100 speed + quality
+uv run modal run benchmarks/bench_probabilistic.py --suite speed
+uv run modal run benchmarks/bench_probabilistic.py::quality
+
+# CPU-only NGBoost comparison (no GPU)
+OPENBOOST_BACKEND=cpu uv run --with ngboost python benchmarks/bench_ngboost_comparison.py
+```
+
+## Quality: UCI vs NGBoost
+
+NGBoost-paper UCI datasets, 20 paired 80/20 splits, shared 500-tree
+budget + patience-50 early stopping on a common val set. NLL, lower is
+better. `delta = OB − NGB` (negative = OpenBoost better). `p` is a paired
+Wilcoxon.
+
+| dataset | OB NLL | NGB NLL | delta | p | cov90 |
+|---|--:|--:|--:|--:|------:|
+| boston | 2.679 | 2.639 | +0.040 | 0.57 | 0.894 |
+| concrete | 3.128 | 3.135 | −0.007 | 0.60 | 0.848 |
+| energy | 1.694 | 1.701 | −0.007 | 0.09 | 0.911 |
+| kin8nm | −0.430 | −0.400 | −0.031 | **2e-6** | 0.851 |
+| protein | 1.930 | 1.943 | −0.013 | **0.002** | 0.940 |
+| wine | 1.028 | 1.031 | −0.003 | 0.13 | 0.879 |
+| yacht | 0.814 | 0.828 | −0.014 | 0.73 | 0.877 |
+| california | 0.584 | 0.596 | −0.011 | **2e-6** | 0.900 |
+
+Tied-or-better on **8 of the 11** datasets in the suite. Significant wins
+on kin8nm, protein, california; no significant loss. boston +0.04 is not
+significant. 90% coverage lands in 0.85–0.94.
+
+**Coverage of the suite is incomplete.** `naval_propulsion_plant`,
+`power`, and `YearPredictionMSD` are **unmeasured**, not neutral: they
+dropped out during an OpenML outage (naval's `data_id` 44898 is a
+deactivated version; the name endpoint was returning 503). The rerun that
+would fill them in has not been completed, so "tied-or-better" covers 8
+datasets and says nothing about the other 3. california is from a local
+baseline run, because Modal's figshare egress returned 403.
+
+```bash
+uv run modal run benchmarks/bench_probabilistic.py::quality
+```
+
+## Capability: FormulaBoost
+
+Sales curve `y = a(z) * x ** sigmoid(b(z) * x)`. Train `x ∈ [0.25, 2.5]`,
+extrap `x ∈ [3, 5]` against the noiseless true curve. 300 rounds, depth 3,
+lr 0.1. n = 200K:
+
+| | test RMSE | extrap RMSE | corr `b` | fit |
+|---|--:|--:|--:|--:|
+| FormulaBoost `full` | 0.139 | **0.183** | **0.877** | 17.1s |
+| FormulaBoost `diag` | 0.130 | 0.181 | 0.730 | 12.0s |
+| FormulaBoost `plain` | 1.410 | 3.931 | 0.652 | 10.3s |
+| global `(a, b)` | 1.641 | 4.396 | n/a | n/a |
+| black-box GBDT | 0.127 | 3.872 | n/a | 0.7s |
+| XGBoost custom (diag Hess) | 0.129 | 0.203 | 0.599 | 19.9s |
+
+Gates (40K and 200K): extrap vs black-box ≥5× (measured **21×**), `full`
+beats `plain`, `full` beats global, `full` not worse than XGBoost-diag.
+
+**Honest reading.** The claim that lands is extrapolation (the formula
+constrains `x`) and recovery of `b(z)`. Full GGN's edge over diag is
+parameter recovery, not test RMSE (`0.181` vs `0.183`). `plain` diverges, so
+GGN is load-bearing. XGBoost's custom-objective API is diagonal-only, so
+it cannot express the off-diagonal term that buys `b(z)`.
+
+```bash
+uv run modal run benchmarks/bench_formula.py
+```
+
+## Capability: Weibull AFT
+
+Both `λ(z)` and `k(z)` vary with covariates. ~35% right-censoring. 300
+rounds, depth 3. n = 200K:
+
+| | C-index | NLL | cov80 | shape corr | fit |
+|---|--:|--:|--:|--:|--:|
+| OpenBoost `WeibullAFT` | **0.680** | **0.761** | 0.803 | **0.997** | 5.4s |
+| XGBoost `survival:aft` (`extreme`) | 0.672 | 0.830 | 0.852 | n/a (global `k = 1.34`) | 11.7s |
+| global constant | 0.500 | 0.896 | 0.810 | n/a | n/a |
+
+C-index is close (ranking follows scale). The NLL gap and the shape
+correlation are the capability: XGBoost holds Weibull shape as one
+hyperparameter. Coverage of the 80% interval is nearer the nominal 0.80.
+
+```bash
+uv run modal run benchmarks/bench_survival.py
+```
+
+## What we are not claiming
+
+- A 3900× speedup at 1M rows. That is a linear extrapolation of NGBoost,
+  not a measurement.
+- That FormulaBoost `full` wins in-sample RMSE. It does not, vs `diag`.
+- That OpenBoost is a faster drop-in for XGBoost/LightGBM mean regression.
+  It is not; those are optimized C++.
+- PGBM as a product competitor. It is a GPU probabilistic reference, not
+  a gate.
+- Quality parity on the full UCI suite. Three datasets (naval, power,
+  YearPredictionMSD) never produced numbers; the claim is 8 of 11.
+
+## JSON artifacts
+
+Reports write to `benchmarks/results/` (gitignored locally; the numbers on
+this page are transcribed from the Modal runs recorded in
+`tasks/todo.md`). Re-run the commands above to refresh them.
diff --git a/docs/cookbook/custom-distribution.md b/docs/cookbook/custom-distribution.md
index 8f94032..2fa076a 100644
--- a/docs/cookbook/custom-distribution.md
+++ b/docs/cookbook/custom-distribution.md
@@ -14,7 +14,7 @@ models block maxima (floods, peak load, worst-case latency). Its NLL is
 
 $$\mathrm{NLL}(y; \mu, \beta) = \log \beta + z + e^{-z}, \qquad z = \frac{y - \mu}{\beta}$$
 
-You only write the NLL — gradients and hessians are derived automatically
+You only write the NLL; gradients and hessians are derived automatically
 (JAX autodiff when `jax` is installed, otherwise exact-chain-rule numerical
 differentiation; no extra dependency required).
 
@@ -71,7 +71,7 @@ assert np.all(lower <= upper)
 `create_custom_distribution` returns an **instance**, which is convenient for
 one-off use. To make the distribution available by name (in `NaturalBoost`,
 `DistributionalGBDT`, and the `OpenBoostDistributionalRegressor` sklearn
-wrapper), register a **class** that constructs with no arguments — the class is
+wrapper), register a **class** that constructs with no arguments. The class is
 instantiated fresh each time the name is resolved:
 
 ```python
@@ -98,7 +98,7 @@ print("'gumbel' listed:", 'gumbel' in ob.list_distributions())
 
 ## Notes
 
-- **Names are case-insensitive** — they are stored lowercased, matching
+- **Names are case-insensitive**: they are stored lowercased, matching
   `ob.get_distribution`. Duplicate names raise `ValueError` unless you pass
   `override=True`.
 - **Fisher information**: custom distributions use an empirical diagonal
diff --git a/docs/cookbook/custom-growth-strategy.md b/docs/cookbook/custom-growth-strategy.md
index b32b509..af9bcc1 100644
--- a/docs/cookbook/custom-growth-strategy.md
+++ b/docs/cookbook/custom-growth-strategy.md
@@ -9,7 +9,7 @@ next to the built-in `'levelwise'`, `'leafwise'`, and `'symmetric'` strategies.
 ## Complete script
 
 The easiest custom strategy wraps a built-in one and rewrites its
-`GrowthConfig`. Here: **depth-capped random-feature growth** — every tree is
+`GrowthConfig`. Here: **depth-capped random-feature growth**, where every tree is
 level-wise but capped at depth 3 and sees only a random 70% of the features,
 regardless of what the model was configured with. (Think of it as an
 extra-randomized, heavily regularized forest layer.)
@@ -42,8 +42,8 @@ class ShallowRandomGrowth(ob.LevelWiseGrowth):
         )
 
 
-# Register once (process-wide); the class must construct with no arguments —
-# it is instantiated fresh each time the name is resolved.
+# Register once (process-wide); the class must construct with no arguments,
+# because it is instantiated fresh each time the name is resolved.
 ob.register_growth_strategy('shallow_random', ShallowRandomGrowth)
 
 # --- Use it by name ----------------------------------------------------------
@@ -79,13 +79,13 @@ class MyGrowth(ob.GrowthStrategy):
         return tree  # an ob.TreeStructure
 ```
 
-- `binned` — binned features, shape `(n_features, n_samples)`, uint8 (bin 255
+- `binned`: binned features, shape `(n_features, n_samples)`, uint8 (bin 255
   is reserved for missing values).
-- `grad`, `hess` — float32 arrays of shape `(n_samples,)`.
-- `config` — an `ob.GrowthConfig` dataclass carrying `max_depth`,
+- `grad`, `hess`: float32 arrays of shape `(n_samples,)`.
+- `config`: an `ob.GrowthConfig` dataclass carrying `max_depth`,
   `max_leaves`, `min_child_weight`, `reg_lambda`, `reg_alpha`, `min_gain`,
   `subsample`, and `colsample_bytree`.
-- `has_missing` / `is_categorical` / `n_categories` — optional per-feature
+- `has_missing` / `is_categorical` / `n_categories`: optional per-feature
   metadata arrays; a minimal strategy may ignore them (numeric, fully-observed
   data), but then must not be used with missing values or categoricals.
 - Return an `ob.TreeStructure` (routing arrays + leaf values), which handles
diff --git a/docs/cookbook/custom-loss.md b/docs/cookbook/custom-loss.md
index 301249e..6adb969 100644
--- a/docs/cookbook/custom-loss.md
+++ b/docs/cookbook/custom-loss.md
@@ -1,7 +1,7 @@
 # Recipe: Custom Loss with a True Loss Value
 
 Register a custom objective under a string name so it works everywhere a
-built-in loss name does — including correct train/validation loss reporting.
+built-in loss name does, including correct train/validation loss reporting.
 
 **You will use:** `ob.register_loss`, the `loss_value_fn` hook,
 `GradientBoosting(loss='<your name>')`.
@@ -10,7 +10,7 @@ built-in loss name does — including correct train/validation loss reporting.
 
 A callable passed as `loss=` only tells OpenBoost the *gradient* and *hessian*.
 When training needs a scalar loss (for `Logger`, `EarlyStopping`, or history),
-OpenBoost falls back to a second-order Taylor proxy `mean(grad² / (2·hess))` —
+OpenBoost falls back to a second-order Taylor proxy `mean(grad² / (2·hess))`,
 which is usually *not* your actual loss. Registering the loss with a
 `loss_value_fn` fixes that, and gives the loss a reusable name.
 
diff --git a/docs/cookbook/device-loss.md b/docs/cookbook/device-loss.md
index ba86c20..9e382ed 100644
--- a/docs/cookbook/device-loss.md
+++ b/docs/cookbook/device-loss.md
@@ -6,7 +6,7 @@ entirely on the GPU, skipping the per-round host round-trip.
 **You will use:** `ob.device_loss`, `GradientBoosting(loss=<callable>)`.
 
 !!! warning "This recipe needs CUDA to show a benefit"
-    `ob.device_loss` is a **no-op on the CPU backend** — the script below runs
+    `ob.device_loss` is a **no-op on the CPU backend**, so the script below runs
     anywhere (and is what CI runs), but the round-trip it eliminates only
     exists on the CUDA backend. The honest claim is: same results everywhere,
     faster only on GPU.
@@ -23,7 +23,7 @@ device. For big datasets that copy can dominate the round.
 Decorating a loss with `@ob.device_loss` sets `fn.__openboost_device__ = True`
 and changes what the CUDA path hands you:
 
-- `pred` is the **device** (CuPy) prediction array — no copy.
+- `pred` is the **device** (CuPy) prediction array: no copy.
 - `y` is the training target **already resident on the device** (moved once
   per `fit` and cached).
 - You **must return device `(grad, hess)` float32 arrays** with the same
@@ -94,7 +94,7 @@ model.fit(X_train, y_train)  # gradients never leave the GPU
 ## Notes
 
 - **Correctness first**: returning host arrays (or non-float32) from a marked
-  loss on the CUDA path violates the contract — if you cannot keep the math on
+  loss on the CUDA path violates the contract. If you cannot keep the math on
   the device, simply leave the loss unmarked and accept the round-trip.
 - Works with named registration too: `ob.register_loss('gpu_logcosh',
   gpu_logcosh)` keeps the device marker, since the callable itself carries it.
diff --git a/docs/cookbook/index.md b/docs/cookbook/index.md
index 2b34fd2..61aae75 100644
--- a/docs/cookbook/index.md
+++ b/docs/cookbook/index.md
@@ -1,9 +1,17 @@
 # Extending OpenBoost
 
-OpenBoost is all-Python, so every extension point is a plain Python object plus
-a registry call — no C++ plugins, no recompilation. Three registries make
-custom components usable **by name**, exactly like the built-ins, and one
-decorator opts a custom loss into GPU-native execution:
+OpenBoost is all-Python, so every extension point is a plain Python object
+plus a registry call: no C++ plugins, no recompilation.
+
+If the thing you want to boost is a **formula** `y = f(θ, x)`, you do not
+need a registry: use [FormulaBoost](../user-guide/formulaboost.md)
+and pass the callable. If it is a **distribution NLL**, use the distribution
+registry below. If it is **right-censored Weibull**, use
+[WeibullAFT](../user-guide/survival.md).
+
+Three registries make custom components usable **by name**, exactly like
+the built-ins, and one decorator opts a custom loss into GPU-native
+execution:
 
 | Extension point | Register with | Then use as |
 |---|---|---|
@@ -18,8 +26,8 @@ Shared registry rules:
   `register_growth_strategy` take a *class* that must construct with no
   arguments; it is instantiated fresh each time the name is resolved. Names
   are stored lowercased (case-insensitive lookup).
-- **No silent replacement**: registering an existing name — including a
-  built-in like `'mse'` or `'levelwise'` — raises `ValueError` unless you pass
+- **No silent replacement**: registering an existing name, including a
+  built-in like `'mse'` or `'levelwise'`, raises `ValueError` unless you pass
   `override=True`.
 - **Process-wide and import-time**: registrations live for the lifetime of the
   Python process. Put them at import time of your own module so saved models
@@ -31,24 +39,25 @@ Shared registry rules:
 Each recipe is a complete, copy-pasteable script (they are executed in CI, so
 they stay runnable):
 
-- **[Custom Loss](custom-loss.md)** — an asymmetric objective registered by
+- **[Custom Loss](custom-loss.md)**: an asymmetric objective registered by
   name, with a `loss_value_fn` so logging and early stopping report the true
   loss instead of a Taylor proxy.
-- **[Custom Distribution](custom-distribution.md)** — a Gumbel distribution
+- **[Custom Distribution](custom-distribution.md)**: a Gumbel distribution
   for NaturalBoost from just its NLL (autodiff or numerical gradients),
   registered so `distribution='gumbel'` works.
-- **[Custom Growth Strategy](custom-growth-strategy.md)** — depth-capped
+- **[Custom Growth Strategy](custom-growth-strategy.md)**: depth-capped
   random-feature growth in ~20 lines, plus the full `GrowthStrategy` contract
   for from-scratch strategies.
-- **[Device-Native Loss (GPU)](device-loss.md)** — the `@ob.device_loss`
+- **[Device-Native Loss (GPU)](device-loss.md)**: the `@ob.device_loss`
   contract for computing gradients entirely on the GPU (honest note: it needs
   CUDA to show any benefit; on CPU it is a no-op).
 
 ## Related pages
 
-- [Custom Loss Functions tutorial](../tutorials/custom-loss.md) — gradient/
+- [FormulaBoost](../user-guide/formulaboost.md): boost parameters of a formula, no registry required.
+- [Custom Loss Functions tutorial](../tutorials/custom-loss.md): gradient/
   hessian derivations for many classic losses.
-- [Custom Distributions guide](../user-guide/naturalboost/custom-distributions.md)
-  — link functions, `init_fn`, JAX vs numerical gradients.
-- [Callbacks](../user-guide/training/callbacks.md) — the training hooks that
+- [Custom Distributions guide](../user-guide/naturalboost/custom-distributions.md):
+  link functions, `init_fn`, JAX vs numerical gradients.
+- [Callbacks](../user-guide/training/callbacks.md): the training hooks that
   registered losses report into.
diff --git a/docs/getting-started/gpu-setup.md b/docs/getting-started/gpu-setup.md
index b748347..63b04e6 100644
--- a/docs/getting-started/gpu-setup.md
+++ b/docs/getting-started/gpu-setup.md
@@ -1,87 +1,103 @@
 # GPU Setup
 
-OpenBoost automatically detects and uses CUDA GPUs when available.
+OpenBoost uses CUDA for histogram building and tree construction.
+NaturalBoost, FormulaBoost, and WeibullAFT share that tree path. Some
+objective math still runs on the host (FormulaBoost GGN today; LogNormal /
+digamma families). Trees are the expensive part at scale.
 
-## Verify GPU Detection
+## Verify detection
 
 ```python
 import openboost as ob
 
-print(f"Backend: {ob.get_backend()}")  # "cuda" or "cpu"
-print(f"Using GPU: {ob.is_cuda()}")    # True if GPU active
+print(ob.get_backend())   # "cuda" or "cpu"
+print(ob.is_cuda())       # True if a GPU is active
 ```
 
-## Manual Backend Selection
+Install the extra first: `pip install --pre "openboost[cuda]"`.
+
+## Pin a backend
 
 ```python
 import openboost as ob
 
-# Force CPU (useful for debugging or comparison)
-ob.set_backend("cpu")
-
-# Force GPU
+ob.set_backend("cpu")    # debug / comparison
 ob.set_backend("cuda")
 
-# Or use environment variable
+# Or:
 # export OPENBOOST_BACKEND=cuda
 ```
 
-## GPU Performance
+`backend_context("cpu")` is a temporary switch that restores the previous
+backend on exit.
+
+## What the A100 numbers actually are
+
+NaturalBoost vs NGBoost, heteroscedastic Normal, 80 features, 500 trees,
+Modal A100. NGBoost has **no GPU implementation**, so this is GPU OpenBoost
+against CPU NGBoost, which is the comparison that exists in the world.
+
+| n_train | OpenBoost (A100) | NGBoost (CPU) | speedup |
+|--------:|-----------------:|--------------:|--------:|
+| 45K | 3.01s | 1414s | 470× |
+| 90K | 2.21s | 2716s | **1229×** |
+| 450K | 5.40s | skipped (hours) | n/a |
+| 900K | 6.94s | skipped | n/a |
+
+NLL is tied at every size that both ran. Full tables:
+[Benchmarks](../benchmarks.md).
 
-GPU acceleration provides significant speedups for larger datasets:
+!!! tip "When GPU helps"
+    Histogram trees win at tens of thousands of rows and up. Below ~5K
+    samples, kernel launch overhead often matches CPU. Use `float32`
+    features. Missing values and categoricals fall back from the
+    GPU-native builder to the hybrid path (warning emitted).
 
-| Dataset Size | Typical Speedup |
-|--------------|-----------------|
-| <5K samples | ~1x (CPU overhead dominates) |
-| 5K-10K | 2-7x |
-| 25K+ | 2-3x |
-| 100K+ | 5-10x |
+On CPU, NaturalBoost and NGBoost are ~parity (0.8–1.3×). Do not quote the
+A100 ratio as a CPU claim.
 
-!!! tip "Best practices for GPU"
-    - Ensure data is `float32` (not `float64`)
-    - Use larger datasets (GPU overhead not worth it for <5K samples)
-    - GPU shows best speedup at 10K+ samples
+## FormulaBoost and WeibullAFT on GPU
 
-## Multi-GPU Training
+Both use GPU trees. FormulaBoost's GGN (finite-difference Jacobian +
+`K×K` solve) currently runs on the host; at 200K rows full GGN still beat
+an XGBoost custom objective (17s vs 20s on A100). WeibullAFT's expected
+Fisher step is cheap (`2×2` per row).
+
+## Multi-GPU
 
 ```python
 import openboost as ob
 
-# Use multiple GPUs with Ray (requires ray[default])
 model = ob.GradientBoosting(n_trees=100, n_gpus=4)
 model.fit(X, y)
 
-# Or specify exact GPU devices
 model = ob.GradientBoosting(n_trees=100, devices=[0, 2])
 model.fit(X, y)
 ```
 
+Requires `pip install --pre "openboost[distributed]"` (Ray).
+NaturalBoost / FormulaBoost / WeibullAFT currently train on one GPU.
+
 ## Requirements
 
-- NVIDIA GPU with CUDA Compute Capability 3.5+
-- CUDA Toolkit 11.0+ or 12.0+
-- Numba 0.60+ (`pip install numba`)
+- NVIDIA GPU, CUDA Compute Capability 3.5+
+- CUDA Toolkit 11 or 12
+- `numba-cuda>=0.23`
 
 ## Troubleshooting
 
-### Training seems slow on GPU
+**Training seems slow on GPU.** Features should be `float32`. Tiny datasets
+do not amortize kernel launch. Confirm `ob.is_cuda()` is True.
+
+**CUDA not detected.**
 
-- Ensure data is `float32` (not `float64`)
-- Use larger datasets (GPU overhead not worth it for <5K samples)
-- GPU shows best speedup at 10K+ samples
+1. `nvidia-smi`
+2. `python -c "from numba import cuda; print(list(cuda.gpus))"`
+3. Reinstall `openboost[cuda]` against the CUDA version on the machine
 
-### Model trained on GPU, loading on CPU machine
+**Trained on GPU, loading on CPU.** Saved models are backend-agnostic.
 
 ```python
-# Models are saved in a backend-agnostic format
 model.save("model.joblib")
-
-# Load on any machine (CPU or GPU)
-loaded = ob.GradientBoosting.load("model.joblib")
+loaded = ob.NaturalBoostNormal.load("model.joblib")  # CPU or GPU
 ```
-
-### CUDA not detected
-
-1. Check CUDA installation: `nvcc --version`
-2. Check Numba can see GPU: `python -c "from numba import cuda; print(cuda.gpus)"`
-3. Ensure compatible CUDA version with Numba
diff --git a/docs/getting-started/installation.md b/docs/getting-started/installation.md
index 6a81837..c7536eb 100644
--- a/docs/getting-started/installation.md
+++ b/docs/getting-started/installation.md
@@ -1,76 +1,92 @@
 # Installation
 
-## Quick Install
+OpenBoost 1.0 is a release candidate. Install with `--pre` until 1.0.0.
+
+## Quick install
 
 === "pip"
 
     ```bash
-    pip install openboost
+    pip install --pre openboost
     ```
 
 === "uv"
 
     ```bash
-    uv add openboost
+    uv add --prerelease=allow openboost
     ```
 
-=== "conda"
-
-    ```bash
-    # Coming soon
-    conda install -c conda-forge openboost
-    ```
+`pip install openboost` without `--pre` still resolves the older stable
+release.
 
-## With GPU Support
+## GPU support
 
-For CUDA GPU acceleration:
+Numba CUDA kernels for histogram trees. Requires an NVIDIA GPU.
 
 === "pip"
 
     ```bash
-    pip install "openboost[cuda]"
+    pip install --pre "openboost[cuda]"
     ```
 
 === "uv"
 
     ```bash
-    uv add "openboost[cuda]"
+    uv add --prerelease=allow "openboost[cuda]"
     ```
 
-## Optional Dependencies
+Then:
+
+```python
+import openboost as ob
+print(ob.get_backend(), ob.is_cuda())   # "cuda" True when a GPU is visible
+```
+
+See [GPU setup](gpu-setup.md) for backend pinning, multi-GPU, and
+troubleshooting.
+
+## Optional extras
 
 | Extra | What it includes | Install |
 |-------|-----------------|---------|
-| `cuda` | CuPy for GPU acceleration | `pip install "openboost[cuda]"` |
-| `sklearn` | scikit-learn integration | `pip install "openboost[sklearn]"` |
-| `distributed` | Ray for multi-GPU training | `pip install "openboost[distributed]"` |
-| `all` | Everything | `pip install "openboost[all]"` |
+| `cuda` | numba-cuda + CuPy for GPU trees | `pip install --pre "openboost[cuda]"` |
+| `sklearn` | scikit-learn wrappers | `pip install --pre "openboost[sklearn]"` |
+| `jax` | autodiff for custom distributions / formulas | `pip install --pre "openboost[jax]"` |
+| `distributed` | Ray for multi-GPU | `pip install --pre "openboost[distributed]"` |
+| `all` | Everything | `pip install --pre "openboost[all]"` |
+
+Finite-difference Jacobians work without JAX. Install `jax` when you want
+autodiff on a custom NLL.
 
 ## Requirements
 
 - Python 3.10+
 - NumPy 1.24+
 - Numba 0.60+
+- SciPy 1.10+
 
-### For GPU Support
+### For GPU
 
-- NVIDIA GPU with CUDA Compute Capability 3.5+
-- CUDA Toolkit 11.0+ or 12.0+
+- NVIDIA GPU, CUDA Compute Capability 3.5+
+- CUDA Toolkit 11 or 12
+- `numba-cuda>=0.23`, `cupy-cuda12x>=13`
 
-## Verify Installation
+## Verify
 
 ```python
 import openboost as ob
 
-print(f"OpenBoost version: {ob.__version__}")
-print(f"Backend: {ob.get_backend()}")  # "cuda" or "cpu"
-print(f"GPU available: {ob.is_cuda()}")
+print(f"OpenBoost {ob.__version__}")
+print(f"Backend: {ob.get_backend()}")   # "cuda" or "cpu"
+print(f"GPU:     {ob.is_cuda()}")
 ```
 
-## Development Installation
+## Development install
 
 ```bash
 git clone https://github.com/jxucoder/openboost.git
 cd openboost
 uv sync --extra dev
+# GPU kernels:
+uv sync --extra cuda --extra dev
 ```
diff --git a/docs/getting-started/quickstart.md b/docs/getting-started/quickstart.md
index 7ef7627..2a6deba 100644
--- a/docs/getting-started/quickstart.md
+++ b/docs/getting-started/quickstart.md
@@ -1,166 +1,130 @@
 # Quickstart
 
-Get up and running with OpenBoost in 5 minutes.
+Distributional regression and a varying-coefficient formula in a few
+minutes. Each example is self-contained.
 
-## Basic Regression
+## NaturalBoost: a distribution, not a point
 
 ```python
 import numpy as np
 import openboost as ob
 
-# Generate sample data
-np.random.seed(42)
-X = np.random.randn(1000, 10).astype(np.float32)
-y = (X[:, 0] * 2 + X[:, 1] + np.random.randn(1000) * 0.1).astype(np.float32)
-
-# Split data
-X_train, X_test = X[:800], X[800:]
-y_train, y_test = y[:800], y[800:]
-
-# Train model
-model = ob.GradientBoosting(
-    n_trees=100,
-    max_depth=6,
-    learning_rate=0.1,
-    loss='mse',
-)
-model.fit(X_train, y_train)
+rng = np.random.default_rng(0)
+X = rng.standard_normal((2000, 8)).astype(np.float32)
+noise = 0.4 + np.abs(X[:, 0])
+y = (2.0 * X[:, 0] + X[:, 1] + noise * rng.standard_normal(2000)).astype(np.float32)
+X_train, X_test = X[:1600], X[1600:]
+y_train, y_test = y[:1600], y[1600:]
 
-# Predict
-predictions = model.predict(X_test)
+model = ob.NaturalBoostNormal(n_trees=200, max_depth=3, learning_rate=0.05)
+model.fit(X_train, y_train)
 
-# Evaluate
-rmse = np.sqrt(np.mean((predictions - y_test) ** 2))
-print(f"RMSE: {rmse:.4f}")
+mean = model.predict(X_test)
+lo, hi = model.predict_interval(X_test, alpha=0.1)
+print(f"90% coverage: {np.mean((y_test >= lo) & (y_test <= hi)):.1%}")
 ```
 
-## Binary Classification
-
-```python
-import openboost as ob
+Pick a family that matches the data (`NaturalBoostLogNormal`, `Gamma`,
+`Poisson`, `Tweedie`, `NegBin`, `StudentT`). See
+[distributions](../user-guide/naturalboost/distributions.md).
 
-model = ob.GradientBoosting(
-    n_trees=100,
-    max_depth=6,
-    loss='logloss',
-)
-model.fit(X_train, y_train)  # y_train: 0 or 1
+## FormulaBoost: boost the parameters of a formula
 
-# Get probabilities
-logits = model.predict(X_test)
-probabilities = 1 / (1 + np.exp(-logits))
-
-# Get class predictions
-predictions = (probabilities > 0.5).astype(int)
-```
-
-## Multi-Class Classification
+Features `Z` learn parameter surfaces; a structural input `x` (spend, dose,
+time) goes into a formula you write.
 
 ```python
+import numpy as np
 import openboost as ob
 
-y_class = np.random.randint(0, 5, size=len(X_train))  # labels 0..4
-
-model = ob.MultiClassGradientBoosting(
-    n_classes=5,
-    n_trees=100,
-    max_depth=6,
+def sales(theta, x):
+    a, b = theta
+    return a * x ** (1.0 / (1.0 + np.exp(-np.clip(b * x, -30.0, 30.0))))
+
+rng = np.random.default_rng(0)
+n = 4000
+Z = rng.uniform(0, 1, (n, 4)).astype(np.float32)
+x = rng.uniform(0.3, 2.5, n)
+a = np.exp(0.4 + 0.8 * Z[:, 0] - 0.5 * Z[:, 1])
+b = 0.6 + 1.5 * Z[:, 2]
+y = sales((a, b), x) + 0.05 * rng.standard_normal(n)
+
+model = ob.FormulaBoost(
+    formula=sales, n_params=2, links=("log", "identity"),
+    param_names=("a", "b"), precond="full",
+    n_trees=80, max_depth=3, learning_rate=0.1,
 )
-model.fit(X_train, y_class)  # labels must be 0, 1, 2, 3, or 4
-
-# Get probabilities
-probabilities = model.predict_proba(X_test)  # Shape: (n_samples, n_classes)
-
-# Get class predictions
-predictions = model.predict(X_test)
+model.fit(Z[:3200], y[:3200], model_input=x[:3200])
+params = model.predict_params(Z[3200:])
+print(params["a"][:5], params["b"][:5])
 ```
 
-## Uncertainty Quantification
+`precond="full"` (default) is the GGN preconditioner. `plain` (raw
+gradient) diverges on this formula. Details:
+[FormulaBoost](../user-guide/formulaboost.md).
+
+## WeibullAFT: survival with a per-row shape
 
 ```python
+import numpy as np
 import openboost as ob
 
-# Train probabilistic model
-model = ob.NaturalBoostNormal(n_trees=100, max_depth=4)
-model.fit(X_train, y_train)
-
-# Point prediction (mean)
-mean = model.predict(X_test)
-
-# 90% prediction interval
-lower, upper = model.predict_interval(X_test, alpha=0.1)
-
-# Sample from predicted distribution
-samples = model.sample(X_test, n_samples=100)
+rng = np.random.default_rng(0)
+n = 3000
+Z = rng.standard_normal((n, 6)).astype(np.float32)
+lam = np.exp(0.5 + 0.8 * Z[:, 0])
+k = np.exp(0.2 + 0.6 * Z[:, 1])
+u = rng.random(n)
+time = lam * ((-np.log(u)) ** (1.0 / k))
+event = (rng.random(n) > 0.3).astype(np.float64)  # ~30% right-censored
+
+model = ob.WeibullAFT(n_trees=150, max_depth=3, learning_rate=0.1)
+model.fit(Z[:2400], time[:2400], event=event[:2400])
+params = model.predict_params(Z[2400:])
+t_hat = model.predict(Z[2400:])            # median time
+s = model.predict_survival(Z[2400:], t=1.0)
+print(params["scale"][:3], params["shape"][:3], t_hat[:3])
 ```
 
-## sklearn-Compatible API
-
-```python
-from openboost import OpenBoostRegressor, OpenBoostClassifier
-from sklearn.model_selection import cross_val_score, GridSearchCV
-
-# Regressor
-reg = OpenBoostRegressor(n_estimators=100, max_depth=6)
-reg.fit(X_train, y_train)
-print(f"R² Score: {reg.score(X_test, y_test):.4f}")
-
-# Cross-validation
-scores = cross_val_score(reg, X, y, cv=5)
-print(f"CV Score: {scores.mean():.4f} ± {scores.std():.4f}")
-
-# Grid search
-param_grid = {
-    'n_estimators': [50, 100, 200],
-    'max_depth': [4, 6, 8],
-    'learning_rate': [0.05, 0.1, 0.2],
-}
-grid = GridSearchCV(reg, param_grid, cv=3)
-grid.fit(X_train, y_train)
-print(f"Best params: {grid.best_params_}")
-```
+XGBoost's `survival:aft` cannot vary Weibull shape with covariates.
+[Weibull AFT](../user-guide/survival.md).
 
-## Callbacks
+## Point-estimate GBDT (also here)
 
 ```python
-import openboost as ob
-from openboost import EarlyStopping, Logger
+model = ob.GradientBoosting(n_trees=100, max_depth=6, loss="mse")
+model.fit(X_train, y_train)
+pred = model.predict(X_test)
+```
 
-model = ob.GradientBoosting(n_trees=500, max_depth=6)
+Binary classification uses `loss="logloss"`. Multi-class uses
+`ob.MultiClassGradientBoosting`. sklearn wrappers
+(`OpenBoostRegressor`, `OpenBoostClassifier`,
+`OpenBoostDistributionalRegressor`) work with `GridSearchCV` and
+`Pipeline`. See [Gradient Boosting](../user-guide/models/gradient-boosting.md)
+and [sklearn integration](../user-guide/sklearn-integration.md).
 
-callbacks = [
-    EarlyStopping(patience=10, min_delta=0.001),
-    Logger(period=10),
-]
+## Callbacks, GPU, persistence
 
+```python
+model = ob.NaturalBoostNormal(n_trees=500, max_depth=3)
 model.fit(
     X_train, y_train,
-    callbacks=callbacks,
     eval_set=[(X_test, y_test)],
+    callbacks=[ob.EarlyStopping(patience=20), ob.Logger(period=20)],
 )
-
-print(f"Stopped at {len(model.trees_)} trees")
+model.save("model.joblib")
+loaded = ob.NaturalBoostNormal.load("model.joblib")
 ```
 
-## Saving and Loading
-
-```python
-import openboost as ob
-
-# Train
-model = ob.GradientBoosting(n_trees=100)
-model.fit(X_train, y_train)
-
-# Save
-model.save('my_model.joblib')
-
-# Load
-loaded_model = ob.GradientBoosting.load('my_model.joblib')
-predictions = loaded_model.predict(X_test)
-```
+GPU trees are automatic when CUDA is installed
+(`pip install --pre "openboost[cuda]"`). Force a backend with
+`ob.set_backend("cuda")` or `OPENBOOST_BACKEND=cuda`. See
+[GPU setup](gpu-setup.md).
 
-## Next Steps
+## Next
 
-- [GPU Setup](gpu-setup.md) - Configure GPU acceleration
-- [Uncertainty Quantification](../tutorials/uncertainty.md) - Deep dive into NaturalBoost
-- [Custom Loss Functions](../tutorials/custom-loss.md) - Define your own objectives
+- [How it works](../user-guide/how-it-works.md): distributional regression and varying-coefficient models
+- [Uncertainty tutorial](../tutorials/uncertainty.md): intervals, sampling, NLL
+- [Benchmarks](../benchmarks.md): speed, quality, capability numbers
+- [Custom distributions](../user-guide/naturalboost/custom-distributions.md)
diff --git a/docs/index.md b/docs/index.md
index 53762a3..e028cde 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,110 +1,141 @@
 # OpenBoost
 
-<p align="center" style="font-size: 1.2em;">
-  <strong>The hackable gradient boosting platform — probabilistic predictions, interpretable GAMs, and custom algorithms in readable Python, with CPU and CUDA tree backends.</strong>
+<p align="center" style="font-size: 1.4em;">
+  <strong>GPU gradient boosting for distributional regression.</strong>
+</p>
+
+<p align="center">
+  Every parameter of <code>F(y | x)</code> gets its own tree ensemble, updated
+  with the full <code>K&times;K</code> natural gradient rather than a diagonal
+  approximation. <code>FormulaBoost</code> extends the same engine to
+  varying-coefficient formulas <code>y = f(θ(z), x)</code>.
 </p>
 
 <p align="center">
   <a href="getting-started/quickstart/">Quickstart</a> •
-  <a href="user-guide/models/gradient-boosting/">Models</a> •
-  <a href="user-guide/naturalboost/overview/">NaturalBoost</a> •
-  <a href="api/openboost/">API Reference</a>
+  <a href="user-guide/how-it-works/">How it works</a> •
+  <a href="user-guide/formulaboost/">FormulaBoost</a> •
+  <a href="user-guide/survival/">Weibull AFT</a> •
+  <a href="benchmarks/">Benchmarks</a> •
+  <a href="api/openboost/">API</a>
 </p>
 
 ---
 
-## Why OpenBoost?
+## Why this exists
 
-For standard GBDT, use XGBoost/LightGBM—they're highly optimized C++.
+[Distributional regression](https://doi.org/10.1111/j.1467-9876.2005.00510.x)
+(GAMLSS) fits every parameter of `F(y | x)` as a function of covariates,
+not just the mean. NGBoost does that with boosting, on CPU, for a few
+built-in families.
 
-For GBDT **variants** (probabilistic predictions, interpretable GAMs, custom algorithms), OpenBoost provides reusable Python primitives and a CUDA tree-building path:
+A [varying-coefficient model](https://doi.org/10.1111/j.2517-6161.1993.tb01939.x)
+does the same for a structural formula `y ≈ f(θ(z), x)`. XGBoost custom
+objectives accept only a **diagonal** Hessian, so they cannot represent
+the off-diagonal coupling in `θ`. XGBoost `survival:aft` holds the
+Weibull shape as one global hyperparameter.
 
-- **NaturalBoost**: full-distribution prediction; comparable to NGBoost on the committed CPU benchmark (0.8-1.3x wall-clock, quality within ~1%)
-- **OpenBoostGAM**: interpretable main effects with an optional GPU path; benchmark it on your own workload with the included harness
+OpenBoost fits both classes, censoring included, with histogram trees, a
+full per-sample Fisher or GGN metric, and a GPU path.
 
-Plus: ~20K lines of readable Python. Modify, extend, and build on—no C++ required.
+| | NGBoost | XGBoost | OpenBoost |
+|---|---|---|---|
+| What varies with covariates | Distribution parameters (fixed catalogue) | The mean, or a diagonal custom objective | Distribution or formula parameters |
+| Metric | Natural gradient | Diagonal Hessian only | Fisher / full GGN |
+| GPU trees | No | Yes (mean / AFT location) | Yes |
+| Covariate-dependent Weibull shape `k(z)` | n/a | No (global hyperparameter) | Yes |
+| Formula off-diagonals | n/a | Inexpressible | Full GGN |
 
-## Quick Example
+For ordinary mean regression, use XGBoost or LightGBM. They are faster
+C++.
 
-```python
-import openboost as ob
+## The models
 
-# Standard gradient boosting
-model = ob.GradientBoosting(n_trees=100, max_depth=6)
-model.fit(X_train, y_train)
-predictions = model.predict(X_test)
+=== "NaturalBoost"
 
-# Probabilistic predictions with uncertainty
-prob_model = ob.NaturalBoostNormal(n_trees=100)
-prob_model.fit(X_train, y_train)
-mean = prob_model.predict(X_test)
-lower, upper = prob_model.predict_interval(X_test, alpha=0.1)  # 90% interval
-```
+    Distributional regression. Predict `F(y | x)`, not a point.
 
-## Features
+    ```python
+    import openboost as ob
 
-### :rocket: GPU Accelerated
+    model = ob.NaturalBoostNormal(n_trees=500, max_depth=3, learning_rate=0.03)
+    model.fit(X_train, y_train)
+    mean = model.predict(X_test)
+    lo, hi = model.predict_interval(X_test, alpha=0.1)   # 90% interval
+    ```
 
-Numba CUDA kernels accelerate histogram building and tree construction. Some
-features and model stages remain CPU-only or deliberately fall back to CPU;
-the model guides document those boundaries.
+=== "FormulaBoost"
 
-### :brain: Probabilistic Predictions
+    Varying-coefficient model. Write `y = f(θ, x)`; trees learn `θ(z)`.
 
-NaturalBoost provides full probability distributions with uncertainty quantification. 8 built-in distributions including Normal, Gamma, Tweedie, and Negative Binomial.
+    ```python
+    import numpy as np
+    import openboost as ob
 
-### :snake: All Python
+    def sales(theta, x):
+        a, b = theta
+        return a * x ** (1.0 / (1.0 + np.exp(-b * x)))
 
-~20K lines of readable, hackable code. No C++ compilation needed. Understand and modify the algorithms.
+    model = ob.FormulaBoost(
+        formula=sales, n_params=2, links=("log", "identity"),
+        param_names=("a", "b"), precond="full",
+    )
+    model.fit(Z_train, y_train, model_input=x_train)
+    params = model.predict_params(Z_test)          # per-row a(z), b(z)
+    yhat = model.predict(Z_test, model_input=x_new)
+    ```
 
-### :gear: sklearn Compatible
+=== "WeibullAFT"
 
-Drop-in replacement for scikit-learn pipelines. Works with GridSearchCV, cross_val_score, and Pipeline.
+    Censored Weibull AFT; both scale and shape vary with `z`.
 
-## Installation
+    ```python
+    import openboost as ob
 
-```bash
-pip install --pre openboost
+    model = ob.WeibullAFT(n_trees=300, max_depth=3)
+    model.fit(Z_train, time_train, event=observed)  # 1 = event, 0 = censored
+    params = model.predict_params(Z_test)           # {scale, shape}
+    t_hat = model.predict(Z_test)                   # median time
+    s = model.predict_survival(Z_test, t=5.0)       # S(5 | z)
+    ```
 
-# With GPU support
-pip install --pre "openboost[cuda]"
-```
-
-Without `--pre`, pip installs the older stable release rather than the current
-1.0 release candidate.
+Same trainer, same tree engine. Mean regression is the single-parameter case
+of it: `GradientBoosting`, `OpenBoostGAM`, DART, linear-leaf models, and
+sklearn wrappers ship in the same package.
 
-## What's Included
+## Benchmarks
 
-| Category | Models |
-|----------|--------|
-| **Standard GBDT** | GradientBoosting, MultiClassGradientBoosting, DART |
-| **Interpretable** | OpenBoostGAM, LinearLeafGBDT |
-| **Probabilistic** | NaturalBoostNormal, LogNormal, Gamma, Poisson, StudentT, Tweedie, NegBin |
+Early and incomplete, so read them as directional rather than settled.
 
-## Performance
+| Gate | Where it stands |
+|---|---|
+| Speed vs NGBoost | Seconds on an A100 at sizes where NGBoost, which is CPU-only, takes most of an hour. On CPU the two are near parity, so this is a claim about the GPU tree path. |
+| Quality vs NGBoost | Tied or better NLL on the 8 UCI datasets measured, over 20 paired splits. Three datasets are still unmeasured. |
+| Capability | FormulaBoost and WeibullAFT recover parameter surfaces that a diagonal-Hessian objective cannot express. |
 
-The repository currently includes one auditable third-party comparison:
-`benchmarks/results/ngboost_comparison_20260720.json`.
+XGBoostLSS and LightGBMLSS are the nearest alternatives and are not yet in the
+comparison. Numbers, caveats, and reproduce commands:
+[Benchmarks](benchmarks.md).
 
-| Benchmark | Result |
-|-----------|--------|
-| NaturalBoost vs NGBoost (CPU) | ~parity: 0.8-1.3x, NLL/CRPS/RMSE within ~1% (`ngboost_comparison_20260720.json`) |
+## Who this is for
 
-GPU benchmark harnesses are included, but exact third-party speedup claims are
-not published until the corresponding raw result artifact and environment
-metadata are committed. For standard GBDT, use XGBoost/LightGBM; OpenBoost's
-value is in research-friendly variants and extensibility.
+- Insurance pricing, energy/demand, credit risk: you need `F(y | x)`, not a point
+- Curve / dose / saturation models where the formula is the product and `θ(z)` is what you ship
+- Survival analysis where the Weibull shape should move with covariates
+- Anyone who has been hand-rolling an XGBoost custom objective and hitting the diagonal-Hessian wall
 
-## Who Is OpenBoost For?
+Not the right tool if you want the fastest possible MSE/logloss GBDT. Use
+XGBoost or LightGBM for that.
 
-- **Kaggle Competitors** - Probabilistic predictions that XGBoost can't do
-- **ML Researchers** - Prototype new algorithms in Python
-- **Product teams** - Prototype interpretable or probabilistic models before production hardening
-- **Students** - Actually understand how gradient boosting works
+## Install
 
-## Roadmap
+```bash
+pip install --pre openboost
+pip install --pre "openboost[cuda]"
+```
 
-**Train-many optimization**: Industry workloads often train many models (hyperparameter tuning, CV, per-segment models). XGBoost optimizes for one model fast. OpenBoost plans to enable native optimization for training many models efficiently.
+Without `--pre`, pip installs the older stable release rather than the
+current 1.0 release candidate. See [Installation](getting-started/installation.md).
 
 ## License
 
diff --git a/docs/migration/from-xgboost.md b/docs/migration/from-xgboost.md
index c1a313b..02189d1 100644
--- a/docs/migration/from-xgboost.md
+++ b/docs/migration/from-xgboost.md
@@ -1,6 +1,12 @@
 # Migrating from XGBoost to OpenBoost
 
-This guide helps you transition from XGBoost to OpenBoost with minimal changes.
+Stay on XGBoost or LightGBM for ordinary mean regression. They are
+faster C++. Switch to OpenBoost for **distributional regression**
+(`F(y | x)`), a **varying-coefficient** formula `y = f(θ(z), x)`, or a
+Weibull AFT whose shape varies with covariates.
+
+This page maps XGBoost APIs onto OpenBoost for the overlap, then shows the
+three things XGBoost cannot express.
 
 ## Parameter Mapping
 
@@ -219,79 +225,92 @@ joblib.dump(model, 'model.joblib')
 loaded = joblib.load('model.joblib')
 ```
 
-## Feature Comparison
+## Feature comparison
 
-| Feature | XGBoost | OpenBoost |
-|---------|---------|-----------|
-| GPU Support | ✅ | ✅ |
-| Custom Loss | ⚠️ (requires Python wrapper) | ✅ (native Python) |
-| Uncertainty | ❌ | ✅ (NaturalBoost) |
-| Interpretable GAM | ❌ | ✅ (OpenBoostGAM) |
-| Linear Leaves | ❌ | ✅ (LinearLeafGBDT) |
-| DART | ✅ | ✅ |
-| Growth Strategies | Level-wise | Level-wise, Leaf-wise, Symmetric |
-| GOSS Sampling | ❌ (use LightGBM) | ✅ |
-| Pure Python | ❌ (C++) | ✅ |
+| | XGBoost | OpenBoost |
+|---|---|---|
+| Point-estimate GBDT | Fast C++, the default choice | Works; not the reason to switch |
+| GPU trees | Yes | Yes |
+| Custom loss | Python `obj` callback; **diagonal Hessian only** | Native `(grad, hess)`; FormulaBoost does **full GGN** |
+| Distributional / NGBoost-style | No | `NaturalBoost*` |
+| Formula `y = f(θ, x)` with coupled params | Diagonal custom obj only | `FormulaBoost(precond="full")` |
+| Survival AFT | Location only; scale is a **global** hyperparameter | `WeibullAFT` boosts `λ(z)` **and** `k(z)` |
+| Interpretable GAM | No (use SHAP) | `OpenBoostGAM` |
+| All Python | No | Yes (~20K lines) |
 
-## What OpenBoost Does Better
+## Where OpenBoost adds something
 
-### 1. Uncertainty Quantification
+### 1. A full distribution
 
 ```python
-# XGBoost: Just point predictions
-pred = xgb_model.predict(X_test)  # Single number
+# XGBoost: a point
+pred = xgb_model.predict(X_test)
 
-# OpenBoost: Full distributions
+# OpenBoost: parameters of a distribution
 model = ob.NaturalBoostNormal(n_trees=100)
 model.fit(X_train, y_train)
 mean = model.predict(X_test)
-lower, upper = model.predict_interval(X_test)  # 90% interval
-samples = model.sample(X_test, n_samples=1000)  # Monte Carlo
+lo, hi = model.predict_interval(X_test, alpha=0.1)
+samples = model.sample(X_test, n_samples=1000)
 ```
 
-### 2. Custom Loss Functions
+On the UCI datasets measured so far, NLL is tied or better vs NGBoost, and on
+an A100 NaturalBoost fits in seconds at sizes where CPU-only NGBoost needs
+most of an hour. See [Benchmarks](../benchmarks.md) for the caveats.
 
-```python
-# XGBoost: Requires Python callback wrapper, tricky to get right
-# OpenBoost: Native Python, just return (grad, hess)
-def my_loss(pred, y):
-    grad = pred - y
-    hess = np.ones_like(pred)
-    return grad.astype(np.float32), hess.astype(np.float32)
+### 2. A formula with off-diagonal GGN
 
-model = ob.GradientBoosting(loss=my_loss)
-```
-
-### 3. Interpretable Models
+XGBoost custom objectives cannot represent the off-diagonal of `JᵀJ`.
+That term is what recovers coupled parameters (`b(z)` on the sales curve:
+corr 0.877 vs 0.599). Black-box XGBoost also cannot extrapolate in the
+structural input `x`. FormulaBoost is ~21x better there.
 
 ```python
-# XGBoost: SHAP values (post-hoc, expensive)
-# OpenBoost: Inherently interpretable GAM
-gam = ob.OpenBoostGAM(n_rounds=500)
-gam.fit(X_train, y_train)
-gam.plot_shape_function(0, feature_name="age")
+def sales(theta, x):
+    a, b = theta
+    return a * x ** (1.0 / (1.0 + np.exp(-b * x)))
+
+model = ob.FormulaBoost(
+    formula=sales, n_params=2, links=("log", "identity"),
+    param_names=("a", "b"), precond="full",
+)
+model.fit(Z, y, model_input=x)
+params = model.predict_params(Z)   # the actual deliverable
 ```
 
-### 4. Code Readability
+### 3. Per-row Weibull shape
 
 ```python
-# XGBoost: 200K+ lines of C++
-# OpenBoost: ~6K lines of Python you can actually read and modify
+# XGBoost AFT: one global scale hyperparameter, same k for every row
+# OpenBoost: both λ(z) and k(z) are ensembles
+model = ob.WeibullAFT(n_trees=300, max_depth=3)
+model.fit(Z, time, event=observed)
+params = model.predict_params(Z)           # {scale, shape}
+s = model.predict_survival(Z, t=12.0)
 ```
 
-## What XGBoost Does Better
+On a varying-shape DGP, shape correlation is 0.997; XGBoost has no
+per-row `k` to correlate. Censored NLL is better (0.761 vs 0.830);
+C-index is close (0.680 vs 0.672).
 
-### 1. Raw Speed on CPU
+### 4. A native Python custom loss (point-estimate)
 
-XGBoost is highly optimized C++ and will be faster on CPU for very large datasets.
+```python
+def my_loss(pred, y):
+    grad = pred - y
+    hess = np.ones_like(pred)
+    return grad.astype(np.float32), hess.astype(np.float32)
 
-### 2. Distributed Training (Spark, Dask)
+model = ob.GradientBoosting(loss=my_loss)
+```
 
-XGBoost has mature distributed training support.
+## What XGBoost does better
 
-### 3. Community and Ecosystem
+- **Point-estimate speed on CPU.** Optimized C++. Use it for MSE/logloss.
+- **Distributed training.** Spark / Dask / dedicated cluster runtimes.
+- **Ecosystem.** More examples, more Stack Overflow, more production war stories.
 
-XGBoost has more examples, tutorials, and community support.
+If the job is "fit a GBDT, get a number," do not migrate.
 
 ## Migration Checklist
 
@@ -311,22 +330,23 @@ You can use both libraries during migration:
 import xgboost as xgb
 import openboost as ob
 
-# Keep XGBoost for existing models
+# Keep XGBoost for existing point-estimate models
 xgb_model = xgb.XGBRegressor()
 xgb_model.fit(X_train, y_train)
+xgb_pred = xgb_model.predict(X_test)
 
-# Try OpenBoost for new features
-ob_model = ob.NaturalBoostNormal()  # Uncertainty!
+# OpenBoost for F(y | x)
+ob_model = ob.NaturalBoostNormal()
 ob_model.fit(X_train, y_train)
-
-# Compare predictions
-xgb_pred = xgb_model.predict(X_test)
 ob_pred = ob_model.predict(X_test)
-print(f"Correlation: {np.corrcoef(xgb_pred, ob_pred)[0,1]:.4f}")
+print(f"Correlation: {np.corrcoef(xgb_pred, ob_pred)[0, 1]:.4f}")
 ```
 
-## Getting Help
+## Getting help
 
-- [Quickstart Guide](../getting-started/quickstart.md) - Get started with OpenBoost
-- [Uncertainty Tutorial](../tutorials/uncertainty.md) - Learn NaturalBoost
-- [Custom Loss Tutorial](../tutorials/custom-loss.md) - Define your own objectives
+- [Quickstart](../getting-started/quickstart.md)
+- [How it works](../user-guide/how-it-works.md)
+- [FormulaBoost](../user-guide/formulaboost.md)
+- [Weibull AFT](../user-guide/survival.md)
+- [Uncertainty tutorial](../tutorials/uncertainty.md)
+- [Benchmarks](../benchmarks.md)
diff --git a/docs/tutorials/custom-loss.md b/docs/tutorials/custom-loss.md
index 1609927..2122655 100644
--- a/docs/tutorials/custom-loss.md
+++ b/docs/tutorials/custom-loss.md
@@ -1,6 +1,9 @@
 # Custom Loss Functions
 
-OpenBoost lets you define any loss function in Python. No C++, no recompilation.
+Point-estimate custom objectives: you return `(grad, hess)` in Python, no
+C++. That Hessian is **diagonal**. If the model is a formula with coupled
+parameters, use [FormulaBoost](../user-guide/formulaboost.md)
+(`precond="full"`) instead of rolling an XGBoost-style custom objective.
 
 ## How Gradient Boosting Works
 
@@ -347,7 +350,7 @@ OpenBoost includes these losses out of the box:
 | Huber | `'huber'` | Outlier-robust |
 | Quantile | `'quantile'` | Quantile regression |
 | LogLoss | `'logloss'` | Binary classification |
-| Softmax | (automatic) | Multi-class — used internally by `MultiClassGradientBoosting`, not a `loss=` string |
+| Softmax | (automatic) | Multi-class, used internally by `MultiClassGradientBoosting`, not a `loss=` string |
 | Poisson | `'poisson'` | Count data |
 | Gamma | `'gamma'` | Positive continuous |
 | Tweedie | `'tweedie'` | Zero-inflated positive |
@@ -361,5 +364,6 @@ model = ob.GradientBoosting(n_trees=100, loss='tweedie', tweedie_rho=1.5)
 
 ## Next Steps
 
-- [Uncertainty Quantification](uncertainty.md) - Probabilistic predictions
-- [Migration from XGBoost](../migration/from-xgboost.md) - Switching from XGBoost
+- [FormulaBoost](../user-guide/formulaboost.md): formula parameters, full GGN
+- [Uncertainty Quantification](uncertainty.md)
+- [Migration from XGBoost](../migration/from-xgboost.md)
diff --git a/docs/tutorials/uncertainty.md b/docs/tutorials/uncertainty.md
index 0073df7..b598169 100644
--- a/docs/tutorials/uncertainty.md
+++ b/docs/tutorials/uncertainty.md
@@ -295,5 +295,8 @@ print(f"P(claim > $10k): {prob_large_claim[0]:.1%}")
 
 ## Next Steps
 
-- [Custom Loss Functions](custom-loss.md) - Define your own objectives
-- [Migration from XGBoost](../migration/from-xgboost.md) - Switching from XGBoost
+- [How it works](../user-guide/how-it-works.md)
+- [FormulaBoost](../user-guide/formulaboost.md): when the model is a formula, not a distribution
+- [Weibull AFT](../user-guide/survival.md)
+- [Custom Loss Functions](custom-loss.md)
+- [Migration from XGBoost](../migration/from-xgboost.md)
diff --git a/docs/user-guide/evaluation/metrics.md b/docs/user-guide/evaluation/metrics.md
index 47eaf0e..dda4515 100644
--- a/docs/user-guide/evaluation/metrics.md
+++ b/docs/user-guide/evaluation/metrics.md
@@ -50,7 +50,8 @@ weighted_auc = ob.roc_auc_score(y_true, y_proba, sample_weight=weights)
 
 ## Probabilistic Metrics
 
-For NaturalBoost and distributional models:
+For NaturalBoost, FormulaBoost (via residuals / parameter recovery), and
+WeibullAFT (censored NLL):
 
 ```python
 import openboost as ob
diff --git a/docs/user-guide/formulaboost.md b/docs/user-guide/formulaboost.md
new file mode 100644
index 0000000..916e616
--- /dev/null
+++ b/docs/user-guide/formulaboost.md
@@ -0,0 +1,138 @@
+# FormulaBoost
+
+A [varying-coefficient model](https://doi.org/10.1111/j.2517-6161.1993.tb01939.x):
+trees learn `θ(Z)`, and a structural input `x` enters only through a
+formula you write.
+
+Covariates `Z` determine parameter surfaces `θ(Z)` via trees. A structural
+input `x` (spend, dose, time, …) goes into a formula you write:
+
+```
+y ≈ f(θ(Z), x)
+```
+
+The trees never see `x` as a split feature unless you put it in `Z`. That
+is the point: the shape in `x` is constrained by `f`, so the model can
+extrapolate in `x` while still letting `θ` vary with `Z`.
+
+## Minimal example
+
+```python
+import numpy as np
+import openboost as ob
+
+def sales(theta, x):
+    a, b = theta
+    sigmoid = 1.0 / (1.0 + np.exp(-np.clip(b * x, -30.0, 30.0)))
+    return a * x ** sigmoid
+
+model = ob.FormulaBoost(
+    formula=sales,
+    n_params=2,
+    links=("log", "identity"),   # a > 0, b unconstrained
+    param_names=("a", "b"),
+    precond="full",              # GGN; default
+    n_trees=200,
+    max_depth=3,
+    learning_rate=0.1,
+)
+model.fit(Z_train, y_train, model_input=x_train)
+
+params = model.predict_params(Z_test)              # dict: a(z), b(z)
+yhat = model.predict(Z_test, model_input=x_test)
+```
+
+`formula(theta, x)` receives `theta` as a tuple of `K` arrays (already
+passed through the links) and `x` as a 1-d array. It must return a 1-d
+prediction the same length as `x`.
+
+## Links
+
+Each parameter has a link that maps unconstrained tree scores to the
+domain the formula expects.
+
+| Link | Constrained range | Typical use |
+|---|---|---|
+| `identity` | ℝ | unconstrained coefficients |
+| `log` | (0, ∞) | scales, rates |
+| `softplus` | (0, ∞) | scales, smoother than `log` |
+| `sigmoid` | (0, 1) | fractions, saturations |
+
+`links` must have length `n_params`.
+
+## Preconditioner
+
+FormulaBoost scores the formula with a finite-difference Jacobian `J` and
+preconditions with the generalized Gauss-Newton matrix `JᵀJ` (plus
+Levenberg–Marquardt damping `damp`).
+
+| `precond` | What it does | When |
+|---|---|---|
+| `full` (default) | Invert the `K×K` GGN including off-diagonals | Production. Recovers coupled parameters. |
+| `diag` | Invert only `diag(JᵀJ)` | Faster; similar predictive RMSE, weaker parameter recovery. |
+| `plain` | Raw gradient, no GGN | Debugging only. Diverges on the sales curve. |
+
+XGBoost's custom-objective API accepts only a **diagonal** Hessian. That
+is why `full` vs XGBoost-diag is the capability comparison, not a speed
+contest.
+
+`damp` (default `1.0`) is added to the GGN diagonal. Increase it if
+updates explode; decrease it if training crawls.
+
+## Fit signature
+
+```python
+model.fit(
+    Z, y, model_input=x,
+    sample_weight=None,
+    eval_set=[(Z_val, y_val, x_val)],
+    callbacks=[ob.EarlyStopping(patience=20)],
+    early_stopping_rounds=20,
+)
+```
+
+`eval_set` entries are always `(X, y, model_input)`, three arrays.
+
+## What to evaluate
+
+On the sales-curve benchmark (200K rows, train `x ∈ [0.25, 2.5]`, extrap
+`x ∈ [3, 5]`):
+
+| | test RMSE | extrap RMSE | corr `b(z)` |
+|---|--:|--:|--:|
+| FormulaBoost `full` | 0.139 | **0.183** | **0.877** |
+| FormulaBoost `diag` | 0.130 | 0.181 | 0.730 |
+| black-box GBDT | 0.127 | 3.872 | n/a |
+| XGBoost custom (diag Hess) | 0.129 | 0.203 | 0.599 |
+| global `(a, b)` | 1.641 | 4.396 | n/a |
+
+Honest reading:
+
+- **Extrapolation** is the product claim (~21× vs black-box). The formula
+  constrains the shape in `x`; a GBDT that splits on `x` does not.
+- **Parameter recovery** of `b(z)` is what `full` GGN buys. If you only
+  care about in-sample RMSE, `diag` is enough.
+- `plain` is not a baseline you should ship (extrap RMSE ~3.9).
+
+Reproduce: `uv run modal run benchmarks/bench_formula.py`. Full notes:
+[Benchmarks](../benchmarks.md).
+
+## Tips
+
+- Keep `K` small. The GGN is `K×K` per row; two to five parameters is the
+  intended range.
+- Put only **covariates** in `Z`. Put the structural axis in
+  `model_input`. Mixing them usually destroys extrapolation.
+- Start with a global fit in your head (what would one `(a, b)` be?) so
+  you can tell whether the surfaces are recovering anything.
+- Shallower trees (`max_depth=3`) and more rounds, same as NaturalBoost.
+
+## Persistence
+
+```python
+model.save("formula.joblib")
+loaded = ob.FormulaBoost.load("formula.joblib")
+```
+
+The formula callable is pickled with the model. Keep it importable under
+the same module path if you move files between save and load.
diff --git a/docs/user-guide/how-it-works.md b/docs/user-guide/how-it-works.md
new file mode 100644
index 0000000..2d28461
--- /dev/null
+++ b/docs/user-guide/how-it-works.md
@@ -0,0 +1,109 @@
+# Distributional and varying-coefficient models
+
+OpenBoost estimates the parameters of a statistical model as functions of
+covariates, by gradient boosting.
+
+That is **distributional regression** when the model is a conditional
+distribution `F(y | x)`: each parameter (location, scale, shape, ...)
+gets its own additive predictor. The statistics literature calls the
+GAM version [GAMLSS](https://doi.org/10.1111/j.1467-9876.2005.00510.x);
+the boosting version is [gamboostLSS](https://doi.org/10.1214/12-AOAS580)
+/ NGBoost. NaturalBoost and WeibullAFT live here.
+
+It is a **varying-coefficient model**
+([Hastie & Tibshirani, 1993](https://doi.org/10.1111/j.2517-6161.1993.tb01939.x))
+when you supply a structural formula `y ≈ f(θ(z), x)`. The trees learn
+`θ(z)`; `x` enters only through `f`. That is FormulaBoost.
+
+```
+z  --trees-->  F (n × K raw scores)  --link-->  θ  --f-->  prediction  --loss-->  scalar
+```
+
+`z` are the covariates the trees split on. `θ = (θ₁, …, θₖ)` is the
+parameter vector of the statistical model. Each `θᵢ` is its own ensemble.
+
+The per-round update is a (damped) generalized Gauss-Newton / natural
+gradient step:
+
+```
+dᵢ = (Jᵢᵀ Mᵢ Jᵢ + λI)⁻¹ gᵢ     # K × K per sample, K small
+```
+
+Models are configurations of `f`, the loss, and the metric `M`.
+
+| Model | Statistical object | Loss | Metric `M` | `K` |
+|---|---|---|---|---|
+| [NaturalBoost](naturalboost/overview.md) | Conditional distribution | NLL | analytic Fisher | # distribution params |
+| [WeibullAFT](survival.md) | Censored Weibull AFT | censored NLL | expected Fisher | 2 (scale, shape) |
+| [FormulaBoost](formulaboost.md) | Varying-coefficient formula | MSE | GGN (`plain` / `diag` / `full`) | # formula params |
+| `GradientBoosting` | Conditional mean | MSE / logloss / … | scalar Hessian | 1 |
+
+`K = 1` is ordinary mean regression. `K > 1` is the point of this library.
+
+## What this is not
+
+Not a faster drop-in for XGBoost/LightGBM mean regression. Those are
+optimized C++ and should stay the default for MSE/logloss.
+
+Not only "faster NGBoost." NGBoost is distributional regression for a
+fixed set of families, with exact-split trees, on CPU. OpenBoost does
+the same job with histogram trees and a GPU path, and also fits
+varying-coefficient formulas and a censored Weibull whose shape depends
+on covariates, neither of which NGBoost can write down.
+
+On the overlap (Normal NLL on UCI) quality is tied-or-better and GPU
+training is much faster: [Benchmarks](../benchmarks.md).
+
+## Why the metric matters
+
+Raw gradients of a multi-parameter likelihood are often badly scaled.
+On the sales-curve formula, `precond="plain"` diverges (extrapolation RMSE
+~3.9). Diagonal GGN is usable. Full GGN recovers the hard coefficient
+(`b(z)` corr 0.877 vs 0.730 diag vs 0.599 for an XGBoost custom objective,
+which cannot represent off-diagonals).
+
+Same story for Weibull AFT: the observed Hessian's scale term explodes
+when `λ` is wrong. The **expected** Fisher does not, and is what
+`WeibullAFT` uses.
+
+`precond="full"` is the default for FormulaBoost. Do not turn it off
+unless you are debugging.
+
+## Shared training API
+
+NaturalBoost, FormulaBoost, and WeibullAFT share one trainer:
+
+```python
+model.fit(
+    X, y,                          # FormulaBoost also needs model_input=x
+    eval_set=[...],                # early stopping / logging
+    callbacks=[ob.EarlyStopping(patience=20)],
+    early_stopping_rounds=20,
+)
+params = model.predict_params(X)   # dict of per-row parameter arrays
+```
+
+- NaturalBoost / WeibullAFT: `eval_set` is `(X_val, y_val)` (plus `event`
+  for AFT).
+- FormulaBoost: `eval_set` is `(X_val, y_val, model_input_val)`.
+
+GPU: histogram trees run on CUDA when `openboost[cuda]` is installed.
+NaturalBoost Normal/Poisson have device gradient kernels. FormulaBoost GGN
+is currently host-side (still faster than XGBoost-diag at 200K on A100).
+
+## Choosing a model
+
+| You have | Use |
+|---|---|
+| `y` should be a distribution given `X` | `NaturalBoost*` |
+| A known curve `y = f(θ, x)` with `θ` depending on other features | `FormulaBoost` |
+| Right-censored times, Weibull, shape should depend on covariates | `WeibullAFT` |
+| A custom NLL that is still a distribution | `NaturalBoost` + [custom distribution](naturalboost/custom-distributions.md) |
+| Ordinary mean regression | `GradientBoosting` (or XGBoost / LightGBM) |
+
+## Next
+
+- [NaturalBoost](naturalboost/overview.md)
+- [FormulaBoost](formulaboost.md)
+- [Weibull AFT](survival.md)
+- [Benchmarks](../benchmarks.md)
diff --git a/docs/user-guide/model-persistence.md b/docs/user-guide/model-persistence.md
index 30579ec..6138390 100644
--- a/docs/user-guide/model-persistence.md
+++ b/docs/user-guide/model-persistence.md
@@ -23,17 +23,16 @@ predictions = loaded_model.predict(X_test)
 
 All models support save/load:
 
-- `GradientBoosting`
-- `MultiClassGradientBoosting`
-- `DART`
-- `OpenBoostGAM`
-- `NaturalBoostNormal`, `NaturalBoostGamma`, etc.
-- `LinearLeafGBDT`
+- `NaturalBoostNormal`, `NaturalBoostGamma`, `NaturalBoostTweedie`, …
+- `FormulaBoost` (pickles the formula callable; keep it importable)
+- `WeibullAFT`
+- `GradientBoosting`, `MultiClassGradientBoosting`, `DART`
+- `OpenBoostGAM`, `LinearLeafGBDT`
 
 ## Using joblib/pickle Directly
 
 > **Security warning:** joblib and pickle deserialization executes arbitrary
-> code. Only load model files from sources you trust — never from untrusted
+> code. Only load model files from sources you trust, never from untrusted
 > uploads or downloads. The same caveat applies to `ob.load()` /
 > `Model.load()`, which emit a runtime `UserWarning` for this reason.
 
diff --git a/docs/user-guide/models/dart.md b/docs/user-guide/models/dart.md
index 5c61149..3d36781 100644
--- a/docs/user-guide/models/dart.md
+++ b/docs/user-guide/models/dart.md
@@ -1,6 +1,9 @@
 # DART
 
-Dropout Additive Regression Trees - a regularization technique that randomly drops trees during training.
+Dropout Additive Regression Trees, a regularization technique that randomly
+drops trees during training. Mean-regression GBDT (`K = 1`). For
+distributional regression or a varying-coefficient formula, see
+[How it works](../how-it-works.md).
 
 ## Why DART?
 
diff --git a/docs/user-guide/models/gam.md b/docs/user-guide/models/gam.md
index 615e374..55f1eaa 100644
--- a/docs/user-guide/models/gam.md
+++ b/docs/user-guide/models/gam.md
@@ -1,6 +1,8 @@
 # OpenBoostGAM
 
-GPU-accelerated Generalized Additive Model - interpretable machine learning with feature-level explanations.
+GPU-accelerated Generalized Additive Model: interpretable main effects
+(and optional pairwise interactions). Point-estimate. For a distribution
+or a formula, see [How it works](../how-it-works.md).
 
 ## Why GAM?
 
@@ -12,7 +14,7 @@ prediction = f₁(x₁) + f₂(x₂) + ... + fₙ(xₙ) + intercept
 
 This means you can visualize exactly how each feature affects the prediction.
 
-**Scope:** By default OpenBoostGAM learns *main effects only* — one shape
+**Scope:** By default OpenBoostGAM learns *main effects only*: one shape
 function per feature. Setting `interactions=k` adds `k` pairwise interaction
 terms (GA2M-style, like InterpretML's EBM), which closes part of the accuracy
 gap to EBM while keeping every term inspectable. Be aware that the interaction
@@ -75,7 +77,7 @@ behavior. `n_trees` is accepted as an alias for `n_rounds`. There is no
 `max_depth` parameter: each round applies a regularized Newton update to every
 bin of every feature's shape function (a per-feature lookup table), so
 smoothness is controlled by `learning_rate`, `n_rounds`, `reg_lambda`, and
-`smoothing` — not tree depth.
+`smoothing`, not tree depth.
 
 ## Pairwise Interactions (GA2M)
 
@@ -96,7 +98,7 @@ the `(bin_i, bin_j)` grid and raises `KeyError` for pairs that were not
 selected. Interaction terms are used automatically by `predict`.
 
 **Honest scope note:** interactions close part of the gap to EBM (which enables
-them by default) but the interaction stage always trains on CPU — expect the
+them by default) but the interaction stage always trains on CPU, so expect the
 GPU speedup story to apply to the main-effects stage only.
 
 ## Smoothing and Monotone Constraints
@@ -113,7 +115,7 @@ gam.fit(X_train, y_train)
 
 - `smoothing` solves a fused-ridge (first-difference penalty) system per round,
   so occupied bins anchor the shape while empty bins interpolate between their
-  neighbors. It applies to ordinal (numeric) bins only — categorical features,
+  neighbors. It applies to ordinal (numeric) bins only. Categorical features,
   the missing-value bin, and 2D interaction tables are not smoothed.
 - `monotone` projects the accumulated shape function onto the constraint after
   every round with count-weighted isotonic regression (PAVA). Useful when
@@ -132,7 +134,7 @@ gam.fit(
     early_stopping_rounds=50,   # sugar for EarlyStopping(patience=50, restore_best=True)
 )
 
-gam.evals_result_     # {'eval_0': {'mse': [...]}} — per-round history per eval set
+gam.evals_result_     # {'eval_0': {'mse': [...]}}: per-round history per eval set
 gam.best_iteration_   # set when early stopping is used
 gam.best_score_
 ```
@@ -140,7 +142,7 @@ gam.best_score_
 Early stopping monitors the **last** eval set; with `restore_best=True` (the
 default via `early_stopping_rounds`) the shape tables are restored to the best
 round. With `interactions > 0`, the main-effect and interaction rounds form
-one monitored sequence — stopping during the main-effect phase skips the
+one monitored sequence, so stopping during the main-effect phase skips the
 interaction phase.
 
 **Backend note:** smoothing, monotone projection, the interaction stage, and
@@ -192,5 +194,5 @@ plt.savefig("gam_explanations.png")
 1. **Use more rounds** with lower learning rate for smoother shape functions
 2. **Normalize features** for easier interpretation
 3. **Check shape functions** for unexpected patterns (data issues)
-4. **Compare against a standard GBDT** — if it beats the GAM by a wide
+4. **Compare against a standard GBDT**: if it beats the GAM by a wide
    margin, your data likely has interactions the GAM cannot capture
diff --git a/docs/user-guide/models/gradient-boosting.md b/docs/user-guide/models/gradient-boosting.md
index aa9b280..b8d5386 100644
--- a/docs/user-guide/models/gradient-boosting.md
+++ b/docs/user-guide/models/gradient-boosting.md
@@ -1,6 +1,9 @@
 # Gradient Boosting
 
-The core gradient boosting model for regression and binary classification.
+Mean-regression GBDT for regression and binary classification (`K = 1`):
+one parameter, scalar Hessian. For distributional regression or a
+varying-coefficient formula, start at
+[How it works](../how-it-works.md) instead.
 
 The examples on this page share this setup:
 
diff --git a/docs/user-guide/models/linear-leaf.md b/docs/user-guide/models/linear-leaf.md
index 3432d45..acf5583 100644
--- a/docs/user-guide/models/linear-leaf.md
+++ b/docs/user-guide/models/linear-leaf.md
@@ -1,6 +1,9 @@
 # Linear Leaf GBDT
 
-Trees with linear models in the leaves instead of constant values. Better for extrapolation and smooth relationships.
+Trees with linear models in the leaves instead of constant values. Better
+for *local* linear extrapolation. If the shape in `x` is a known formula
+and you want parameter surfaces `θ(z)`, use
+[FormulaBoost](../formulaboost.md) instead.
 
 ## Why Linear Leaves?
 
@@ -91,7 +94,7 @@ model.fit(
     early_stopping_rounds=20,   # sugar for EarlyStopping(patience=20, restore_best=True)
 )
 
-model.evals_result_     # {'eval_0': {'mse': [...]}} — per-round history per eval set
+model.evals_result_     # {'eval_0': {'mse': [...]}}: per-round history per eval set
 model.best_iteration_   # set when early stopping is used
 model.best_score_
 ```
diff --git a/docs/user-guide/models/multiclass.md b/docs/user-guide/models/multiclass.md
index f088838..a0acb21 100644
--- a/docs/user-guide/models/multiclass.md
+++ b/docs/user-guide/models/multiclass.md
@@ -1,6 +1,8 @@
 # Multi-class Classification
 
-For classification problems with more than 2 classes.
+Softmax GBDT for 3+ classes. Point-estimate (`K = n_classes` logits, not
+a parametric family). For a distribution over a real-valued target, use
+[NaturalBoost](../naturalboost/overview.md).
 
 ## Basic Usage
 
diff --git a/docs/user-guide/naturalboost/custom-distributions.md b/docs/user-guide/naturalboost/custom-distributions.md
index 830b5a9..d055381 100644
--- a/docs/user-guide/naturalboost/custom-distributions.md
+++ b/docs/user-guide/naturalboost/custom-distributions.md
@@ -90,7 +90,7 @@ lower, upper = model.predict_interval(X_test, alpha=0.1)
 
 ## Gradient Computation: JAX or Numerical
 
-You do not provide gradients — they are computed automatically from `nll_fn`:
+You do not provide gradients; they are computed automatically from `nll_fn`:
 
 - **JAX** (used automatically when `jax` is installed): exact autodiff of the
   NLL composed with the link functions, vectorized with `jax.vmap`.
@@ -142,3 +142,8 @@ print(ob.list_distributions())
 normal = ob.get_distribution('normal')
 gamma = ob.get_distribution('gamma')
 ```
+
+If the object you want to boost is a formula `y = f(θ, x)` rather than a
+probability distribution, skip this page and use
+[FormulaBoost](../formulaboost.md). For right-censored Weibull,
+use [WeibullAFT](../survival.md).
diff --git a/docs/user-guide/naturalboost/distributions.md b/docs/user-guide/naturalboost/distributions.md
index 96801de..77c3d25 100644
--- a/docs/user-guide/naturalboost/distributions.md
+++ b/docs/user-guide/naturalboost/distributions.md
@@ -1,6 +1,9 @@
 # Distributions
 
-NaturalBoost supports multiple probability distributions for different data types.
+NaturalBoost families. Pick one that matches the support of `y`. For a
+formula that is not a distribution, use
+[FormulaBoost](../formulaboost.md). For censored survival, use
+[WeibullAFT](../survival.md).
 
 ## Choosing a Distribution
 
diff --git a/docs/user-guide/naturalboost/overview.md b/docs/user-guide/naturalboost/overview.md
index b5e34f2..7e6a08a 100644
--- a/docs/user-guide/naturalboost/overview.md
+++ b/docs/user-guide/naturalboost/overview.md
@@ -1,6 +1,16 @@
-# NaturalBoost Overview
+# NaturalBoost
 
-NaturalBoost predicts full probability distributions instead of just point estimates, giving you uncertainty bounds on your predictions.
+NaturalBoost is **distributional regression** via boosting: each
+parameter of `F(y | x)` (`loc`, `scale`, …) is its own ensemble, stepped
+with natural gradient. That is the
+[GAMLSS](https://doi.org/10.1111/j.1467-9876.2005.00510.x) / NGBoost
+model class, with a GPU histogram-tree path.
+
+For a structural formula that is not a distribution, use
+[FormulaBoost](../formulaboost.md) (varying-coefficient). For
+right-censored Weibull survival, use
+[WeibullAFT](../survival.md). Shared engine:
+[How it works](../how-it-works.md).
 
 ## Why Uncertainty Matters
 
@@ -81,33 +91,22 @@ prob_exceed = np.mean(samples > threshold, axis=0)  # P(Y > 10)
 q90 = np.percentile(samples, 90, axis=0)
 ```
 
-## Performance vs NGBoost
+## vs NGBoost
 
-Measured head-to-head against NGBoost 0.5.11 (both with a Normal distribution,
-natural gradient, and an identical budget: 500 boosting rounds, learning rate 0.03,
-depth-3 trees, seed 42, same train/test splits). This is a CPU-vs-CPU comparison —
-NGBoost is CPU-only, and OpenBoost's GPU tree path is deliberately **not** measured
-here. NLL and CRPS use the same closed-form Gaussian formulas for both models.
-Full configs, metrics, and library versions are committed in
-`benchmarks/results/ngboost_comparison_20260720.json`. Reproduce with:
+On GPU (Modal A100, 90K rows, heteroscedastic Normal) NaturalBoost fits in
+2.21s against NGBoost's 2716s, with NLL tied (2.108 vs 2.102). NGBoost has no
+GPU implementation. This is one configuration from an early benchmark run, not
+a settled result.
 
-```bash
-OPENBOOST_BACKEND=cpu uv run --with ngboost python benchmarks/bench_ngboost_comparison.py
-```
+On the NGBoost-paper UCI suite (20 paired splits, same budget) OpenBoost is
+tied-or-better on every dataset that completed; significant NLL wins on
+kin8nm, protein, and california; no significant loss.
 
-| Dataset | Fit time OB / NGB | Test NLL OB / NGB | CRPS OB / NGB | RMSE OB / NGB |
-|---------|-------------------|-------------------|---------------|---------------|
-| Synthetic heteroscedastic, 10K | **16.4s** / 18.8s (1.15x) | **2.124** / 2.134 | **1.169** / 1.174 | **2.145** / 2.149 |
-| Synthetic heteroscedastic, 50K | **74.1s** / 95.3s (1.29x) | 2.122 / **2.116** | 1.184 / **1.177** | 2.165 / **2.152** |
-| California Housing, 20.6K | 30.6s / **25.0s** (0.82x) | **0.572** / 0.575 | **0.255** / 0.256 | **0.518** / 0.521 |
+On CPU the two libraries are ~parity (0.8–1.3×). The speed claim is the
+GPU tree path, not a faster CPU NGBoost clone.
 
-The honest read: on CPU the two libraries are comparable. NaturalBoost's
-histogram-based trees are modestly faster on the larger synthetic dataset (1.29x at
-50K samples), while NGBoost was faster on California Housing (0.82x). Prediction
-quality is essentially tied — NGBoost slightly wins NLL/CRPS/RMSE on the 50K
-synthetic dataset; NaturalBoost slightly wins on the other two. NaturalBoost's main
-differentiators are the GPU tree path and the wider distribution/custom-distribution
-support, not raw CPU speed.
+Full tables, caveats, and reproduce commands:
+[Benchmarks](../../benchmarks.md).
 
 ## Best Practices
 
@@ -122,3 +121,11 @@ model = ob.NaturalBoostNormal(
     learning_rate=0.05,  # Lower LR
 )
 ```
+
+## Related
+
+- [How it works](../how-it-works.md): the shared engine
+- [FormulaBoost](../formulaboost.md): when `f` is a formula, not a distribution
+- [Weibull AFT](../survival.md): censored survival
+- [Custom distributions](custom-distributions.md)
+- [Benchmarks](../../benchmarks.md)
diff --git a/docs/user-guide/sklearn-integration.md b/docs/user-guide/sklearn-integration.md
index 6d9b3a7..1dbfcd8 100644
--- a/docs/user-guide/sklearn-integration.md
+++ b/docs/user-guide/sklearn-integration.md
@@ -1,6 +1,8 @@
 # sklearn Integration
 
-OpenBoost provides sklearn-compatible wrappers for seamless integration with scikit-learn pipelines.
+OpenBoost provides sklearn-compatible wrappers for mean-regression GBDT
+and NaturalBoost. `FormulaBoost` and `WeibullAFT` are core APIs (no
+sklearn estimator yet); use them directly.
 
 ## Available Wrappers
 
diff --git a/docs/user-guide/survival.md b/docs/user-guide/survival.md
new file mode 100644
index 0000000..bf7f1b5
--- /dev/null
+++ b/docs/user-guide/survival.md
@@ -0,0 +1,107 @@
+# Weibull AFT
+
+Boosted Weibull accelerated-failure-time model with right censoring.
+
+Both the scale `λ(z)` and the shape `k(z)` are boosting ensembles over
+covariates `z`, trained on a right-censored negative log-likelihood.
+That is the capability NGBoost and XGBoost AFT do not have:
+
+- NGBoost has no censored likelihood.
+- XGBoost `survival:aft` learns a location and holds the distribution
+  scale as **one global hyperparameter**, so every row shares the same
+  Weibull shape.
+
+## Minimal example
+
+```python
+import openboost as ob
+
+model = ob.WeibullAFT(n_trees=300, max_depth=3, learning_rate=0.1)
+model.fit(Z_train, time_train, event=observed)   # event: 1 seen, 0 censored
+
+params = model.predict_params(Z_test)            # {scale, shape} per row
+t_hat = model.predict(Z_test)                    # median survival time
+q90 = model.predict_quantile(Z_test, q=0.9)
+s = model.predict_survival(Z_test, t=12.0)       # S(12 | z)
+nll = model.nll(Z_test, time_test, event=observed_test)
+```
+
+`y` is observed time and must be strictly positive. `event` is 1 for an
+observed failure and 0 for right-censored. If `event` is omitted, every
+row is treated as observed.
+
+## Likelihood and metric
+
+Weibull survival function:
+
+```
+S(t | z) = exp( − (t / λ(z))^{k(z)} )
+```
+
+Trees produce unconstrained scores `(u, v)`; links map them to
+`λ = exp(u)`, `k = exp(v)`. The trainer steps in that unconstrained space
+using the **expected Fisher information** of the censored Weibull, not
+the observed Hessian. The observed scale term `k² z` blows up when `λ` is
+wrong and freezes the update. The expected Fisher does not.
+
+`damp` (default `1.0`) is Levenberg–Marquardt damping on that `2×2`
+matrix. Increase it if NLL spikes in the first rounds.
+
+## Fit signature
+
+```python
+model.fit(
+    Z, time, event=observed,
+    eval_set=[(Z_val, time_val, event_val)],
+    callbacks=[ob.EarlyStopping(patience=20)],
+    early_stopping_rounds=20,
+)
+```
+
+`eval_set` entries are `(X, y)` or `(X, y, event)`. The logged metric is
+censored NLL.
+
+## Predictions
+
+| Method | Returns |
+|---|---|
+| `predict_params(X)` | `dict` with `scale` (`λ`) and `shape` (`k`), shape `(n,)` |
+| `predict(X)` | median time, `λ (ln 2)^{1/k}` |
+| `predict_median(X)` | same as `predict` |
+| `predict_quantile(X, q)` | time at which `P(T ≤ t) = q` |
+| `predict_survival(X, t)` | `S(t \| z)` for a scalar or per-row `t` |
+| `nll(X, y, event=)` | mean censored negative log-likelihood |
+
+## vs XGBoost `survival:aft`
+
+Synthetic DGP where **both** `λ(z)` and `k(z)` vary, ~35% right-censoring,
+200K rows, 300 rounds:
+
+| | C-index | NLL | 80% coverage | shape corr | fit |
+|---|--:|--:|--:|--:|--:|
+| OpenBoost `WeibullAFT` | **0.680** | **0.761** | 0.803 | **0.997** | 5.4s |
+| XGBoost `survival:aft` (`extreme`) | 0.672 | 0.830 | 0.852 | n/a (global `k=1.34`) | 11.7s |
+| global constant | 0.500 | 0.896 | 0.810 | n/a | n/a |
+
+C-index is close, since ranking mostly follows the scale. The NLL gap and the
+shape correlation are the capability: OpenBoost recovers `k(z)`, XGBoost
+cannot represent it. Coverage of the 80% interval is nearer the nominal
+0.80 (XGBoost over-covers).
+
+Reproduce: `uv run modal run benchmarks/bench_survival.py`. Notes:
+[Benchmarks](../benchmarks.md).
+
+## Tips
+
+- Times must be `> 0`. Shift or clip before fitting.
+- `event` dtype does not matter as long as it is 0/1.
+- Shallower trees (`max_depth=3`) and more rounds, same as NaturalBoost.
+- If NLL diverges, raise `damp` or lower `learning_rate` before adding
+  trees.
+
+## Persistence
+
+```python
+model.save("aft.joblib")
+loaded = ob.WeibullAFT.load("aft.joblib")
+```
diff --git a/docs/user-guide/training/callbacks.md b/docs/user-guide/training/callbacks.md
index 3013e01..407ee02 100644
--- a/docs/user-guide/training/callbacks.md
+++ b/docs/user-guide/training/callbacks.md
@@ -2,6 +2,11 @@
 
 Control training with callbacks for early stopping, logging, checkpointing, and more.
 
+NaturalBoost, FormulaBoost, WeibullAFT, and the mean-regression models
+all accept `callbacks=` and `eval_set=` on `fit()`.
+`FormulaBoost` eval tuples are `(X, y, model_input)`; `WeibullAFT` tuples
+are `(X, y[, event])`.
+
 ## Available Callbacks
 
 | Callback | Purpose |
diff --git a/examples/README.md b/examples/README.md
index c3f3d55..077eb49 100644
--- a/examples/README.md
+++ b/examples/README.md
@@ -1,213 +1,94 @@
 # OpenBoost Examples
 
-This directory contains runnable examples demonstrating OpenBoost's capabilities.
-
-## Quick Start
+Runnable scripts. Distributional regression first; mean-regression GBDT after.
 
 ```bash
-# Run any example
-uv run python examples/basic_regression.py
-
-# Or with standard Python
-python examples/basic_regression.py
+uv run python examples/uncertainty_quantification.py
 ```
 
-## Examples Overview
-
-| Example | Description | Key Features |
-|---------|-------------|--------------|
-| [basic_regression.py](basic_regression.py) | Standard gradient boosting for regression | `GradientBoosting`, callbacks, feature importance |
-| [binary_classification.py](binary_classification.py) | Binary classification with probability outputs | `OpenBoostClassifier`, ROC AUC, calibration |
-| [multiclass_classification.py](multiclass_classification.py) | Multi-class classification with softmax | `MultiClassGradientBoosting`, confusion matrix |
-| [uncertainty_quantification.py](uncertainty_quantification.py) | Probabilistic predictions with uncertainty | `NaturalBoostNormal`, prediction intervals, CRPS |
-| [kaggle_insurance.py](kaggle_insurance.py) | Insurance claims with Tweedie distribution | `NaturalBoostTweedie`, zero-inflated data |
-| [kaggle_sales.py](kaggle_sales.py) | Sales forecasting with Negative Binomial | `NaturalBoostNegBin`, overdispersed counts |
-| [custom_loss.py](custom_loss.py) | Custom loss functions | Quantile, Huber, asymmetric losses |
-| [gpu_training.py](gpu_training.py) | GPU acceleration guide | Backend selection, benchmarking |
-| [gam_explainability.py](gam_explainability.py) | Interpretable GAM models | `OpenBoostGAM`, shape functions |
-| [sklearn_pipeline.py](sklearn_pipeline.py) | sklearn Pipeline integration | `Pipeline`, `GridSearchCV`, preprocessing |
-| [model_persistence.py](model_persistence.py) | Saving and loading models | `save()`, `load()`, checkpointing |
-
-## Detailed Descriptions
-
-### Basic Regression (`basic_regression.py`)
-
-Learn the fundamentals of OpenBoost with a standard regression task.
-
-**Topics covered:**
-- Training `GradientBoosting` with various hyperparameters
-- Using callbacks (`EarlyStopping`, `Logger`)
-- Computing feature importances
-- sklearn-compatible API with `OpenBoostRegressor`
-- Cross-validation utilities
-
-### Binary Classification (`binary_classification.py`)
-
-Train a binary classifier with probability calibration analysis.
-
-**Topics covered:**
-- Binary classification with `logloss` objective
-- `OpenBoostClassifier` sklearn wrapper
-- ROC AUC, precision, recall, F1 metrics
-- Calibration analysis (Brier score, ECE)
-- Out-of-fold probability predictions
+## Overview
 
-### Uncertainty Quantification (`uncertainty_quantification.py`)
+| Example | Description | Key APIs |
+|---------|-------------|----------|
+| [uncertainty_quantification.py](uncertainty_quantification.py) | Full distribution, intervals, CRPS | `NaturalBoostNormal` |
+| [kaggle_insurance.py](kaggle_insurance.py) | Zero-inflated claims | `NaturalBoostTweedie` |
+| [kaggle_sales.py](kaggle_sales.py) | Overdispersed counts | `NaturalBoostNegBin` |
+| [basic_regression.py](basic_regression.py) | Point-estimate regression | `GradientBoosting` |
+| [binary_classification.py](binary_classification.py) | Binary classifier | `OpenBoostClassifier` |
+| [multiclass_classification.py](multiclass_classification.py) | Softmax multi-class | `MultiClassGradientBoosting` |
+| [custom_loss.py](custom_loss.py) | Quantile / Huber / asymmetric | custom `loss=` |
+| [gpu_training.py](gpu_training.py) | Backend selection | `set_backend` |
+| [gam_explainability.py](gam_explainability.py) | Interpretable main effects | `OpenBoostGAM` |
+| [sklearn_pipeline.py](sklearn_pipeline.py) | `Pipeline` / `GridSearchCV` | `OpenBoostRegressor` |
+| [model_persistence.py](model_persistence.py) | `save` / `load` | `PersistenceMixin` |
 
-The power of NaturalBoost: full probability distributions, not just point estimates!
+FormulaBoost and WeibullAFT do not have example scripts yet; copy from the
+docs:
 
-**Topics covered:**
-- Training `NaturalBoostNormal` for probabilistic predictions
-- Prediction intervals (90%, 80%, 50%)
-- Quantile predictions
-- Sampling from predicted distributions
-- Proper scoring rules (CRPS, NLL)
-- Heteroscedastic uncertainty
+- [FormulaBoost](https://jxucoder.github.io/openboost/user-guide/formulaboost/)
+- [Weibull AFT](https://jxucoder.github.io/openboost/user-guide/survival/)
+- [Quickstart](https://jxucoder.github.io/openboost/getting-started/quickstart/)
 
-### Kaggle Insurance (`kaggle_insurance.py`)
+## Uncertainty (`uncertainty_quantification.py`)
 
-Tweedie distribution for insurance claim prediction (like Porto Seguro, Allstate).
+NaturalBoost: intervals, quantiles, sampling, CRPS / NLL, heteroscedastic
+noise.
 
-**Topics covered:**
-- `NaturalBoostTweedie` for zero-inflated positive continuous data
-- Risk segmentation analysis
-- Probability of large claims
-- Individual risk assessment
-- Comparison with simple MSE model
+## Insurance (`kaggle_insurance.py`)
 
-### Kaggle Sales (`kaggle_sales.py`)
+`NaturalBoostTweedie` for zero-inflated positive claims. Risk
+segmentation, P(large claim), vs a plain MSE GBDT.
 
-Negative Binomial for sales/demand forecasting (like Rossmann, Bike Sharing).
+## Sales (`kaggle_sales.py`)
 
-**Topics covered:**
-- `NaturalBoostNegBin` for overdispersed count data
-- Inventory planning (service levels)
-- Day-of-week and promotional effects
-- Probability of high demand
-- Comparison with Poisson model
+`NaturalBoostNegBin` for overdispersed counts. Service levels, vs Poisson.
 
-### Custom Loss Functions (`custom_loss.py`)
+## Point-estimate GBDT
 
-Build any loss function you need!
+`basic_regression.py`, `binary_classification.py`,
+`multiclass_classification.py` cover `GradientBoosting` / sklearn
+wrappers, callbacks, and feature importance. Use these when you want a
+number, not `F(y | x)`. For production MSE/logloss at scale,
+XGBoost / LightGBM are still faster C++.
 
-**Topics covered:**
-- Quantile regression for different percentiles
-- Huber loss for outlier robustness
-- Asymmetric loss for business costs
-- Log-cosh smooth approximation
-- How to write custom loss functions
+## Custom loss (`custom_loss.py`)
 
-### GPU Training (`gpu_training.py`)
+Quantile, Huber, asymmetric, log-cosh. Each returns `(grad, hess)`. For a
+**formula** with coupled parameters, this is the wrong tool; use
+FormulaBoost (full GGN), not a diagonal custom objective.
 
-Get the most out of GPU acceleration.
+## GPU (`gpu_training.py`)
 
-**Topics covered:**
-- Automatic GPU detection
-- Manual backend selection
-- Performance benchmarking
-- Best practices for GPU training
-- Multi-GPU training overview
+Backend detection, pinning, and a small wall-clock check. Headline A100
+numbers live in [benchmarks](https://jxucoder.github.io/openboost/benchmarks/).
 
-### GAM Explainability (`gam_explainability.py`)
+## GAM (`gam_explainability.py`)
 
-Interpretable machine learning with `OpenBoostGAM`.
-
-**Topics covered:**
-- Training interpretable GAM models
-- Visualizing shape functions
-- Per-feature contribution analysis
-- Explaining individual predictions
-- Trade-offs vs black-box models
+Main-effect shape functions, per-feature contributions.
 
 ## Requirements
 
-All examples work with the base OpenBoost installation:
-
 ```bash
-pip install openboost
-```
-
-Some examples benefit from optional dependencies:
-
-```bash
-# For sklearn integration examples
-pip install scikit-learn
-
-# For visualization
-pip install matplotlib
-
-# For GPU examples
-pip install numba  # CUDA support included
-```
-
-## Running Examples
-
-### Local Development
-
-```bash
-# From the repository root
-cd openboost
-
-# Run with uv
-uv run python examples/basic_regression.py
-
-# Or standard Python
-python examples/basic_regression.py
-```
-
-### In a Notebook
-
-```python
-# Copy-paste code from examples into Jupyter/Colab cells
-import openboost as ob
-
-model = ob.GradientBoosting(n_trees=100)
-model.fit(X_train, y_train)
-```
-
-### On Cloud (Modal, etc.)
-
-```python
-# Examples work on cloud GPU instances
-import modal
-
-app = modal.App()
-
-@app.function(gpu="A100")
-def train_model():
-    import openboost as ob
-    # ... example code ...
+pip install --pre openboost
+pip install scikit-learn matplotlib   # sklearn + plots
+pip install --pre "openboost[cuda]"   # GPU example
 ```
 
 ## Tips
 
-1. **Start simple**: Begin with `basic_regression.py` to understand the API
-2. **Check GPU**: Run `gpu_training.py` to verify GPU setup
-3. **Explore uncertainty**: `uncertainty_quantification.py` shows NaturalBoost's unique value
-4. **For Kaggle**: `kaggle_insurance.py` and `kaggle_sales.py` are ready-to-adapt templates
-5. **Custom needs**: `custom_loss.py` shows how to extend OpenBoost
+1. Start with `uncertainty_quantification.py` if you care about intervals.
+2. Confirm CUDA with `gpu_training.py` (`ob.is_cuda()`).
+3. Insurance / sales examples are the distribution-family templates.
+4. `custom_loss.py` is for point-estimate objectives; FormulaBoost is for
+   `y = f(θ, x)`.
 
 ## Troubleshooting
 
-**Example won't run?**
-- Ensure OpenBoost is installed: `pip install openboost`
-- For sklearn examples: `pip install scikit-learn`
-
-**GPU not detected?**
-- Check CUDA installation: `nvidia-smi`
-- Ensure numba is installed: `pip install numba`
-- See `gpu_training.py` for debugging tips
-
-**Plots not showing?**
-- Install matplotlib: `pip install matplotlib`
-- In headless environments, plots save to files
-
-## Contributing
+**Import errors.** `pip install --pre openboost` (without `--pre` you get
+the older stable release).
 
-Have a cool example to share? PRs welcome!
+**GPU not detected.** `nvidia-smi`, then
+`python -c "from numba import cuda; print(list(cuda.gpus))"`.
 
-Guidelines:
-- Self-contained (generates synthetic data or uses sklearn datasets)
-- Well-commented
-- Demonstrates a clear use case
-- Follows existing style
+**Headless plots.** Examples write figures to files when no display is
+available.
diff --git a/mkdocs.yml b/mkdocs.yml
index 11e3b6f..08c45df 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -1,5 +1,5 @@
 site_name: OpenBoost
-site_description: GPU-native, all-Python platform for tree-based machine learning
+site_description: GPU gradient boosting for distributional regression. Every parameter of F(y|x) gets its own tree ensemble, trained with the full natural gradient.
 site_url: https://jxucoder.github.io/openboost
 repo_url: https://github.com/jxucoder/openboost
 repo_name: jxucoder/openboost
@@ -81,16 +81,19 @@ nav:
       - Quickstart: getting-started/quickstart.md
       - GPU Setup: getting-started/gpu-setup.md
   - User Guide:
-      - Models:
+      - Distributional & varying-coefficient:
+          - How it works: user-guide/how-it-works.md
+          - NaturalBoost: user-guide/naturalboost/overview.md
+          - Distributions: user-guide/naturalboost/distributions.md
+          - Custom Distributions: user-guide/naturalboost/custom-distributions.md
+          - FormulaBoost: user-guide/formulaboost.md
+          - Weibull AFT: user-guide/survival.md
+      - Mean regression:
           - Gradient Boosting: user-guide/models/gradient-boosting.md
           - Multi-class Classification: user-guide/models/multiclass.md
           - DART: user-guide/models/dart.md
           - OpenBoostGAM: user-guide/models/gam.md
           - Linear Leaf GBDT: user-guide/models/linear-leaf.md
-      - NaturalBoost:
-          - Overview: user-guide/naturalboost/overview.md
-          - Distributions: user-guide/naturalboost/distributions.md
-          - Custom Distributions: user-guide/naturalboost/custom-distributions.md
       - Training:
           - Callbacks: user-guide/training/callbacks.md
           - Large-Scale Training: user-guide/training/large-scale.md
@@ -107,6 +110,7 @@ nav:
       - Custom Distribution: cookbook/custom-distribution.md
       - Custom Growth Strategy: cookbook/custom-growth-strategy.md
       - Device-Native Loss (GPU): cookbook/device-loss.md
+  - Benchmarks: benchmarks.md
   - API Reference:
       - openboost: api/openboost.md
       - Models: api/models.md
diff --git a/pyproject.toml b/pyproject.toml
index ae350e9..4666fac 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -5,14 +5,25 @@ build-backend = "hatchling.build"
 [project]
 name = "openboost"
 version = "1.0.0rc1"
-description = "Hackable gradient boosting platform: probabilistic predictions, interpretable GAMs, and custom algorithms in readable Python, with CPU and CUDA backends"
+description = "GPU gradient boosting for distributional regression and varying-coefficient models"
 readme = "README.md"
 license = "Apache-2.0"
 requires-python = ">=3.10"
 authors = [
     { name = "J. Xu" }
 ]
-keywords = ["gradient-boosting", "gbdt", "gpu", "cuda", "machine-learning"]
+keywords = [
+    "distributional-regression",
+    "varying-coefficient",
+    "gamlss",
+    "natural-gradient",
+    "probabilistic",
+    "survival",
+    "gbdt",
+    "gpu",
+    "cuda",
+    "machine-learning",
+]
 classifiers = [
     "Development Status :: 4 - Beta",
     "Intended Audience :: Developers",
diff --git a/src/openboost/__init__.py b/src/openboost/__init__.py
index 472bba7..8bfd74b 100644
--- a/src/openboost/__init__.py
+++ b/src/openboost/__init__.py
@@ -1,35 +1,21 @@
-"""OpenBoost: The PyTorch of Gradient Boosting.
+"""OpenBoost: GPU gradient boosting for distributional regression.
 
-Train-many optimized, research-friendly, GPU-accelerated gradient boosting.
+Every parameter of F(y|x) gets its own ensemble of histogram trees, updated
+with the full K x K natural gradient: GAMLSS-style distributional regression
+(NaturalBoost, plus WeibullAFT for right-censored data) and varying-coefficient
+models (FormulaBoost).
 
-Quick Start (Batched Training):
     >>> import openboost as ob
-    >>>
-    >>> # Simple scikit-learn-like API
-    >>> model = ob.GradientBoosting(n_trees=100, loss='mse')
+    >>> model = ob.NaturalBoostNormal(n_trees=100)
     >>> model.fit(X_train, y_train)
-    >>> predictions = model.predict(X_test)
+    >>> mean = model.predict(X_test)
+    >>> lo, hi = model.predict_interval(X_test, alpha=0.1)
 
-Custom Loss Functions:
-    >>> def quantile_loss(pred, y, tau=0.5):
-    ...     residual = y - pred
-    ...     grad = np.where(residual > 0, -tau, 1 - tau)
-    ...     hess = np.ones_like(pred)
-    ...     return grad, hess
-    >>> model = ob.GradientBoosting(n_trees=100, loss=quantile_loss)
-    >>> model.fit(X_train, y_train)
+    >>> model = ob.FormulaBoost(formula=f, n_params=2, links=("log", "identity"))
+    >>> model.fit(Z, y, model_input=x)
 
-Low-Level API (Full Control):
-    >>> # Bin data once, reuse everywhere
-    >>> X_binned = ob.array(X_train)
-    >>>
-    >>> # You own the training loop
-    >>> pred = np.zeros(len(y_train))
-    >>> for round in range(100):
-    ...     grad = 2 * (pred - y_train)  # Your loss, your gradients
-    ...     hess = np.ones_like(grad) * 2
-    ...     tree = ob.fit_tree(X_binned, grad, hess)
-    ...     pred = pred + 0.1 * tree(X_binned)
+    >>> model = ob.WeibullAFT()
+    >>> model.fit(Z, time, event=observed)
 """
 
 import warnings as _warnings
diff --git a/src/openboost/_models/_formula.py b/src/openboost/_models/_formula.py
index e017f8e..faf1ca5 100644
--- a/src/openboost/_models/_formula.py
+++ b/src/openboost/_models/_formula.py
@@ -1,6 +1,6 @@
 """FormulaBoost: boost every parameter of a user formula.
 
-Varying-coefficient / semi-parametric boosting. Features ``Z`` determine
+Varying-coefficient model. Features ``Z`` determine
 parameter surfaces ``theta(Z)`` via trees; a user formula
 ``y ≈ f(theta, x)`` consumes those parameters and a structural input ``x``.
 """
@@ -36,11 +36,13 @@ class FormulaBoost(PersistenceMixin):
             ``sigmoid``), length K.
         loss: Training loss. Currently ``mse`` only.
         precond: GGN preconditioner: ``full`` (default), ``diag``, or
-            ``plain`` (raw gradient — usually a bad idea).
+            ``plain`` (raw gradient, usually a bad idea).
         damp: Levenberg–Marquardt damping added to the GGN matrix.
         param_names: Optional names for the K parameters. Defaults to
             ``theta_0``, ``theta_1``, ...
-        n_trees, max_depth, learning_rate, ...: standard tree knobs.
+
+    Tree knobs (``n_trees``, ``max_depth``, ``learning_rate``, ...) match
+    ``GradientBoosting``.
 
     Example:
         ```python
diff --git a/src/openboost/_models/_survival.py b/src/openboost/_models/_survival.py
index 00c9fdf..ca29c81 100644
--- a/src/openboost/_models/_survival.py
+++ b/src/openboost/_models/_survival.py
@@ -31,7 +31,9 @@ class WeibullAFT(PersistenceMixin):
 
     Args:
         damp: Levenberg-Marquardt damping on the per-sample 2x2 information.
-        n_trees, max_depth, learning_rate, ...: standard tree knobs.
+
+    Tree knobs (``n_trees``, ``max_depth``, ``learning_rate``, ...) match
+    ``GradientBoosting``.
 
     Example:
         ```python