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