diff --git a/README.md b/README.md
index e6cb7af5..16773d07 100644
--- a/README.md
+++ b/README.md
@@ -7,15 +7,21 @@
 [![Docs](https://img.shields.io/badge/docs-ovvo--financial.github.io-blue)](https://ovvo-financial.github.io/NNS-python/)
 [![License](https://img.shields.io/badge/license-GPL--3.0--only-blue)](https://github.com/OVVO-Financial/NNS-python/blob/main/LICENSE)
 
-`ovvo-nns` brings Nonlinear Nonparametric Statistics to Python as the `nns` import package. It is a parity-focused port of the R `NNS` 13.0+ package, designed for real-world data that violate symmetry, linearity, or distributional assumptions.
-
-NNS is built around partial moments, the lower and upper components of variance, and uses them across nonlinear dependence, correlation, causation, regression, classification, forecasting, stochastic dominance, stochastic superiority, Monte Carlo simulation, and numerical differentiation workflows.
-
-
-> NNS was created by Fred Viole as the companion R package to Viole, F. and Nawrocki, D. (2013), *Nonlinear Nonparametric Statistics: Using Partial Moments*. **Book (2nd Edition):** https://ovvo-financial.github.io/NNS/book/
->   
->  **Implementation:** For a direct quantitative finance implementation of NNS, see [OVVO Labs](https://www.ovvolabs.com)
-
+`ovvo-nns` brings Nonlinear Nonparametric Statistics to Python as the `nns`
+import package. It is a parity-focused port of the R `NNS` 13.0+ package,
+designed for real-world data that violate symmetry, linearity, or
+distributional assumptions.
+
+The R package is the reference implementation and the source of truth for the
+statistical behavior, terminology, and canonical example curriculum. Python is
+native and does not call R at runtime.
+
+> NNS was created by Fred Viole as the companion R package to Viole, F. and
+> Nawrocki, D. (2013), *Nonlinear Nonparametric Statistics: Using Partial
+> Moments*. **Book (2nd Edition):** https://ovvo-financial.github.io/NNS/book/
+>
+> **Implementation:** For a direct quantitative finance implementation of NNS,
+> see [OVVO Labs](https://www.ovvolabs.com).
 
 ## Package at a glance
 
@@ -25,24 +31,18 @@ NNS is built around partial moments, the lower and upper components of variance,
 | Import package | `nns` |
 | Current version | `1.4.0` |
 | Python | `>=3.11` |
-| Required runtime dependencies | NumPy, SciPy |
+| Required runtime dependencies | NumPy, SciPy, Matplotlib |
 | R required at runtime | No |
 | Native acceleration | Private, optional `nns._nnscore` kernels where available |
 | Public API status | Stable, parity-focused |
 | License | GPL-3.0-only |
 
-The public package is Python-native and does not call R at runtime. Some core kernels can use the private `_nnscore` extension when it is present, while public functions keep Python implementations and explicit fallback behavior.
-
 ## Install
 
 ```bash
 pip install ovvo-nns
 ```
 
-This includes the matplotlib plotting API (`nns.plotting`); matplotlib is a
-regular dependency and is imported lazily, so `import nns` stays light. See
-the [plot parity policy](https://ovvo-financial.github.io/NNS-python/plot_parity_policy/).
-
 Use the package as `nns`:
 
 ```python
@@ -51,7 +51,8 @@ import nns
 print(nns.__version__)
 ```
 
-Source builds use `scikit-build-core` and `nanobind` for the optional native extension. Published wheels should be preferred when available.
+Published wheels are preferred. Source builds use `scikit-build-core` and
+`nanobind` for the optional native extension.
 
 ## Quick start
 
@@ -59,58 +60,15 @@ Source builds use `scikit-build-core` and `nanobind` for the optional native ext
 import numpy as np
 from nns import lpm, nns_dep, nns_reg, upm
 
-x = np.array([-2.0, -1.0, 0.5, 3.0], dtype=np.float64)
-
-lower = lpm(degree=2, target=0.0, x=x)
-upper = upm(degree=2, target=0.0, x=x)
-
-print("lower partial moment:", lower)
-print("upper partial moment:", upper)
-```
-
-Measure nonlinear dependence:
-
-```python
-import numpy as np
-from nns import nns_cor, nns_dep
-
-grid = np.linspace(-2.0, 2.0, 80, dtype=np.float64)
-y = grid**2
-
-print("NNS dependence:", nns_dep(grid, y))
-print("NNS correlation:", nns_cor(grid, y))
-```
-
-Fit a nonlinear regression and estimate new points:
-
-```python
-import numpy as np
-from nns import nns_reg
-
-x = np.linspace(-3.0, 3.0, 80, dtype=np.float64)
-y = np.sin(x) + 0.2 * x
-points = np.array([-1.5, 0.0, 1.5], dtype=np.float64)
+x = np.array([-2.0, -1.0, 0.5, 3.0])
+print("LPM2:", lpm(2, 0.0, x))
+print("UPM2:", upm(2, 0.0, x))
 
-fit = nns_reg(x, y, point_est=points, confidence_interval=None)
+grid = np.linspace(-2.0, 2.0, 80)
+print("nonlinear dependence:", nns_dep(grid, grid**2))
 
-print("R2:", fit["R2"])
-print(np.column_stack((points, fit["Point.est"])))
-```
-
-Forecast a univariate series:
-
-```python
-import numpy as np
-from nns import nns_arma, nns_seas
-
-t = np.arange(1, 60, dtype=np.float64)
-series = 10.0 + np.sin(t / 3.0) + 0.05 * t
-
-seasonality = nns_seas(series, modulo=[3, 4, 6], mod_only=True)
-forecast = nns_arma(series, h=3, seasonal_factor=4, method="lin")
-
-print("best seasonal period:", seasonality["best.period"])
-print("forecast:", forecast)
+fit = nns_reg(grid, np.sin(grid), point_est=np.array([-1.0, 0.0, 1.0]))
+print("point estimates:", fit["Point.est"])
 ```
 
 ## Main API areas
@@ -127,67 +85,69 @@ print("forecast:", forecast)
 | Stochastic dominance | `fsd`, `ssd`, `tsd`, `nns_sd_cluster`, `sd_efficient_set` |
 | Stochastic superiority and simulation | `nns_ss`, `nns_mc`, `nns_meboot` |
 | Differentiation | `nns_diff`, `dy_dx`, `dy_d` |
-| Categorical helpers | `encode_factor_codes`, `factor_2_dummy`, `factor_2_dummy_fr`, `prepare_factor_predictors` |
-
-See [API status](https://ovvo-financial.github.io/NNS-python/api_status/) for implemented, partial, guarded, and known-gap paths.
-
-## Design boundaries
 
-NNS Python prioritizes stable public behavior from installed R NNS 13.0+, not private helper parity. The package returns NumPy arrays and plain dictionaries rather than R `data.table` objects, uses explicit Python errors for several unsafe R coercions, and generally ignores plotting side effects.
+See [API status](https://ovvo-financial.github.io/NNS-python/api_status/) for
+implemented, partial, guarded, and known-gap paths.
 
-Important boundaries:
+## Canonical examples
 
-- R is used only for parity tests and local cache regeneration, not normal runtime use.
-- Stochastic exact stream parity is not expected because Python paths use NumPy random generation.
-- Factor and class ordering should be passed explicitly when ordering matters.
-- Direct raw-factor `nns_m_reg(..., factor_2_dummy=True)` is intentionally guarded. Use `prepare_factor_predictors(...)` before `nns_m_reg(...)`.
-- Compute functions still return values, not figures; passing `plot=True` (where R has it) additionally renders a Matplotlib figure as a side effect via the `nns.plotting` layer, which is color/element-faithful to R but not pixel-diffed. The plot functions can also be called directly on a computed result.
+The R package's nine numbered vignettes define the canonical NNS curriculum.
+Python follows the same numbering, topic names, and statistical intent:
 
-See [behavior conventions](https://ovvo-financial.github.io/NNS-python/conventions/) for detailed compatibility notes.
+| # | Topic | Python entry point |
+|---|---|---|
+| 01 | Overview | [`01_overview.py`](examples/vignettes/01_overview.py) |
+| 02 | Partial Moments | [`02_partial_moments.py`](examples/vignettes/02_partial_moments.py) |
+| 03 | Correlation and Dependence | [`03_correlation_and_dependence.py`](examples/vignettes/03_correlation_and_dependence.py) |
+| 04 | Normalization and Rescaling | [`04_normalization_and_rescaling.py`](examples/vignettes/04_normalization_and_rescaling.py) |
+| 05 | Sampling and Simulation | [`05_sampling_and_simulation.py`](examples/vignettes/05_sampling_and_simulation.py) |
+| 06 | Comparing Distributions | [`06_comparing_distributions.py`](examples/vignettes/06_comparing_distributions.py) |
+| 07 | Clustering and Regression | [`07_clustering_and_regression.py`](examples/vignettes/07_clustering_and_regression.py) |
+| 08 | Classification | [`08_classification.py`](examples/vignettes/08_classification.py) |
+| 09 | Forecasting | [`09_forecasting.py`](examples/vignettes/09_forecasting.py) |
 
-## Examples
+The mapping is recorded in
+[`examples/vignettes/manifest.yml`](examples/vignettes/manifest.yml) and
+validated in CI. The older unnumbered scripts remain as focused examples and
+backward-compatible entry points.
 
-Runnable, self-checking example scripts live in
-[`examples/vignettes`](https://github.com/OVVO-Financial/NNS-python/tree/main/examples/vignettes), mirroring the R NNS vignettes. They
-are exercised in CI by `tests/docs/test_vignette_examples.py`, so they stay in
-sync with the package.
-
-| Topic | Script |
-|---|---|
-| Overview | [`overview.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/overview.py) |
-| Partial moments | [`partial_moments.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/partial_moments.py) |
-| Descriptive and distributional tools | [`descriptive_distributional_tools.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/descriptive_distributional_tools.py) |
-| Dependence and nonlinear association | [`dependence_nonlinear_association.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/dependence_nonlinear_association.py) |
-| Normalization and rescaling | [`normalization_rescaling.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/normalization_rescaling.py) |
-| Hypothesis, ANOVA and stochastic superiority | [`hypothesis_anova_stochastic_superiority.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/hypothesis_anova_stochastic_superiority.py) |
-| Regression, boosting, stacking and causality | [`regression_boosting_stacking_causality.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/regression_boosting_stacking_causality.py) |
-| Time series forecasting | [`time_series_forecasting.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/time_series_forecasting.py) |
-| Simulation, bootstrap and risk-neutral | [`simulation_bootstrap_riskneutral.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/simulation_bootstrap_riskneutral.py) |
-| Portfolio and stochastic dominance | [`portfolio_stochastic_dominance.py`](https://github.com/OVVO-Financial/NNS-python/blob/main/examples/vignettes/portfolio_stochastic_dominance.py) |
-
-Run one example:
+Run one canonical example:
 
 ```bash
-uv run python examples/vignettes/partial_moments.py
+uv run python examples/vignettes/02_partial_moments.py
 ```
 
-Run all of them with a PASS/FAIL summary:
+Run all nine in R curriculum order:
 
 ```bash
 uv run python examples/run_all_vignettes.py
 ```
 
-## Documentation
+## Design boundaries
+
+NNS Python prioritizes stable public behavior from installed R NNS 13.0+, not
+private helper parity. The package returns NumPy arrays and plain dictionaries
+rather than R `data.table` objects and uses explicit Python errors for unsafe R
+coercions.
 
-The full documentation site is hosted at
-**<https://ovvo-financial.github.io/NNS-python/>**.
+- R is used for parity tests and local cache regeneration, not normal runtime.
+- Exact stochastic stream parity is not expected for every randomized path.
+- Factor and class ordering should be supplied explicitly when it matters.
+- Classification codes follow the R contract and start at 1.
+- Compute functions return values; `plot=True` adds Matplotlib rendering as a
+  side effect without changing the statistical result.
 
-- [API reference manual](https://ovvo-financial.github.io/NNS-python/api_reference/)
+See [behavior conventions](https://ovvo-financial.github.io/NNS-python/conventions/)
+for detailed compatibility notes.
+
+## Documentation
+
+- [API reference](https://ovvo-financial.github.io/NNS-python/api_reference/)
 - [API status and known gaps](https://ovvo-financial.github.io/NNS-python/api_status/)
-- [Behavior conventions and intentional divergences](https://ovvo-financial.github.io/NNS-python/conventions/)
-- [Parity target, cache regeneration, and automation](https://ovvo-financial.github.io/NNS-python/parity/)
+- [Behavior conventions](https://ovvo-financial.github.io/NNS-python/conventions/)
+- [Parity policy and cache regeneration](https://ovvo-financial.github.io/NNS-python/parity/)
 - [Benchmarks](https://ovvo-financial.github.io/NNS-python/benchmarks/)
-- [Examples](https://github.com/OVVO-Financial/NNS-python/tree/main/examples/vignettes)
+- [Canonical examples](examples/vignettes/README.md)
 
 ## Development
 
@@ -204,11 +164,9 @@ Run benchmark tests explicitly:
 uv run pytest -n0 -m benchmark --benchmark-enable tests/benchmarks/
 ```
 
-The default parity suite is cache-backed and does not require `Rscript`. `Rscript` and the R `NNS` package are needed only when regenerating parity caches or running live R comparison scripts.
-
-## Benchmarks
-
-Benchmarks compare selected Python paths with installed R NNS 13.0+ baselines. Many core operations are faster in Python, while some large stochastic-dominance workloads remain faster in R because the R package uses compiled kernels for those paths. See [benchmarks](https://ovvo-financial.github.io/NNS-python/benchmarks/) for current measurements and commands.
+The default parity suite is cache-backed and does not require `Rscript`.
+`Rscript` and the R `NNS` package are needed only when regenerating parity
+caches or running live R comparison scripts.
 
 ## Authors and contributors
 
@@ -218,4 +176,5 @@ Benchmarks compare selected Python paths with installed R NNS 13.0+ baselines. M
 
 ## Attribution
 
-Upstream R package and reference implementation: [OVVO-Financial/NNS](https://github.com/OVVO-Financial/NNS)
+Upstream R package and reference implementation:
+[OVVO-Financial/NNS](https://github.com/OVVO-Financial/NNS)
diff --git a/examples/run_all_vignettes.py b/examples/run_all_vignettes.py
index 1599c616..852dd67b 100644
--- a/examples/run_all_vignettes.py
+++ b/examples/run_all_vignettes.py
@@ -1,21 +1,14 @@
-"""Run every NNS Python vignette end to end and print its output.
+"""Run the nine canonical NNS Python vignettes in R curriculum order.
 
-This is a single, IDLE-friendly driver for the vignette example scripts in
-``examples/vignettes/``. Open it in IDLE and press **F5** (or run
-``python examples/run_all_vignettes.py`` from a terminal) to execute all
-vignettes in the documented order and print each one's output, so you can
-compare it against the R NNS vignettes PDF.
-
-Each vignette also self-checks with assertions, so this driver reports PASS/FAIL
-per vignette and a final summary, and exits non-zero if any vignette fails.
-
-Requires the package to be importable (``pip install -e .`` from the repo root).
+The R package is the source of truth. The numbered Python entry points under
+``examples/vignettes`` follow the same 01-09 topic sequence and self-check with
+assertions. This driver reports PASS/FAIL and exits non-zero on failure.
 """
-
 from __future__ import annotations
 
 import importlib.util
 import os
+import sys
 import time
 import traceback
 from pathlib import Path
@@ -23,32 +16,20 @@
 REPO_ROOT = Path(__file__).resolve().parents[1]
 VIGNETTE_DIR = REPO_ROOT / "examples" / "vignettes"
 
-# Ordered to match the R NNS vignettes PDF.
-# (number, title, example-script stem)
 VIGNETTES = [
-    ("00", "Overview", "overview"),
-    ("01", "Partial moments", "partial_moments"),
-    ("02", "Descriptive & distributional tools",
-     "descriptive_distributional_tools"),
-    ("03", "Dependence & nonlinear association",
-     "dependence_nonlinear_association"),
-    ("04", "Normalization & rescaling",
-     "normalization_rescaling"),
-    ("05", "Hypothesis, ANOVA & stochastic superiority",
-     "hypothesis_anova_stochastic_superiority"),
-    ("06", "Regression, boosting, stacking & causality",
-     "regression_boosting_stacking_causality"),
-    ("07", "Time series forecasting",
-     "time_series_forecasting"),
-    ("08", "Simulation, bootstrap & risk-neutral",
-     "simulation_bootstrap_riskneutral"),
-    ("09", "Portfolio & stochastic dominance",
-     "portfolio_stochastic_dominance"),
+    ("01", "Overview", "01_overview"),
+    ("02", "Partial Moments", "02_partial_moments"),
+    ("03", "Correlation and Dependence", "03_correlation_and_dependence"),
+    ("04", "Normalization and Rescaling", "04_normalization_and_rescaling"),
+    ("05", "Sampling and Simulation", "05_sampling_and_simulation"),
+    ("06", "Comparing Distributions", "06_comparing_distributions"),
+    ("07", "Clustering and Regression", "07_clustering_and_regression"),
+    ("08", "Classification", "08_classification"),
+    ("09", "Forecasting", "09_forecasting"),
 ]
 
 
-def _load_vignette_main(stem):
-    """Import an example script by path and return its ``main`` callable."""
+def _load_vignette_main(stem: str):
     path = VIGNETTE_DIR / f"{stem}.py"
     spec = importlib.util.spec_from_file_location(f"nns_vignette_{stem}", path)
     if spec is None or spec.loader is None:
@@ -59,10 +40,12 @@ def _load_vignette_main(stem):
 
 
 def run() -> int:
-    # Match the cwd the test suite uses so any relative paths resolve.
     os.chdir(REPO_ROOT)
+    # Numbered compatibility entry points import maintained sibling scripts.
+    if str(VIGNETTE_DIR) not in sys.path:
+        sys.path.insert(0, str(VIGNETTE_DIR))
 
-    results = []
+    results: list[tuple[str, str, bool, float]] = []
     for number, title, stem in VIGNETTES:
         banner = f"  Vignette {number}: {title}  "
         rule = "=" * max(len(banner), 60)
@@ -78,11 +61,10 @@ def run() -> int:
         except Exception:
             ok = False
             print(traceback.format_exc())
-        elapsed = time.perf_counter() - start
-        results.append((number, title, ok, elapsed))
+        results.append((number, title, ok, time.perf_counter() - start))
 
     print("\n" + "=" * 60)
-    print("  Vignette verification summary")
+    print("  Canonical vignette verification summary")
     print("=" * 60)
     passed = 0
     for number, title, ok, elapsed in results:
@@ -91,8 +73,6 @@ def run() -> int:
         print(f"  [{status}] {number}  {title}  ({elapsed:.2f}s)")
     print("-" * 60)
     print(f"  {passed}/{len(results)} vignettes passed")
-    if passed != len(results):
-        print("  Some vignettes FAILED — see tracebacks above.")
     return 0 if passed == len(results) else 1
 
 
diff --git a/examples/vignettes/01_overview.py b/examples/vignettes/01_overview.py
new file mode 100644
index 00000000..55de7a7a
--- /dev/null
+++ b/examples/vignettes/01_overview.py
@@ -0,0 +1,12 @@
+"""Canonical vignette 01: Overview.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_01_Overview.Rmd
+
+The maintained implementation lives in ``overview.py``. This numbered entry
+point keeps the Python example catalog aligned with the canonical R sequence.
+"""
+from overview import main
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/02_partial_moments.py b/examples/vignettes/02_partial_moments.py
new file mode 100644
index 00000000..18a0af10
--- /dev/null
+++ b/examples/vignettes/02_partial_moments.py
@@ -0,0 +1,80 @@
+"""Canonical vignette 02: Partial moments.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_02_Partial_Moments.Rmd
+
+Covers the mean and variance identities, covariance reconstruction, empirical
+CDFs, a joint CDF/Bayes identity, inverse-CDF sampling, and numerical
+integration using the public Python API.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from nns import (
+    co_lpm,
+    co_upm,
+    d_lpm,
+    d_upm,
+    lpm,
+    lpm_ratio,
+    lpm_var,
+    nns_mode,
+    nns_moments,
+    pm_matrix,
+    upm,
+)
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    x = rng.normal(size=100)
+    y = rng.normal(size=100)
+    n = x.size
+
+    mean_pm = float(upm(1, 0.0, x) - lpm(1, 0.0, x))
+    variance_pm = float(upm(2, x.mean(), x) + lpm(2, x.mean(), x)) * n / (n - 1)
+    assert np.isclose(mean_pm, x.mean())
+    assert np.isclose(variance_pm, x.var(ddof=1))
+
+    covariance_pm = (
+        co_lpm(1, x, y, x.mean(), y.mean())
+        + co_upm(1, x, y, x.mean(), y.mean())
+        - d_lpm(1, 1, x, y, x.mean(), y.mean())
+        - d_upm(1, 1, x, y, x.mean(), y.mean())
+    ) * n / (n - 1)
+    assert np.isclose(float(covariance_pm), np.cov(x, y)[0, 1])
+
+    pm = pm_matrix(1, 1, "mean", np.column_stack((x, y)), True, names=["x", "y"])
+    reconstructed = pm["clpm"] + pm["cupm"] - pm["dlpm"] - pm["dupm"]
+    np.testing.assert_allclose(reconstructed, np.cov(x, y), atol=1e-8)
+
+    targets = np.array([0.0, 1.0])
+    cdf_pm = np.array([float(lpm(0, t, x)) for t in targets])
+    np.testing.assert_allclose(cdf_pm, [np.mean(x <= t) for t in targets])
+
+    # P(X > 0 | Y > 0) = P(X > 0, Y > 0) / P(Y > 0).
+    bayes = float(co_upm(0, x, y, 0.0, 0.0) / upm(0, 0.0, y))
+    assert 0.0 <= bayes <= 1.0
+
+    percentiles = np.linspace(0.05, 0.95, 19)
+    samples = np.asarray([lpm_var(p, 0.0, x) for p in percentiles], dtype=float)
+    recovered = np.asarray([lpm_ratio(0, t, x) for t in samples], dtype=float)
+    assert samples.shape == recovered.shape
+
+    grid = np.linspace(0.0, 1.0, 1001)
+    integral = float(upm(1, 0.0, grid**2) - lpm(1, 0.0, grid**2))
+    assert np.isclose(integral, 1.0 / 3.0, atol=2e-3)
+
+    print("mean via partial moments:", round(mean_pm, 6))
+    print("sample variance via partial moments:", round(variance_pm, 6))
+    print("sample covariance via co-partial moments:", round(float(covariance_pm), 6))
+    print("moments:", nns_moments(x))
+    print("mode:", nns_mode(x))
+    print("CDF at [0, 1]:", np.round(cdf_pm, 4))
+    print("P(X > 0 | Y > 0):", round(bayes, 4))
+    print("integral of x^2 on [0, 1]:", round(integral, 6))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/03_correlation_and_dependence.py b/examples/vignettes/03_correlation_and_dependence.py
new file mode 100644
index 00000000..499d8c99
--- /dev/null
+++ b/examples/vignettes/03_correlation_and_dependence.py
@@ -0,0 +1,48 @@
+"""Canonical vignette 03: Correlation and dependence.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_03_Correlation_and_Dependence.Rmd
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_copula, nns_dep
+
+
+def main() -> None:
+    x = np.arange(0.0, 3.01, 0.01)
+    linear = nns_dep(x, 2.0 * x)
+    nonlinear = nns_dep(x, x**10)
+
+    cyclic_x = np.arange(0.0, 12.0 * np.pi, np.pi / 100.0)
+    cyclic_y = np.sin(cyclic_x)
+    cyclic = nns_dep(cyclic_x, cyclic_y)
+    asym_xy = nns_dep(cyclic_x, cyclic_y, asym=True)["Dependence"]
+    asym_yx = nns_dep(cyclic_y, cyclic_x, asym=True)["Dependence"]
+
+    rng = np.random.default_rng(123)
+    points = rng.uniform(-1.0, 1.0, size=(10000, 2))
+    radius2 = np.sum(points**2, axis=1)
+    ring = points[(radius2 <= 1.0) & (radius2 >= 0.95)]
+    ring_dep = nns_dep(ring[:, 0], ring[:, 1])
+
+    frame = rng.normal(size=(1000, 3))
+    copula = float(nns_copula(frame, continuous=True))
+
+    assert linear["Correlation"] > 0.99 and linear["Dependence"] > 0.99
+    assert nonlinear["Dependence"] >= abs(nonlinear["Correlation"])
+    assert cyclic["Dependence"] > abs(cyclic["Correlation"])
+    assert 0.0 <= float(ring_dep["Dependence"]) <= 1.0
+    assert 0.0 <= copula <= 1.0
+
+    print("linear:", linear)
+    print("x^10:", nonlinear)
+    print("sin(x):", cyclic)
+    print("asymmetric D(y|x), D(x|y):", round(float(asym_xy), 4), round(float(asym_yx), 4))
+    print("ring dependence:", round(float(ring_dep["Dependence"]), 4))
+    print("three-variable copula dependence:", round(copula, 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/04_normalization_and_rescaling.py b/examples/vignettes/04_normalization_and_rescaling.py
new file mode 100644
index 00000000..4b25b3ef
--- /dev/null
+++ b/examples/vignettes/04_normalization_and_rescaling.py
@@ -0,0 +1,11 @@
+"""Canonical vignette 04: Normalization and rescaling.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_04_Normalization_and_Rescaling.Rmd
+
+The maintained implementation lives in ``normalization_rescaling.py``.
+"""
+from normalization_rescaling import main
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/05_sampling_and_simulation.py b/examples/vignettes/05_sampling_and_simulation.py
new file mode 100644
index 00000000..019ade5f
--- /dev/null
+++ b/examples/vignettes/05_sampling_and_simulation.py
@@ -0,0 +1,47 @@
+"""Canonical vignette 05: Sampling and simulation.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_05_Sampling.Rmd
+
+Combines partial-moment CDF/inverse-CDF sampling with maximum-entropy bootstrap
+and dependence-targeted Monte Carlo examples.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from nns import lpm_ratio, lpm_var, nns_mc, nns_meboot
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    x = rng.normal(size=100)
+
+    targets = np.sort(x)
+    empirical = np.asarray([np.mean(x <= t) for t in targets])
+    pm_cdf = np.asarray([lpm_ratio(0, t, x) for t in targets], dtype=float)
+    np.testing.assert_allclose(pm_cdf, empirical)
+
+    percentiles = np.linspace(0.01, 0.99, 99)
+    samples = {
+        degree: np.asarray([lpm_var(p, degree, x) for p in percentiles], dtype=float)
+        for degree in (0.0, 0.25, 0.5, 1.0, 2.0)
+    }
+    np.testing.assert_allclose(samples[0.0], np.quantile(x, percentiles, method="linear"))
+
+    series = np.cumsum(rng.normal(scale=0.7, size=80))
+    meboot = nns_meboot(series, reps=10, rho=0.95, random_seed=1)
+    mc = nns_mc(series, reps=1, lower_rho=-1.0, upper_rho=1.0, by=0.5, random_seed=1)
+    assert np.asarray(meboot["ensemble"]).shape == series.shape
+    assert np.asarray(mc["ensemble"]).shape == series.shape
+
+    print("CDF parity max error:", float(np.max(np.abs(pm_cdf - empirical))))
+    print("inverse-CDF sample heads:")
+    for degree, values in samples.items():
+        print(f"  degree={degree:g}:", np.round(values[:5], 4))
+    print("meboot ensemble shape:", np.asarray(meboot["ensemble"]).shape)
+    print("MC rho groups:", list(mc["replicates"]))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/06_comparing_distributions.py b/examples/vignettes/06_comparing_distributions.py
new file mode 100644
index 00000000..0d864a83
--- /dev/null
+++ b/examples/vignettes/06_comparing_distributions.py
@@ -0,0 +1,62 @@
+"""Canonical vignette 06: Comparing distributions.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_06_Comparing_Distributions.Rmd
+
+Covers NNS ANOVA certainty, stochastic superiority, pairwise stochastic
+dominance, the stochastic-dominance efficient set, and SD clustering.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from nns import (
+    fsd_uni,
+    nns_anova,
+    nns_sd_cluster,
+    nns_ss,
+    sd_efficient_set,
+    ssd_uni,
+    tsd_uni,
+)
+
+
+def main() -> None:
+    rng = np.random.default_rng(123)
+    x = rng.normal(0.0, 1.0, size=1000)
+    equal_mean = rng.normal(0.0, 2.0, size=1000)
+    shifted = x + 1.0
+
+    equal = nns_anova(x, equal_mean, means_only=True, random_seed=1)
+    unequal = nns_anova(x, shifted, means_only=True, random_seed=1)
+    superiority = nns_ss(x, shifted)
+
+    assert 0.0 <= float(equal["Certainty"]) <= 1.0
+    assert 0.0 <= float(unequal["Certainty"]) <= 1.0
+    assert fsd_uni(shifted, x) == 1
+    assert ssd_uni(shifted, x) == 1
+    assert tsd_uni(shifted, x) == 1
+
+    base = [rng.normal(size=500) for _ in range(4)]
+    panel = np.column_stack(
+        [item for pair in ((z, z + 1.0) for z in base) for item in pair]
+    )
+    efficient = sd_efficient_set(panel, degree=1)
+    names = [f"x{i + 1}" for i in range(panel.shape[1])]
+    clusters = nns_sd_cluster(panel, degree=1, names=names)
+
+    print("ANOVA certainty, equal means:", round(float(equal["Certainty"]), 4))
+    print("ANOVA certainty, shifted means:", round(float(unequal["Certainty"]), 4))
+    print("stochastic superiority:", superiority)
+    print(
+        "FSD/SSD/TSD shifted over base:",
+        fsd_uni(shifted, x),
+        ssd_uni(shifted, x),
+        tsd_uni(shifted, x),
+    )
+    print("SD efficient set:", efficient)
+    print("SD clusters:", clusters["Clusters"])
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/07_clustering_and_regression.py b/examples/vignettes/07_clustering_and_regression.py
new file mode 100644
index 00000000..be465a8c
--- /dev/null
+++ b/examples/vignettes/07_clustering_and_regression.py
@@ -0,0 +1,45 @@
+"""Canonical vignette 07: Clustering and regression.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_07_Clustering_and_Regression.Rmd
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_part, nns_reg
+
+
+def main() -> None:
+    x = np.arange(-5.0, 5.05, 0.05)
+    y = x**3
+
+    full = nns_part(x, y, order=4, obs_req=0)
+    x_only = nns_part(x, y, type="XONLY", order=4, obs_req=0)
+    assert full["order"] == 4 and x_only["order"] == 4
+
+    # X-only partition paths may contain multiple levels (for example q1121),
+    # but every branch after the root must be a left/right label: 1 or 2.
+    quadrant_paths = np.asarray(x_only["dt"]["quadrant"], dtype=str)
+    assert all(set(path.removeprefix("q")) <= {"1", "2"} for path in quadrant_paths)
+
+    points = np.array([-6.0, -2.0, 0.0, 2.0, 6.0])
+    univariate = nns_reg(x, y, point_est=points, confidence_interval=None)
+    assert np.asarray(univariate["Point.est"]).shape == points.shape
+
+    rng = np.random.default_rng(123)
+    design = rng.uniform(-2.0, 2.0, size=(160, 2))
+    target = design[:, 0] ** 3 + 3.0 * design[:, 1] - design[:, 1] ** 3 - 3.0 * design[:, 0]
+    multivariate = nns_reg(design, target, point_est=design[:10], order="max")
+    smoothed = nns_reg(x, y + rng.normal(scale=2.0, size=x.size), point_est=points, smooth=True)
+
+    print("joint partition regression points:", len(full["regression.points"]["x"]))
+    print("X-only partition regression points:", len(x_only["regression.points"]["x"]))
+    print("univariate R2:", round(float(univariate["R2"]), 4))
+    print("univariate point estimates:", np.round(univariate["Point.est"], 4))
+    print("multivariate point estimates:", np.round(multivariate["Point.est"], 4))
+    print("smoothed point estimates:", np.round(smoothed["Point.est"], 4))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/08_classification.py b/examples/vignettes/08_classification.py
new file mode 100644
index 00000000..3d1f3e07
--- /dev/null
+++ b/examples/vignettes/08_classification.py
@@ -0,0 +1,73 @@
+"""Canonical vignette 08: Classification.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_08_Classification.Rmd
+
+Demonstrates the same three public classification paths as R: ``nns_reg`` as
+the base learner, ``nns_boost`` as the resampled ensemble, and ``nns_stack``
+as the cross-validated regression/dimension-reduction ensemble. Class codes
+start at 1, matching the R contract.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from nns import nns_boost, nns_reg, nns_stack
+
+
+def _three_class_data(seed: int = 123) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    rng = np.random.default_rng(seed)
+    centers = np.array([[-3.0, -2.0, 0.0, 1.0], [0.0, 2.5, 3.0, -1.0], [3.0, -1.0, -2.5, 2.0]])
+    train = np.vstack([rng.normal(center, 0.45, size=(35, 4)) for center in centers])
+    test = np.vstack([rng.normal(center, 0.45, size=(5, 4)) for center in centers])
+    y_train = np.repeat(np.arange(1, 4), 35).astype(float)
+    y_test = np.repeat(np.arange(1, 4), 5).astype(float)
+    return train, y_train, test, y_test
+
+
+def main() -> None:
+    x_train, y_train, x_test, y_test = _three_class_data()
+
+    reg = nns_reg(x_train, y_train, type="CLASS", point_est=x_test)
+    boost = nns_boost(
+        x_train,
+        y_train,
+        x_test,
+        type="CLASS",
+        epochs=5,
+        learner_trials=10,
+        cv_size=0.25,
+        balance=True,
+        status=False,
+        seed=123,
+    )
+    stack = nns_stack(
+        x_train,
+        y_train,
+        x_test,
+        type="CLASS",
+        balance=True,
+        folds=1,
+        cv_size=0.25,
+        status=False,
+        seed=123,
+    )
+
+    predictions = {
+        "reg": np.asarray(reg["Point.est"], dtype=float),
+        "boost": np.asarray(boost["results"], dtype=float),
+        "stack": np.asarray(stack["stack"], dtype=float),
+    }
+    for values in predictions.values():
+        assert values.shape == y_test.shape
+        assert set(np.unique(values)) <= {1.0, 2.0, 3.0}
+
+    for name, values in predictions.items():
+        print(f"{name} accuracy:", round(float(np.mean(values == y_test)), 4))
+        print(f"{name} predictions:", values.astype(int))
+    print("boost feature weights:", boost["feature.weights"])
+    print("stack weights:", stack["weights"])
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/09_forecasting.py b/examples/vignettes/09_forecasting.py
new file mode 100644
index 00000000..c36c155a
--- /dev/null
+++ b/examples/vignettes/09_forecasting.py
@@ -0,0 +1,11 @@
+"""Canonical vignette 09: Forecasting.
+
+Source of truth:
+    NNS/vignettes/NNSvignette_09_Forecasting.Rmd
+
+The maintained implementation lives in ``time_series_forecasting.py``.
+"""
+from time_series_forecasting import main
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/vignettes/README.md b/examples/vignettes/README.md
new file mode 100644
index 00000000..f5f0a1ac
--- /dev/null
+++ b/examples/vignettes/README.md
@@ -0,0 +1,34 @@
+# Canonical NNS examples
+
+The R package is the source of truth for the NNS example curriculum. Its nine
+numbered vignettes define the statistical narrative, datasets, section order,
+and interpretation. Python implements the same public workflows with
+language-appropriate containers and syntax.
+
+| # | Canonical topic | R source | Python entry point |
+|---|---|---|---|
+| 01 | Overview | `NNSvignette_01_Overview.Rmd` | `01_overview.py` |
+| 02 | Partial Moments | `NNSvignette_02_Partial_Moments.Rmd` | `02_partial_moments.py` |
+| 03 | Correlation and Dependence | `NNSvignette_03_Correlation_and_Dependence.Rmd` | `03_correlation_and_dependence.py` |
+| 04 | Normalization and Rescaling | `NNSvignette_04_Normalization_and_Rescaling.Rmd` | `04_normalization_and_rescaling.py` |
+| 05 | Sampling and Simulation | `NNSvignette_05_Sampling.Rmd` | `05_sampling_and_simulation.py` |
+| 06 | Comparing Distributions | `NNSvignette_06_Comparing_Distributions.Rmd` | `06_comparing_distributions.py` |
+| 07 | Clustering and Regression | `NNSvignette_07_Clustering_and_Regression.Rmd` | `07_clustering_and_regression.py` |
+| 08 | Classification | `NNSvignette_08_Classification.Rmd` | `08_classification.py` |
+| 09 | Forecasting | `NNSvignette_09_Forecasting.Rmd` | `09_forecasting.py` |
+
+The unnumbered scripts remain available as focused implementation examples and
+backward-compatible entry points. New canonical material should be added to R
+first and then ported to the corresponding numbered Python example.
+
+Run one example:
+
+```bash
+uv run python examples/vignettes/02_partial_moments.py
+```
+
+Run the full example suite:
+
+```bash
+uv run python examples/run_all_vignettes.py
+```
diff --git a/examples/vignettes/manifest.yml b/examples/vignettes/manifest.yml
new file mode 100644
index 00000000..d948f50f
--- /dev/null
+++ b/examples/vignettes/manifest.yml
@@ -0,0 +1,43 @@
+canonical_source:
+  repository: OVVO-Financial/NNS
+  branch: NNS-Beta-Version
+  directory: vignettes
+  policy: R defines the statistical narrative, datasets, and section order.
+
+examples:
+  - id: '01'
+    topic: Overview
+    r: NNSvignette_01_Overview.Rmd
+    python: 01_overview.py
+  - id: '02'
+    topic: Partial Moments
+    r: NNSvignette_02_Partial_Moments.Rmd
+    python: 02_partial_moments.py
+  - id: '03'
+    topic: Correlation and Dependence
+    r: NNSvignette_03_Correlation_and_Dependence.Rmd
+    python: 03_correlation_and_dependence.py
+  - id: '04'
+    topic: Normalization and Rescaling
+    r: NNSvignette_04_Normalization_and_Rescaling.Rmd
+    python: 04_normalization_and_rescaling.py
+  - id: '05'
+    topic: Sampling and Simulation
+    r: NNSvignette_05_Sampling.Rmd
+    python: 05_sampling_and_simulation.py
+  - id: '06'
+    topic: Comparing Distributions
+    r: NNSvignette_06_Comparing_Distributions.Rmd
+    python: 06_comparing_distributions.py
+  - id: '07'
+    topic: Clustering and Regression
+    r: NNSvignette_07_Clustering_and_Regression.Rmd
+    python: 07_clustering_and_regression.py
+  - id: '08'
+    topic: Classification
+    r: NNSvignette_08_Classification.Rmd
+    python: 08_classification.py
+  - id: '09'
+    topic: Forecasting
+    r: NNSvignette_09_Forecasting.Rmd
+    python: 09_forecasting.py
diff --git a/tests/docs/test_canonical_vignette_manifest.py b/tests/docs/test_canonical_vignette_manifest.py
new file mode 100644
index 00000000..3359bf0f
--- /dev/null
+++ b/tests/docs/test_canonical_vignette_manifest.py
@@ -0,0 +1,35 @@
+from __future__ import annotations
+
+import re
+from pathlib import Path
+
+ROOT = Path(__file__).resolve().parents[2]
+EXAMPLE_DIR = ROOT / "examples" / "vignettes"
+EXPECTED = {
+    "01": ("NNSvignette_01_Overview.Rmd", "01_overview.py"),
+    "02": ("NNSvignette_02_Partial_Moments.Rmd", "02_partial_moments.py"),
+    "03": ("NNSvignette_03_Correlation_and_Dependence.Rmd", "03_correlation_and_dependence.py"),
+    "04": ("NNSvignette_04_Normalization_and_Rescaling.Rmd", "04_normalization_and_rescaling.py"),
+    "05": ("NNSvignette_05_Sampling.Rmd", "05_sampling_and_simulation.py"),
+    "06": ("NNSvignette_06_Comparing_Distributions.Rmd", "06_comparing_distributions.py"),
+    "07": ("NNSvignette_07_Clustering_and_Regression.Rmd", "07_clustering_and_regression.py"),
+    "08": ("NNSvignette_08_Classification.Rmd", "08_classification.py"),
+    "09": ("NNSvignette_09_Forecasting.Rmd", "09_forecasting.py"),
+}
+
+
+def test_canonical_vignette_set_is_complete_and_ordered() -> None:
+    numbered = sorted(path.name for path in EXAMPLE_DIR.glob("[0-9][0-9]_*.py"))
+    assert numbered == [python for _, python in EXPECTED.values()]
+
+
+def test_each_python_entrypoint_names_its_r_source() -> None:
+    for r_source, python_name in EXPECTED.values():
+        text = (EXAMPLE_DIR / python_name).read_text(encoding="utf-8")
+        assert r_source in text
+
+
+def test_manifest_matches_the_canonical_pairs() -> None:
+    text = (EXAMPLE_DIR / "manifest.yml").read_text(encoding="utf-8")
+    pairs = re.findall(r"\n\s+r: (\S+)\n\s+python: (\S+)", text)
+    assert pairs == list(EXPECTED.values())