Add Python vignettes for NNS R 13.0 curriculum

.github/workflows/native-backend-ci.yml

@@ -40,6 +40,9 @@ jobs:
       - name: Run parity from committed R cache
         run: NNS_R_CACHE_ONLY=1 python -m pytest -q tests/parity
 
+      - name: Run vignette examples
+        run: python -m pytest -q tests/docs/test_vignette_examples.py
+
       - name: Run ruff
         run: ruff check .
 

docs/vignettes/00_overview.md

@@ -0,0 +1,60 @@
+# 00 — Overview
+
+NNS (Nonlinear Nonparametric Statistics) builds its entire toolkit on **partial
+moments** — the pieces of variance that lie above and below a target. Because
+partial moments make no assumption of symmetry, linearity, or a parametric
+distribution, the same primitives reconstruct classical statistics (variance,
+covariance, the CDF) *and* extend naturally to nonlinear dependence,
+regression, forecasting, and stochastic dominance.
+
+This curriculum follows the R NNS vignettes:
+
+1. Partial moments — the foundational LPM/UPM primitives.
+2. Descriptive and distributional tools.
+3. Dependence and nonlinear association.
+4. Normalization and rescaling.
+5. Hypothesis testing: ANOVA and stochastic superiority.
+6. Regression, boosting, stacking, and causality.
+7. Time series forecasting.
+8. Simulation, bootstrap, and risk-neutral sampling.
+9. Portfolios and stochastic dominance.
+
+## Quick import
+
+```python
+import numpy as np
+from nns import (
+    lpm, upm, lpm_ratio,
+    nns_moments, nns_dep, nns_copula, pm_matrix,
+)
+```
+
+## A one-screen tour
+
+```python
+rng = np.random.default_rng(42)
+
+# Variance is the sum of second-degree partial moments about the mean.
+y = rng.normal(size=3000)
+mu = float(np.mean(y))
+n = y.size
+pm_variance = (lpm(2, mu, y) + upm(2, mu, y)) * (n / (n - 1))
+assert np.isclose(pm_variance, np.var(y, ddof=1))
+
+# The empirical CDF is LPM.ratio with degree 0.
+assert np.isclose(lpm_ratio(0, 0.0, y), np.mean(y <= 0.0))
+
+# Pearson correlation misses y = x**2; partial-moment dependence does not.
+x = rng.uniform(-1, 1, size=2000)
+yq = x**2 + rng.normal(scale=0.05, size=2000)
+print("Pearson:", np.corrcoef(x, yq)[0, 1])      # ~0
+print("Dependence:", nns_dep(x, yq)["Dependence"])  # clearly positive
+```
+
+Run the full tour:
+
+```bash
+python examples/vignettes/overview.py
+```
+
+The remaining vignettes unpack each of these ideas in turn.

docs/vignettes/01_partial_moments.md

@@ -0,0 +1,74 @@
+# 01 — Partial moments
+
+Partial moments split a distribution at a target `t`. The **lower partial
+moment** `lpm(degree, t, x)` accumulates deviations below `t`; the **upper
+partial moment** `upm(degree, t, x)` accumulates deviations above it. Their
+ratios give probabilities, and their inverses give quantiles. Everything else
+in NNS is built from these.
+
+```python
+import numpy as np
+from nns import lpm, upm, lpm_ratio, upm_ratio, lpm_var, upm_var
+
+rng = np.random.default_rng(123)
+x = rng.normal(size=100)
+mu = float(np.mean(x))
+n = x.size
+```
+
+## The mean as a partial-moment balance point
+
+The first-degree upper and lower partial moments about 0 balance at the mean:
+
+```python
+mean_via_pm = upm(1, 0.0, x) - lpm(1, 0.0, x)
+assert np.isclose(mean_via_pm, np.mean(x))
+```
+
+## Variance decomposition around the mean
+
+Second-degree partial moments about the mean sum to the **population**
+variance; multiply by `n / (n - 1)` for the sample variance:
+
+```python
+population_variance = upm(2, mu, x) + lpm(2, mu, x)
+sample_variance = population_variance * (n / (n - 1))
+assert np.isclose(sample_variance, np.var(x, ddof=1))
+```
+
+This is the central NNS idea: variance is not a monolithic quantity but the sum
+of an upside and a downside piece, each measurable on its own.
+
+## Empirical CDF via `lpm_ratio(0, t, x)`
+
+The degree-0 lower partial moment ratio is the proportion of mass at or below
+`t` — exactly the empirical CDF. `upm_ratio` is the complementary survival
+function:
+
+```python
+for t in (-1.0, 0.0, 1.0):
+    assert np.isclose(lpm_ratio(0, t, x), np.mean(x <= t))
+    assert np.isclose(upm_ratio(0, t, x), 1.0 - lpm_ratio(0, t, x))
+```
+
+## Value-at-risk quantiles
+
+`lpm_var(p, 0, x)` inverts the degree-0 CDF, returning the `p`-quantile;
+`upm_var(p, 0, x)` returns the right-tail `(1 - p)` quantile:
+
+```python
+p = np.array([0.05, 0.25, 0.5, 0.75, 0.95])
+left = np.array([lpm_var(q, 0.0, x) for q in p])
+assert np.allclose(left, np.quantile(x, p, method="linear"))
+
+right = np.array([upm_var(q, 0.0, x) for q in p])
+assert np.allclose(right, np.quantile(x, 1.0 - p, method="linear"))
+```
+
+For integer degrees 1–4, `lpm_var`/`upm_var` solve the exact partial-moment
+ratio inversion (the polynomial root-finder ported from R NNS 13.0 in PR #3),
+producing continuous VaR estimates rather than raw order statistics.
+
+```bash
+python examples/vignettes/partial_moments.py
+```

docs/vignettes/02_descriptive_distributional_tools.md

@@ -0,0 +1,64 @@
+# 02 — Descriptive and distributional tools
+
+Partial moments give a full descriptive toolkit: moment summaries, robust
+modes, covariance matrices, and quantile tables — all without distributional
+assumptions.
+
+```python
+import numpy as np
+from nns import nns_moments, nns_mode, pm_matrix, lpm_var, lpm_ratio
+
+rng = np.random.default_rng(123)
+x = rng.normal(size=200)
+y = rng.normal(size=200)
+```
+
+## Moment summaries
+
+`nns_moments` returns mean, variance, skewness, and kurtosis. The `population`
+flag toggles the `n/(n-1)` rescaling:
+
+```python
+nns_moments(x, population=True)   # {'mean', 'variance', 'skewness', 'kurtosis'}
+nns_moments(x, population=False)  # sample variance is larger
+```
+
+## Modes (continuous and discrete)
+
+`nns_mode` estimates a continuous mode by default, or returns discrete /
+multiple modes:
+
+```python
+nns_mode(x)                                                   # continuous estimate
+nns_mode(np.array([1, 2, 2, 3, 3, 4, 4, 5], dtype=float),
+         discrete=True, multi=True)                            # several modes
+```
+
+## Covariance reconstruction from a partial moment matrix
+
+`pm_matrix` returns the four co-partial-moment blocks (`clpm`, `cupm`, `dlpm`,
+`dupm`). The covariance matrix is recovered as `clpm + cupm - dlpm - dupm`:
+
+```python
+pm = pm_matrix(1, 1, "mean", np.column_stack((x, y)), True, names=["x", "y"])
+reconstructed = pm["clpm"] + pm["cupm"] - pm["dlpm"] - pm["dupm"]
+assert np.allclose(reconstructed, np.cov(x, y))
+```
+
+This mirrors the R vignette's covariance-matrix reassembly and shows that the
+classical covariance is just a difference of co-partial moments.
+
+## Quantile table via `lpm_var`
+
+A quantile table is a sweep of `lpm_var` over percentiles; `lpm_ratio` recovers
+the CDF at each threshold as a round-trip check:
+
+```python
+p = np.arange(0.05, 0.96, 0.1)
+thresholds = np.array([lpm_var(q, 0.0, x) for q in p])
+recovered = np.array([lpm_ratio(0, t, x) for t in thresholds])  # equals p
+```
+
+```bash
+python examples/vignettes/descriptive_distributional_tools.py
+```

docs/vignettes/03_dependence_nonlinear_association.md

@@ -0,0 +1,64 @@
+# 03 — Dependence and nonlinear association
+
+Pearson correlation measures *linear* co-movement. When a relationship is
+nonlinear, correlation can collapse toward zero even though the variables are
+perfectly dependent. NNS measures dependence directly from partial moments, so
+it sees structure the linear coefficient misses.
+
+```python
+import numpy as np
+from nns import nns_dep, nns_copula, pm_matrix
+```
+
+## Linear baseline
+
+For `y = 2x`, both correlation and dependence are ~1:
+
+```python
+x = np.arange(0.0, 3.01, 0.01)
+lin = nns_dep(x, 2.0 * x)
+# lin["Correlation"] ~ 1, lin["Dependence"] ~ 1
+```
+
+## Where Pearson collapses
+
+For `y = sin(x)` over many periods, Pearson correlation is weak while
+partial-moment dependence stays high:
+
+```python
+xs = np.arange(0.0, 12.0 * np.pi, np.pi / 100.0)
+ys = np.sin(xs)
+sine = nns_dep(xs, ys)
+# sine["Correlation"] ~ 0.20 (weak), sine["Dependence"] ~ 0.81 (strong)
+assert sine["Dependence"] > 3.0 * abs(sine["Correlation"])
+```
+
+The partial moment dependence vs Pearson correlation contrast is the whole
+point: correlation answers "how linear?", dependence answers "how related?".
+
+## Asymmetric dependence
+
+Dependence need not be symmetric — `D(y | x)` can differ from `D(x | y)`:
+
+```python
+asym_xy = nns_dep(xs, ys, asym=True)["Dependence"]
+asym_yx = nns_dep(ys, xs, asym=True)["Dependence"]
+```
+
+## Multivariate dependence and copulas
+
+`pm_matrix` exposes the co-partial-moment blocks for a frame, and `nns_copula`
+summarizes the joint dependence structure (near 0.5 for independent columns):
+
+```python
+rng = np.random.default_rng(123)
+frame = np.column_stack((rng.normal(size=1000),
+                         rng.normal(size=1000),
+                         rng.normal(size=1000)))
+pm = pm_matrix(1, 1, "mean", frame, True, names=["a", "b", "c"])
+nns_copula(frame, continuous=True)
+```
+
+```bash
+python examples/vignettes/dependence_nonlinear_association.py
+```

docs/vignettes/04_normalization_rescaling.md

@@ -0,0 +1,68 @@
+# 04 — Normalization and rescaling
+
+NNS provides two complementary transforms: `nns_norm` aligns variables onto a
+common scale (linearly or nonlinearly), and `nns_rescale` maps a vector onto an
+explicit interval or a risk-neutral target.
+
+```python
+import numpy as np
+from nns import nns_norm, nns_rescale
+```
+
+## `nns_norm` — linear and nonlinear
+
+Given columns with wildly different means and spreads, linear normalization
+brings them onto a shared mean scale:
+
+```python
+rng = np.random.default_rng(123)
+X = np.column_stack((
+    rng.normal(0, 1, 100),
+    rng.normal(0, 5, 100),
+    rng.normal(10, 1, 100),
+    rng.normal(10, 10, 100),
+))
+
+linear = nns_norm(X, linear=True)     # columns share a common mean
+nonlinear = nns_norm(X, linear=False) # partial-moment normalization
+```
+
+After linear normalization every column has the same mean — the precondition
+for comparing series on one axis.
+
+## `nns_rescale` — min-max
+
+Map a vector onto an explicit `[a, b]` interval:
+
+```python
+raw = np.array([-2.5, 0.2, 1.1, 3.7, 5.0])
+scaled = nns_rescale(raw, a=5.0, b=10.0, method="minmax")
+# scaled.min() == 5.0, scaled.max() == 10.0
+```
+
+## `nns_rescale` — risk-neutral
+
+The `"riskneutral"` method rescales a price path so its mean matches a
+risk-neutral target. With `type="Terminal"` the rescaled mean equals the
+forward `S0 * exp(r * T)`; with `type="Discounted"` it equals `S0`:
+
+```python
+s0, r, t = 100.0, 0.03, 1.0
+prices = s0 * np.exp(np.cumsum(rng.normal(0.0005, 0.02, 250)))
+
+terminal = nns_rescale(prices, a=s0, b=r, method="riskneutral",
+                       time_to_maturity=t, type="Terminal")
+assert np.isclose(terminal.mean(), s0 * np.exp(r * t))
+
+discounted = nns_rescale(prices, a=s0, b=r, method="riskneutral",
+                         time_to_maturity=t, type="Discounted")
+assert np.isclose(discounted.mean(), s0)
+```
+
+> Note: the R vignette also illustrates these transforms with overlaid
+> histograms. Plotting is optional and omitted from the docs tests; the
+> numeric invariants above are what matter for parity.
+
+```bash
+python examples/vignettes/normalization_rescaling.py
+```