From 2243d7b12048b6bf4ba484652b87a574bbb86a8f Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:32:49 -0400 Subject: [PATCH 01/17] Add shared support for instructional vignettes --- examples/_vignette_support.py | 180 ++++++++++++++++++++++++++++++++++ 1 file changed, 180 insertions(+) create mode 100644 examples/_vignette_support.py diff --git a/examples/_vignette_support.py b/examples/_vignette_support.py new file mode 100644 index 0000000..0f532ca --- /dev/null +++ b/examples/_vignette_support.py @@ -0,0 +1,180 @@ +"""Shared, non-vignette support for the canonical Python ports. + +The numbered files in ``examples/vignettes`` are intentionally readable as +Jupytext-style instructional scripts. This module keeps repetitive concerns +(dataset loading, output folders, compact result rendering, and figure saving) +out of the statistical narrative. +""" +from __future__ import annotations + +import csv +import os +from pathlib import Path +from pprint import pformat +from typing import Any, Iterable, Sequence + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +EXAMPLES_DIR = Path(__file__).resolve().parent +DATA_DIR = EXAMPLES_DIR / "data" +OUTPUT_ROOT = EXAMPLES_DIR / "vignettes" / "output" + + +def fast_mode() -> bool: + """Return whether CI requested reduced repetitions for expensive sections.""" + return os.environ.get("NNS_VIGNETTE_FAST", "0") == "1" + + +def output_dir(script_file: str) -> Path: + path = OUTPUT_ROOT / Path(script_file).stem + path.mkdir(parents=True, exist_ok=True) + return path + + +def section(title: str) -> None: + line = "=" * len(title) + print(f"\n{title}\n{line}") + + +def subsection(title: str) -> None: + print(f"\n--- {title} ---") + + +def gap(message: str) -> None: + """Make a genuine R/Python API difference impossible to overlook.""" + print(f"\nPYTHON API GAP: {message}") + + +def note(message: str) -> None: + print(f"\nNOTE: {message}") + + +def show(label: str, value: Any, *, precision: int = 6) -> None: + """Display important numerical structures without dumping huge arrays.""" + print(f"\n{label}:") + if isinstance(value, np.ndarray): + print(np.array2string(value, precision=precision, threshold=40, edgeitems=5)) + return + if isinstance(value, dict): + compact = {key: _compact(item) for key, item in value.items()} + print(pformat(compact, width=100, sort_dicts=False)) + return + print(value) + + +def _compact(value: Any) -> Any: + if isinstance(value, np.ndarray): + if value.size <= 20: + return np.round(value, 6).tolist() + return { + "shape": list(value.shape), + "head": np.round(value.reshape(-1)[:5], 6).tolist(), + "tail": np.round(value.reshape(-1)[-5:], 6).tolist(), + } + if isinstance(value, dict): + return {key: _compact(item) for key, item in value.items()} + if isinstance(value, (np.floating, np.integer)): + return value.item() + return value + + +def table(headers: Sequence[str], rows: Iterable[Sequence[Any]], *, precision: int = 5) -> None: + rendered = [] + for row in rows: + rendered.append( + [f"{item:.{precision}f}" if isinstance(item, (float, np.floating)) else str(item) for item in row] + ) + widths = [len(str(header)) for header in headers] + for row in rendered: + widths = [max(width, len(cell)) for width, cell in zip(widths, row, strict=True)] + print(" ".join(str(header).ljust(width) for header, width in zip(headers, widths, strict=True))) + print(" ".join("-" * width for width in widths)) + for row in rendered: + print(" ".join(cell.ljust(width) for cell, width in zip(row, widths, strict=True))) + + +def save_figure(fig: plt.Figure, directory: Path, filename: str) -> Path: + path = directory / filename + fig.tight_layout() + fig.savefig(path, dpi=160, bbox_inches="tight") + plt.close(fig) + print(f"Saved figure: {path}") + return path + + +def load_iris() -> tuple[np.ndarray, np.ndarray, list[str]]: + rows = _read_csv("iris.csv") + feature_names = ["Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"] + x = np.asarray([[float(row[name]) for name in feature_names] for row in rows], dtype=float) + levels = ["setosa", "versicolor", "virginica"] + y = np.asarray([levels.index(row["Species"]) + 1 for row in rows], dtype=float) + return x, y, levels + + +def load_mtcars() -> dict[str, np.ndarray]: + rows = _read_csv("mtcars.csv") + numeric = ["mpg", "cyl", "disp", "hp", "drat", "wt", "qsec", "vs", "am", "gear", "carb"] + return {name: np.asarray([float(row[name]) for row in rows], dtype=float) for name in numeric} + + +def load_air_passengers() -> np.ndarray: + rows = _read_csv("AirPassengers.csv") + return np.asarray([float(row["value"]) for row in rows], dtype=float) + + +def _read_csv(filename: str) -> list[dict[str, str]]: + with (DATA_DIR / filename).open(newline="", encoding="utf-8") as handle: + return list(csv.DictReader(handle)) + + +def empirical_cdf(values: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + ordered = np.sort(np.asarray(values, dtype=float)) + probs = np.arange(1, ordered.size + 1, dtype=float) / ordered.size + return ordered, probs + + +def partition_scatter(ax: plt.Axes, result: dict[str, Any], title: str) -> None: + dt = result["dt"] + labels = np.asarray(dt["quadrant"], dtype=str) + unique, codes = np.unique(labels, return_inverse=True) + ax.scatter(dt["x"], dt["y"], c=codes, cmap="tab20", s=14, alpha=0.7) + rp = result["regression.points"] + ax.scatter(rp["x"], rp["y"], marker="x", s=55, linewidths=1.6) + ax.set_title(f"{title}\n{len(unique)} terminal paths; order={result['order']}") + ax.set_xlabel("x") + ax.set_ylabel("y") + + +def regression_scatter( + ax: plt.Axes, + x: np.ndarray, + y: np.ndarray, + result: dict[str, Any], + title: str, +) -> None: + fitted = result.get("Fitted.xy") + ax.scatter(x, y, s=12, alpha=0.45, label="observed") + if isinstance(fitted, dict): + order = np.argsort(np.asarray(fitted["x"], dtype=float)) + ax.plot(np.asarray(fitted["x"])[order], np.asarray(fitted["y.hat"])[order], label="NNS fit") + if "conf.int.neg" in fitted and "conf.int.pos" in fitted: + ax.fill_between( + np.asarray(fitted["x"])[order], + np.asarray(fitted["conf.int.neg"])[order], + np.asarray(fitted["conf.int.pos"])[order], + alpha=0.2, + label="confidence interval", + ) + rp = result.get("regression.points") + if isinstance(rp, dict): + ax.scatter(rp["x"], rp["y"], marker="x", s=45, label="regression points") + ax.set_title(title) + ax.legend() + + +def figure_grid(nrows: int, ncols: int, *, width: float = 5.0, height: float = 4.0): + return plt.subplots(nrows, ncols, figsize=(width * ncols, height * nrows), squeeze=False) From 670671b6955651eb339beaffc5a806c03a0cd41b Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:33:12 -0400 Subject: [PATCH 02/17] Embed canonical R iris dataset --- examples/data/iris.csv | 151 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 151 insertions(+) create mode 100644 examples/data/iris.csv diff --git a/examples/data/iris.csv b/examples/data/iris.csv new file mode 100644 index 0000000..6bdc84d --- /dev/null +++ b/examples/data/iris.csv @@ -0,0 +1,151 @@ +rownames,Sepal.Length,Sepal.Width,Petal.Length,Petal.Width,Species +1,5.1,3.5,1.4,0.2,setosa +2,4.9,3,1.4,0.2,setosa +3,4.7,3.2,1.3,0.2,setosa +4,4.6,3.1,1.5,0.2,setosa +5,5,3.6,1.4,0.2,setosa +6,5.4,3.9,1.7,0.4,setosa +7,4.6,3.4,1.4,0.3,setosa +8,5,3.4,1.5,0.2,setosa +9,4.4,2.9,1.4,0.2,setosa +10,4.9,3.1,1.5,0.1,setosa +11,5.4,3.7,1.5,0.2,setosa +12,4.8,3.4,1.6,0.2,setosa +13,4.8,3,1.4,0.1,setosa +14,4.3,3,1.1,0.1,setosa +15,5.8,4,1.2,0.2,setosa +16,5.7,4.4,1.5,0.4,setosa +17,5.4,3.9,1.3,0.4,setosa +18,5.1,3.5,1.4,0.3,setosa +19,5.7,3.8,1.7,0.3,setosa +20,5.1,3.8,1.5,0.3,setosa +21,5.4,3.4,1.7,0.2,setosa +22,5.1,3.7,1.5,0.4,setosa +23,4.6,3.6,1,0.2,setosa +24,5.1,3.3,1.7,0.5,setosa +25,4.8,3.4,1.9,0.2,setosa +26,5,3,1.6,0.2,setosa +27,5,3.4,1.6,0.4,setosa +28,5.2,3.5,1.5,0.2,setosa +29,5.2,3.4,1.4,0.2,setosa +30,4.7,3.2,1.6,0.2,setosa +31,4.8,3.1,1.6,0.2,setosa +32,5.4,3.4,1.5,0.4,setosa +33,5.2,4.1,1.5,0.1,setosa +34,5.5,4.2,1.4,0.2,setosa +35,4.9,3.1,1.5,0.2,setosa +36,5,3.2,1.2,0.2,setosa +37,5.5,3.5,1.3,0.2,setosa +38,4.9,3.6,1.4,0.1,setosa +39,4.4,3,1.3,0.2,setosa +40,5.1,3.4,1.5,0.2,setosa +41,5,3.5,1.3,0.3,setosa +42,4.5,2.3,1.3,0.3,setosa +43,4.4,3.2,1.3,0.2,setosa +44,5,3.5,1.6,0.6,setosa +45,5.1,3.8,1.9,0.4,setosa +46,4.8,3,1.4,0.3,setosa +47,5.1,3.8,1.6,0.2,setosa +48,4.6,3.2,1.4,0.2,setosa +49,5.3,3.7,1.5,0.2,setosa +50,5,3.3,1.4,0.2,setosa +51,7,3.2,4.7,1.4,versicolor +52,6.4,3.2,4.5,1.5,versicolor +53,6.9,3.1,4.9,1.5,versicolor +54,5.5,2.3,4,1.3,versicolor +55,6.5,2.8,4.6,1.5,versicolor +56,5.7,2.8,4.5,1.3,versicolor +57,6.3,3.3,4.7,1.6,versicolor +58,4.9,2.4,3.3,1,versicolor +59,6.6,2.9,4.6,1.3,versicolor +60,5.2,2.7,3.9,1.4,versicolor +61,5,2,3.5,1,versicolor +62,5.9,3,4.2,1.5,versicolor +63,6,2.2,4,1,versicolor +64,6.1,2.9,4.7,1.4,versicolor +65,5.6,2.9,3.6,1.3,versicolor +66,6.7,3.1,4.4,1.4,versicolor +67,5.6,3,4.5,1.5,versicolor +68,5.8,2.7,4.1,1,versicolor +69,6.2,2.2,4.5,1.5,versicolor +70,5.6,2.5,3.9,1.1,versicolor +71,5.9,3.2,4.8,1.8,versicolor +72,6.1,2.8,4,1.3,versicolor +73,6.3,2.5,4.9,1.5,versicolor +74,6.1,2.8,4.7,1.2,versicolor +75,6.4,2.9,4.3,1.3,versicolor +76,6.6,3,4.4,1.4,versicolor +77,6.8,2.8,4.8,1.4,versicolor +78,6.7,3,5,1.7,versicolor +79,6,2.9,4.5,1.5,versicolor +80,5.7,2.6,3.5,1,versicolor +81,5.5,2.4,3.8,1.1,versicolor +82,5.5,2.4,3.7,1,versicolor +83,5.8,2.7,3.9,1.2,versicolor +84,6,2.7,5.1,1.6,versicolor +85,5.4,3,4.5,1.5,versicolor +86,6,3.4,4.5,1.6,versicolor +87,6.7,3.1,4.7,1.5,versicolor +88,6.3,2.3,4.4,1.3,versicolor +89,5.6,3,4.1,1.3,versicolor +90,5.5,2.5,4,1.3,versicolor +91,5.5,2.6,4.4,1.2,versicolor +92,6.1,3,4.6,1.4,versicolor +93,5.8,2.6,4,1.2,versicolor +94,5,2.3,3.3,1,versicolor +95,5.6,2.7,4.2,1.3,versicolor +96,5.7,3,4.2,1.2,versicolor +97,5.7,2.9,4.2,1.3,versicolor +98,6.2,2.9,4.3,1.3,versicolor +99,5.1,2.5,3,1.1,versicolor +100,5.7,2.8,4.1,1.3,versicolor +101,6.3,3.3,6,2.5,virginica +102,5.8,2.7,5.1,1.9,virginica +103,7.1,3,5.9,2.1,virginica +104,6.3,2.9,5.6,1.8,virginica +105,6.5,3,5.8,2.2,virginica +106,7.6,3,6.6,2.1,virginica +107,4.9,2.5,4.5,1.7,virginica +108,7.3,2.9,6.3,1.8,virginica +109,6.7,2.5,5.8,1.8,virginica +110,7.2,3.6,6.1,2.5,virginica +111,6.5,3.2,5.1,2,virginica +112,6.4,2.7,5.3,1.9,virginica +113,6.8,3,5.5,2.1,virginica +114,5.7,2.5,5,2,virginica +115,5.8,2.8,5.1,2.4,virginica +116,6.4,3.2,5.3,2.3,virginica +117,6.5,3,5.5,1.8,virginica +118,7.7,3.8,6.7,2.2,virginica +119,7.7,2.6,6.9,2.3,virginica +120,6,2.2,5,1.5,virginica +121,6.9,3.2,5.7,2.3,virginica +122,5.6,2.8,4.9,2,virginica +123,7.7,2.8,6.7,2,virginica +124,6.3,2.7,4.9,1.8,virginica +125,6.7,3.3,5.7,2.1,virginica +126,7.2,3.2,6,1.8,virginica +127,6.2,2.8,4.8,1.8,virginica +128,6.1,3,4.9,1.8,virginica +129,6.4,2.8,5.6,2.1,virginica +130,7.2,3,5.8,1.6,virginica +131,7.4,2.8,6.1,1.9,virginica +132,7.9,3.8,6.4,2,virginica +133,6.4,2.8,5.6,2.2,virginica +134,6.3,2.8,5.1,1.5,virginica +135,6.1,2.6,5.6,1.4,virginica +136,7.7,3,6.1,2.3,virginica +137,6.3,3.4,5.6,2.4,virginica +138,6.4,3.1,5.5,1.8,virginica +139,6,3,4.8,1.8,virginica +140,6.9,3.1,5.4,2.1,virginica +141,6.7,3.1,5.6,2.4,virginica +142,6.9,3.1,5.1,2.3,virginica +143,5.8,2.7,5.1,1.9,virginica +144,6.8,3.2,5.9,2.3,virginica +145,6.7,3.3,5.7,2.5,virginica +146,6.7,3,5.2,2.3,virginica +147,6.3,2.5,5,1.9,virginica +148,6.5,3,5.2,2,virginica +149,6.2,3.4,5.4,2.3,virginica +150,5.9,3,5.1,1.8,virginica From ad58fcbacd56df85520929e65cc5148c506d77c8 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:33:25 -0400 Subject: [PATCH 03/17] Embed canonical R mtcars dataset --- examples/data/mtcars.csv | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 examples/data/mtcars.csv diff --git a/examples/data/mtcars.csv b/examples/data/mtcars.csv new file mode 100644 index 0000000..6bb1f90 --- /dev/null +++ b/examples/data/mtcars.csv @@ -0,0 +1,33 @@ +rownames,mpg,cyl,disp,hp,drat,wt,qsec,vs,am,gear,carb +Mazda RX4,21,6,160,110,3.9,2.62,16.46,0,1,4,4 +Mazda RX4 Wag,21,6,160,110,3.9,2.875,17.02,0,1,4,4 +Datsun 710,22.8,4,108,93,3.85,2.32,18.61,1,1,4,1 +Hornet 4 Drive,21.4,6,258,110,3.08,3.215,19.44,1,0,3,1 +Hornet Sportabout,18.7,8,360,175,3.15,3.44,17.02,0,0,3,2 +Valiant,18.1,6,225,105,2.76,3.46,20.22,1,0,3,1 +Duster 360,14.3,8,360,245,3.21,3.57,15.84,0,0,3,4 +Merc 240D,24.4,4,146.7,62,3.69,3.19,20,1,0,4,2 +Merc 230,22.8,4,140.8,95,3.92,3.15,22.9,1,0,4,2 +Merc 280,19.2,6,167.6,123,3.92,3.44,18.3,1,0,4,4 +Merc 280C,17.8,6,167.6,123,3.92,3.44,18.9,1,0,4,4 +Merc 450SE,16.4,8,275.8,180,3.07,4.07,17.4,0,0,3,3 +Merc 450SL,17.3,8,275.8,180,3.07,3.73,17.6,0,0,3,3 +Merc 450SLC,15.2,8,275.8,180,3.07,3.78,18,0,0,3,3 +Cadillac Fleetwood,10.4,8,472,205,2.93,5.25,17.98,0,0,3,4 +Lincoln Continental,10.4,8,460,215,3,5.424,17.82,0,0,3,4 +Chrysler Imperial,14.7,8,440,230,3.23,5.345,17.42,0,0,3,4 +Fiat 128,32.4,4,78.7,66,4.08,2.2,19.47,1,1,4,1 +Honda Civic,30.4,4,75.7,52,4.93,1.615,18.52,1,1,4,2 +Toyota Corolla,33.9,4,71.1,65,4.22,1.835,19.9,1,1,4,1 +Toyota Corona,21.5,4,120.1,97,3.7,2.465,20.01,1,0,3,1 +Dodge Challenger,15.5,8,318,150,2.76,3.52,16.87,0,0,3,2 +AMC Javelin,15.2,8,304,150,3.15,3.435,17.3,0,0,3,2 +Camaro Z28,13.3,8,350,245,3.73,3.84,15.41,0,0,3,4 +Pontiac Firebird,19.2,8,400,175,3.08,3.845,17.05,0,0,3,2 +Fiat X1-9,27.3,4,79,66,4.08,1.935,18.9,1,1,4,1 +Porsche 914-2,26,4,120.3,91,4.43,2.14,16.7,0,1,5,2 +Lotus Europa,30.4,4,95.1,113,3.77,1.513,16.9,1,1,5,2 +Ford Pantera L,15.8,8,351,264,4.22,3.17,14.5,0,1,5,4 +Ferrari Dino,19.7,6,145,175,3.62,2.77,15.5,0,1,5,6 +Maserati Bora,15,8,301,335,3.54,3.57,14.6,0,1,5,8 +Volvo 142E,21.4,4,121,109,4.11,2.78,18.6,1,1,4,2 From bc419a960f0c210a6b8795a00f1e4089499a12bc Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:33:40 -0400 Subject: [PATCH 04/17] Embed canonical R AirPassengers dataset --- examples/data/AirPassengers.csv | 145 ++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) create mode 100644 examples/data/AirPassengers.csv diff --git a/examples/data/AirPassengers.csv b/examples/data/AirPassengers.csv new file mode 100644 index 0000000..ebb4c40 --- /dev/null +++ b/examples/data/AirPassengers.csv @@ -0,0 +1,145 @@ +rownames,time,value +1,1949,112 +2,1949.08333333333,118 +3,1949.16666666667,132 +4,1949.25,129 +5,1949.33333333333,121 +6,1949.41666666667,135 +7,1949.5,148 +8,1949.58333333333,148 +9,1949.66666666667,136 +10,1949.75,119 +11,1949.83333333333,104 +12,1949.91666666667,118 +13,1950,115 +14,1950.08333333333,126 +15,1950.16666666667,141 +16,1950.25,135 +17,1950.33333333333,125 +18,1950.41666666667,149 +19,1950.5,170 +20,1950.58333333333,170 +21,1950.66666666667,158 +22,1950.75,133 +23,1950.83333333333,114 +24,1950.91666666667,140 +25,1951,145 +26,1951.08333333333,150 +27,1951.16666666667,178 +28,1951.25,163 +29,1951.33333333333,172 +30,1951.41666666667,178 +31,1951.5,199 +32,1951.58333333333,199 +33,1951.66666666667,184 +34,1951.75,162 +35,1951.83333333333,146 +36,1951.91666666667,166 +37,1952,171 +38,1952.08333333333,180 +39,1952.16666666667,193 +40,1952.25,181 +41,1952.33333333333,183 +42,1952.41666666667,218 +43,1952.5,230 +44,1952.58333333333,242 +45,1952.66666666667,209 +46,1952.75,191 +47,1952.83333333333,172 +48,1952.91666666667,194 +49,1953,196 +50,1953.08333333333,196 +51,1953.16666666667,236 +52,1953.25,235 +53,1953.33333333333,229 +54,1953.41666666667,243 +55,1953.5,264 +56,1953.58333333333,272 +57,1953.66666666667,237 +58,1953.75,211 +59,1953.83333333333,180 +60,1953.91666666667,201 +61,1954,204 +62,1954.08333333333,188 +63,1954.16666666667,235 +64,1954.25,227 +65,1954.33333333333,234 +66,1954.41666666667,264 +67,1954.5,302 +68,1954.58333333333,293 +69,1954.66666666667,259 +70,1954.75,229 +71,1954.83333333333,203 +72,1954.91666666667,229 +73,1955,242 +74,1955.08333333334,233 +75,1955.16666666667,267 +76,1955.25,269 +77,1955.33333333334,270 +78,1955.41666666667,315 +79,1955.5,364 +80,1955.58333333334,347 +81,1955.66666666667,312 +82,1955.75,274 +83,1955.83333333334,237 +84,1955.91666666667,278 +85,1956,284 +86,1956.08333333334,277 +87,1956.16666666667,317 +88,1956.25,313 +89,1956.33333333334,318 +90,1956.41666666667,374 +91,1956.5,413 +92,1956.58333333334,405 +93,1956.66666666667,355 +94,1956.75,306 +95,1956.83333333334,271 +96,1956.91666666667,306 +97,1957,315 +98,1957.08333333334,301 +99,1957.16666666667,356 +100,1957.25,348 +101,1957.33333333334,355 +102,1957.41666666667,422 +103,1957.5,465 +104,1957.58333333334,467 +105,1957.66666666667,404 +106,1957.75,347 +107,1957.83333333334,305 +108,1957.91666666667,336 +109,1958,340 +110,1958.08333333334,318 +111,1958.16666666667,362 +112,1958.25,348 +113,1958.33333333334,363 +114,1958.41666666667,435 +115,1958.5,491 +116,1958.58333333334,505 +117,1958.66666666667,404 +118,1958.75,359 +119,1958.83333333334,310 +120,1958.91666666667,337 +121,1959,360 +122,1959.08333333334,342 +123,1959.16666666667,406 +124,1959.25,396 +125,1959.33333333334,420 +126,1959.41666666667,472 +127,1959.5,548 +128,1959.58333333334,559 +129,1959.66666666667,463 +130,1959.75,407 +131,1959.83333333334,362 +132,1959.91666666667,405 +133,1960,417 +134,1960.08333333334,391 +135,1960.16666666667,419 +136,1960.25,461 +137,1960.33333333334,472 +138,1960.41666666667,535 +139,1960.5,622 +140,1960.58333333334,606 +141,1960.66666666667,508 +142,1960.75,461 +143,1960.83333333334,390 +144,1960.91666666667,432 From 95a86cd5cb2699bd1c356378027a4a64f01aa250 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:40:16 -0400 Subject: [PATCH 05/17] Port the partial moments vignette section by section --- examples/vignettes/02_partial_moments.py | 203 +++++++++++++++++------ 1 file changed, 154 insertions(+), 49 deletions(-) diff --git a/examples/vignettes/02_partial_moments.py b/examples/vignettes/02_partial_moments.py index 18a0af1..c8aec00 100644 --- a/examples/vignettes/02_partial_moments.py +++ b/examples/vignettes/02_partial_moments.py @@ -1,16 +1,25 @@ -"""Canonical vignette 02: Partial moments. +# %% [markdown] +# # 02. Getting Started with NNS: Partial Moments +# +# This is a section-for-section Python port of +# `NNSvignette_02_Partial_Moments.Rmd`. The R vignette remains the source of +# truth for the statistical narrative. NumPy uses a different seeded normal +# stream than R, so synthetic draws are distributionally equivalent rather +# than numerically identical. -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 sys +from pathlib import Path + +import matplotlib.pyplot as plt import numpy as np +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) + +from examples._vignette_support import empirical_cdf, output_dir, save_figure, section, show, subsection from nns import ( co_lpm, co_upm, @@ -18,62 +27,158 @@ d_upm, lpm, lpm_ratio, - lpm_var, + nns_cdf, nns_mode, nns_moments, pm_matrix, upm, ) +OUT = output_dir(__file__) + def main() -> None: + section("Partial Moments") + print( + "Parsing variance into lower, upper, co- and divergent partial moments " + "reveals distributional information unavailable from a single variance statistic." + ) + rng = np.random.default_rng(123) x = rng.normal(size=100) y = rng.normal(size=100) n = x.size + # %% [markdown] + # ## Mean + subsection("Mean") + mean_standard = float(np.mean(x)) 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) - + show("NumPy mean(x)", mean_standard) + show("UPM(1, 0, x) - LPM(1, 0, x)", mean_pm) + + # %% [markdown] + # ## Variance + subsection("Variance") + sample_variance = float(np.var(x, ddof=1)) + population_variance = float(upm(2, mean_standard, x) + lpm(2, mean_standard, x)) + sample_variance_pm = population_variance * n / (n - 1) + covariance_of_self = float( + co_lpm(1, x, x, mean_standard, mean_standard) + + co_upm(1, x, x, mean_standard, mean_standard) + - d_lpm(1, 1, x, x, mean_standard, mean_standard) + - d_upm(1, 1, x, x, mean_standard, mean_standard) + ) + show("Sample variance (NumPy)", sample_variance) + show("Sample variance from partial moments", sample_variance_pm) + show("Population adjustment of sample variance", sample_variance * (n - 1) / n) + show("Population variance from partial moments", population_variance) + show("Variance as covariance with itself", covariance_of_self) + + # %% [markdown] + # ## Standard Deviation + subsection("Standard Deviation") + show("Sample standard deviation", float(np.std(x, ddof=1))) + show("Standard deviation from partial moments", sample_variance_pm**0.5) + + # %% [markdown] + # ## First 4 Moments + subsection("First 4 Moments") + show("NNS moments (population)", nns_moments(x)) + show("NNS moments (sample)", nns_moments(x, population=False)) + + # %% [markdown] + # ## Statistical Mode of a Continuous Distribution + subsection("Statistical Mode of a Continuous Distribution") + show("Continuous NNS mode", nns_mode(x)) + show( + "Discrete multiple modes", + nns_mode(np.array([1, 2, 2, 3, 3, 4, 4, 5], dtype=float), discrete=True, multi=True), + ) + + # %% [markdown] + # ## Covariance + subsection("Covariance") + covariance_standard = float(np.cov(x, y, ddof=1)[0, 1]) + covariance_pm = float( + ( + co_lpm(1, x, y, np.mean(x), np.mean(y)) + + co_upm(1, x, y, np.mean(x), np.mean(y)) + - d_lpm(1, 1, x, y, np.mean(x), np.mean(y)) + - d_upm(1, 1, x, y, np.mean(x), np.mean(y)) + ) + * n + / (n - 1) + ) + show("Sample covariance (NumPy)", covariance_standard) + show("Sample covariance from co-partial moments", covariance_pm) + + # %% [markdown] + # ## Covariance Elements and Covariance Matrix + subsection("Covariance Elements and Covariance Matrix") + print("Sigma = CLPM + CUPM - DLPM - DUPM") + cov_mtx = pm_matrix(1, 1, "mean", np.column_stack((x, y)), pop_adj=True, names=["x", "y"]) + show("PM.matrix equivalent", cov_mtx) + reconstructed = cov_mtx["clpm"] + cov_mtx["cupm"] - cov_mtx["dlpm"] - cov_mtx["dupm"] + show("Reassembled covariance matrix", reconstructed) + show("Standard covariance matrix", np.cov(np.column_stack((x, y)), rowvar=False, ddof=1)) + + # %% [markdown] + # ## Pearson Correlation + subsection("Pearson Correlation") + sd_x = sample_variance_pm**0.5 + y_pop_var = float(upm(2, np.mean(y), y) + lpm(2, np.mean(y), y)) + sd_y = (y_pop_var * n / (n - 1)) ** 0.5 + show("Pearson correlation (NumPy)", float(np.corrcoef(x, y)[0, 1])) + show("Pearson correlation reconstructed from partial moments", covariance_pm / (sd_x * sd_y)) + + # %% [markdown] + # ## CDFs (Discrete and Continuous) + subsection("CDFs (Discrete and Continuous)") 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)) + show("Empirical CDF at 0 and 1", np.array([np.mean(x <= target) for target in targets])) + show("LPM degree-0 CDF at 0 and 1", lpm(0, targets, x)) + show("Joint CDF at (0, 0)", co_lpm(0, x, y, 0.0, 0.0)) + show("Vectorized joint CDF at (0,0) and (1,1)", co_lpm(0, x, y, targets, targets)) + + ux = np.asarray(lpm_ratio(0, x, x), dtype=float) + uy = np.asarray(lpm_ratio(0, y, y), dtype=float) + show("Copula at (0.5, 0.5)", co_lpm(0, ux, uy, 0.5, 0.5)) + show("Continuous NNS.CDF at 1", nns_cdf(x, degree=1, target=1.0)) + show("Continuous NNS.CDF at 1 with mean target", nns_cdf(x, degree=1, target=float(np.mean(x)))) + show("Survival function at 1", nns_cdf(x, degree=1, target=1.0, type="survival")) + + ordered, ecdf_values = empirical_cdf(x) + fig, ax = plt.subplots(figsize=(7, 4.5)) + ax.step(ordered, ecdf_values, where="post", label="empirical CDF") + ax.scatter(ordered, lpm(0, ordered, x), s=13, label="LPM degree-0 CDF") + ax.set_title("Empirical CDF and the degree-0 lower partial moment CDF") + ax.set_xlabel("x") + ax.set_ylabel("F(x)") + ax.legend() + save_figure(fig, OUT, "cdf_equivalence.png") + + # %% [markdown] + # ## Numerical Integration + subsection("Numerical Integration") + grid = np.arange(0.0, 1.0001, 0.001) + values = grid**2 + integral = float((upm(1, 0.0, values) - lpm(1, 0.0, values)) * (1.0 - 0.0)) + total_area = float((upm(1, 0.0, values) + lpm(1, 0.0, values)) * (1.0 - 0.0)) + show("Partial-moment approximation to integral_0^1 x^2 dx", integral) + show("Total unsigned area", total_area) + + # %% [markdown] + # ## Bayes' Theorem + subsection("Bayes' Theorem") + print("P(A > 0 | B > 0) = Co.UPM(0, A, B, 0, 0) / UPM(0, 0, B)") + conditional_probability = float(co_upm(0, x, y, 0.0, 0.0) / upm(0, 0.0, y)) + show("P(x > 0 | y > 0)", conditional_probability) + + np.testing.assert_allclose(mean_pm, mean_standard) + np.testing.assert_allclose(sample_variance_pm, sample_variance) + np.testing.assert_allclose(covariance_pm, covariance_standard) + np.testing.assert_allclose(reconstructed, np.cov(np.column_stack((x, y)), rowvar=False, ddof=1)) if __name__ == "__main__": From aba344ae7e812bfe6ad880ceaa2f68135d39ed6b Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:40:44 -0400 Subject: [PATCH 06/17] Port the dependence vignette with figures and explicit gaps --- .../03_correlation_and_dependence.py | 169 ++++++++++++++---- 1 file changed, 136 insertions(+), 33 deletions(-) diff --git a/examples/vignettes/03_correlation_and_dependence.py b/examples/vignettes/03_correlation_and_dependence.py index 499d8c9..cb276bc 100644 --- a/examples/vignettes/03_correlation_and_dependence.py +++ b/examples/vignettes/03_correlation_and_dependence.py @@ -1,47 +1,150 @@ -"""Canonical vignette 03: Correlation and dependence. +# %% [markdown] +# # 03. Getting Started with NNS: Correlation and Dependence +# +# Section-for-section port of the R vignette. Figures are saved beside the +# script under `output/03_correlation_and_dependence/`. -Source of truth: - NNS/vignettes/NNSvignette_03_Correlation_and_Dependence.Rmd -""" from __future__ import annotations +import sys +from pathlib import Path + +import matplotlib.pyplot as plt import numpy as np -from nns import nns_copula, nns_dep +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) + +from examples._vignette_support import gap, output_dir, partition_scatter, save_figure, section, show, subsection +from nns import nns_copula, nns_dep, nns_part + +OUT = output_dir(__file__) + + +def relationship(name: str, x: np.ndarray, y: np.ndarray, *, obs_req: int = 8) -> None: + subsection(name) + show("Pearson correlation", float(np.corrcoef(x, y)[0, 1])) + show("NNS correlation and dependence", nns_dep(x, y)) + part = nns_part(x, y, order=3, obs_req=obs_req) + fig, axes = plt.subplots(1, 2, figsize=(11, 4.3)) + axes[0].scatter(x, y, s=12, alpha=0.65) + axes[0].set_title(name) + axes[0].set_xlabel("x") + axes[0].set_ylabel("y") + partition_scatter(axes[1], part, "NNS partial-moment partition") + save_figure(fig, OUT, name.lower().replace(" ", "_") + ".png") 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) + section("Correlation and Dependence") + print( + "Pearson correlation measures linear association. NNS separates signed " + "correlation from unsigned signal-to-noise dependence so nonlinear relationships " + "need not disappear merely because their global slope is small." + ) + gap( + "R's NNS.part(..., Voronoi=TRUE) draws Voronoi regions internally. The Python " + "partition API returns the same quadrant paths and regression points but has no " + "Voronoi flag, so the figures color the returned partition memberships directly." + ) + + # %% [markdown] + # ## Linear Equivalence + x = np.arange(0.0, 3.001, 0.01) + relationship("Linear Equivalence", x, 2.0 * x) - 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"] + # %% [markdown] + # ## Nonlinear Relationship + relationship("Nonlinear Relationship", x, x**10) + + # %% [markdown] + # ## Cyclic Relationship + x_cycle = np.arange(0.0, 12.0 * np.pi + np.pi / 200.0, np.pi / 100.0) + y_cycle = np.sin(x_cycle) + relationship("Cyclic Relationship", x_cycle, y_cycle, obs_req=0) + + # %% [markdown] + # ## Asymmetrical Analysis + subsection("Asymmetrical Analysis") + show("Pearson cor(x, y)", float(np.corrcoef(x_cycle, y_cycle)[0, 1])) + show("NNS.dep(x, y, asym=True)", nns_dep(x_cycle, y_cycle, asym=True)) + show("Pearson cor(y, x)", float(np.corrcoef(y_cycle, x_cycle)[0, 1])) + show("NNS.dep(y, x, asym=True)", nns_dep(y_cycle, x_cycle, asym=True)) + + # %% [markdown] + # ## Dependence + subsection("Dependence") + rng = np.random.default_rng(123) + cloud = rng.uniform(-1.0, 1.0, size=(10000, 2)) + radius2 = np.sum(cloud**2, axis=1) + ring = cloud[(radius2 <= 1.0) & (radius2 >= 0.95)] + show("NNS dependence for annulus", nns_dep(ring[:, 0], ring[:, 1])) + ring_part = nns_part(ring[:, 0], ring[:, 1], order=3, obs_req=0) + fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) + axes[0].scatter(ring[:, 0], ring[:, 1], s=9, alpha=0.6) + axes[0].set_aspect("equal") + axes[0].set_title("Annulus: near-zero Pearson correlation, strong dependence") + partition_scatter(axes[1], ring_part, "NNS partition of annulus") + axes[1].set_aspect("equal") + save_figure(fig, OUT, "annulus_dependence.png") + + # %% [markdown] + # # p-values for NNS.dep() + subsection("p-values for NNS.dep()") + gap( + "Python nns_dep currently exposes the deterministic coefficient only; R also " + "accepts p.value=TRUE and print.map=TRUE. This vignette reproduces that executable " + "permutation analysis transparently rather than omitting it." + ) + rng = np.random.default_rng(123) + x_perm = np.arange(-5.0, 5.001, 0.1) + y_perm = x_perm**2 + rng.normal(size=x_perm.size) + observed = nns_dep(x_perm, y_perm) + repetitions = 30 if __import__("os").environ.get("NNS_VIGNETTE_FAST") == "1" else 100 + null_corr = np.empty(repetitions) + null_dep = np.empty(repetitions) + for index in range(repetitions): + permuted = rng.permutation(y_perm) + result = nns_dep(x_perm, permuted) + null_corr[index] = result["Correlation"] + null_dep[index] = result["Dependence"] + p_corr = (1.0 + np.sum(np.abs(null_corr) >= abs(observed["Correlation"]))) / (repetitions + 1.0) + p_dep = (1.0 + np.sum(null_dep >= observed["Dependence"])) / (repetitions + 1.0) + show("Observed NNS.dep", observed) + show("Permutation p-value for correlation", p_corr) + show("Permutation p-value for dependence", p_dep) + fig, axes = plt.subplots(1, 2, figsize=(11, 4.2)) + axes[0].hist(null_corr, bins=15) + axes[0].axvline(observed["Correlation"], linewidth=2) + axes[0].set_title("Permutation null: NNS correlation") + axes[1].hist(null_dep, bins=15) + axes[1].axvline(observed["Dependence"], linewidth=2) + axes[1].set_title("Permutation null: NNS dependence") + save_figure(fig, OUT, "permutation_map.png") + # %% [markdown] + # # Multivariate Dependence NNS.copula() + subsection("Multivariate Dependence NNS.copula()") 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)) + independent = rng.normal(size=(1000, 3)) + dependence = nns_copula(independent) + show("Three-variable NNS copula dependence", dependence) + gap( + "R's NNS.copula(..., plot=TRUE, independence.overlay=TRUE) renders its own map. " + "Python returns the scalar multivariate dependence; the figure below supplies a " + "pairwise projection and an independently simulated reference overlay." + ) + reference = rng.normal(size=(1000, 3)) + fig, axes = plt.subplots(1, 3, figsize=(13, 4)) + pairs = [(0, 1), (0, 2), (1, 2)] + for ax, (left, right) in zip(axes, pairs, strict=True): + ax.scatter(independent[:, left], independent[:, right], s=8, alpha=0.35, label="sample") + ax.scatter(reference[:, left], reference[:, right], s=8, alpha=0.15, label="independence overlay") + ax.set_title(f"X{left + 1} vs X{right + 1}") + axes[0].legend() + fig.suptitle(f"Multivariate NNS copula dependence = {dependence:.4f}") + save_figure(fig, OUT, "multivariate_copula.png") if __name__ == "__main__": From 4f99974af5a96e6b8f2c5aa03635ac8da791a80a Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:42:57 -0400 Subject: [PATCH 07/17] Port normalization and rescaling vignette --- .../04_normalization_and_rescaling.py | 183 +++++++++++++++++- 1 file changed, 177 insertions(+), 6 deletions(-) diff --git a/examples/vignettes/04_normalization_and_rescaling.py b/examples/vignettes/04_normalization_and_rescaling.py index 4b25b3e..9d405c2 100644 --- a/examples/vignettes/04_normalization_and_rescaling.py +++ b/examples/vignettes/04_normalization_and_rescaling.py @@ -1,11 +1,182 @@ -"""Canonical vignette 04: Normalization and rescaling. +"""04. Getting Started with NNS: Normalization and Rescaling. -Source of truth: - NNS/vignettes/NNSvignette_04_Normalization_and_Rescaling.Rmd - -The maintained implementation lives in ``normalization_rescaling.py``. +This is an instructional, section-for-section Python port of +``NNSvignette_04_Normalization_and_Rescaling.Rmd``. It preserves the R +vignette's distinction between multivariate normalization and univariate +rescaling, prints the important numerical structures, and writes equivalent +figures to ``examples/vignettes/output/04_normalization_and_rescaling``. """ -from normalization_rescaling import main +from __future__ import annotations + +import math +import sys +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +sys.path.insert(0, str(Path(__file__).resolve().parents[1])) +from _vignette_support import gap, output_dir, save_figure, section, show, subsection, table +from nns import nns_norm, nns_rescale + +OUT = output_dir(__file__) + + +def describe_columns(values: np.ndarray) -> list[tuple[str, float, float]]: + return [ + (f"column_{index + 1}", float(np.mean(values[:, index])), float(np.std(values[:, index], ddof=1))) + for index in range(values.shape[1]) + ] + + +def quantile_normalize(values: np.ndarray) -> np.ndarray: + """Small transparent implementation used only for the vignette comparison.""" + order = np.argsort(values, axis=0) + sorted_values = np.take_along_axis(values, order, axis=0) + rank_means = np.mean(sorted_values, axis=1) + output = np.empty_like(values) + for column in range(values.shape[1]): + output[order[:, column], column] = rank_means + return output + + +def main() -> None: + section("Overview") + print( + "NNS.norm rescales multiple variables to a common magnitude while preserving " + "their individual ordering and distributional form. NNS.rescale applies a " + "one-dimensional affine transformation to a requested range or expectation." + ) + + section("NNS.norm(): Normalize Multiple Variables") + print( + "For means m_j, the ratio grid is RG_ij = m_i / m_j. Linear normalization " + "averages RG by column; nonlinear normalization weights that grid by absolute " + "correlation for fewer than ten variables and NNS dependence otherwise." + ) + + subsection("Basic multivariate example") + rng = np.random.default_rng(123) + a = rng.normal(0, 1, 100) + b = rng.normal(0, 5, 100) + c = rng.normal(10, 1, 100) + d = rng.normal(10, 10, 100) + x = np.column_stack((a, b, c, d)) + + linear = np.asarray(nns_norm(x, linear=True), dtype=float) + nonlinear = np.asarray(nns_norm(x, linear=False), dtype=float) + + show("First six rows, linear normalization", linear[:6]) + table(["variable", "mean", "sample sd"], describe_columns(linear)) + print("Linear mode equalizes the column means exactly, as the R proof predicts.") + + show("First six rows, nonlinear normalization", nonlinear[:6]) + table(["variable", "mean", "sample sd"], describe_columns(nonlinear)) + print( + "Nonlinear means differ because each is a dependence-weighted average of the " + "original means rather than the single grand mean." + ) + + subsection("Unequal vector lengths") + vectors = [rng.normal(0, 1, 10), rng.normal(5, 5, 5), rng.normal(10, 10, 8)] + unequal = nns_norm(vectors) + show("Normalized unequal-length vectors", {f"x_{i + 1}": v for i, v in enumerate(unequal)}) + print( + "As in R, unequal lengths force the linear scaling path because a dependence " + "matrix requires aligned observations." + ) + + subsection("Quantile normalization comparison") + qnorm = quantile_normalize(x) + print( + "Quantile normalization makes every sorted column identical. NNS normalization " + "does the opposite: it aligns magnitudes while retaining each variable's shape." + ) + show("Sorted quantile-normalized columns (first five ranks)", np.sort(qnorm, axis=0)[:5]) + + section("NNS.rescale(): Distribution Rescaling") + subsection("Min-max scaling") + raw = np.array([-2.5, 0.2, 1.1, 3.7, 5.0]) + scaled = nns_rescale(raw, a=5, b=10, method="minmax") + table(["raw", "scaled"], zip(raw, scaled, strict=True)) + show("Scaled range", np.array([np.min(scaled), np.max(scaled)])) + + subsection("Risk-neutral terminal scaling") + rng = np.random.default_rng(123) + s0, rate, maturity = 100.0, 0.05, 1.0 + prices = s0 * np.exp(np.cumsum(rng.normal(0.0005, 0.02, 250))) + terminal = nns_rescale( + prices, + a=s0, + b=rate, + method="riskneutral", + time_to_maturity=maturity, + type="Terminal", + ) + target = s0 * math.exp(rate * maturity) + table( + ["quantity", "value"], + [ + ("mean_original", float(np.mean(prices))), + ("mean_rescaled", float(np.mean(terminal))), + ("target", target), + ], + ) + + subsection("Risk-neutral discounted scaling") + discounted = nns_rescale( + prices, + a=s0, + b=rate, + method="riskneutral", + time_to_maturity=maturity, + type="Discounted", + ) + table( + ["quantity", "value"], + [ + ("mean_returned_series", float(np.mean(discounted))), + ("target_discounted_mean", s0), + ], + ) + gap( + "Python spells the R argument T as time_to_maturity. The statistical operation " + "is the same; only the public parameter name differs." + ) + + section("Distribution-shape comparison") + rng = np.random.default_rng(123) + x1 = rng.normal(5, 2, 1000) + x2 = rng.gamma(3, 1, 1000) + pair = np.column_stack((x1, x2)) + pair_lin = np.asarray(nns_norm(pair, linear=True)) + pair_nonlin = np.asarray(nns_norm(pair, linear=False)) + pair_minmax = np.column_stack( + ((pair[:, i] - pair[:, i].min()) / np.ptp(pair[:, i]) for i in range(pair.shape[1])) + ) + + fig, axes = plt.subplots(2, 2, figsize=(11, 8)) + panels = [ + (pair, "Original Variables"), + (pair_lin, "NNS.norm(linear=True)"), + (pair_nonlin, "NNS.norm(linear=False)"), + (pair_minmax, "Standard Min-Max"), + ] + for ax, (data, title) in zip(axes.ravel(), panels, strict=True): + shared = np.histogram_bin_edges(data.ravel(), bins=15) + ax.hist(data[:, 0], bins=shared, alpha=0.45, label="Normal") + ax.hist(data[:, 1], bins=shared, alpha=0.45, label="Gamma") + ax.set_title(title) + ax.legend() + save_figure(fig, OUT, "normalization_comparison.png") + + section("Conceptual summary") + print( + "NNS.norm is multivariate and dependence-aware. NNS.rescale is univariate and " + "range- or expectation-targeted. Both are monotone affine transformations and " + "therefore preserve rank structure for subsequent copula and dependence work." + ) + if __name__ == "__main__": main() From 8f511164d1849103aa1823f4e2ed472e72340939 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:57:30 -0400 Subject: [PATCH 08/17] Port clustering and regression vignette section by section --- .../vignettes/07_clustering_and_regression.py | 46 +++---------------- 1 file changed, 6 insertions(+), 40 deletions(-) diff --git a/examples/vignettes/07_clustering_and_regression.py b/examples/vignettes/07_clustering_and_regression.py index be465a8..6c3b724 100644 --- a/examples/vignettes/07_clustering_and_regression.py +++ b/examples/vignettes/07_clustering_and_regression.py @@ -1,45 +1,11 @@ -"""Canonical vignette 07: Clustering and regression. +"""07. Getting Started with NNS: Clustering and Regression. -Source of truth: +Instructional Python port of: NNS/vignettes/NNSvignette_07_Clustering_and_Regression.Rmd + +Run directly to print the important returned structures and write the figures to +``examples/vignettes/output/07_clustering_and_regression``. """ 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() +import \ No newline at end of file From 75045e4785e8af1794f3879d500946e626520d17 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:58:34 -0400 Subject: [PATCH 09/17] Make instructional vignettes directly executable --- examples/vignettes/sitecustomize.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) create mode 100644 examples/vignettes/sitecustomize.py diff --git a/examples/vignettes/sitecustomize.py b/examples/vignettes/sitecustomize.py new file mode 100644 index 0000000..169f13b --- /dev/null +++ b/examples/vignettes/sitecustomize.py @@ -0,0 +1,14 @@ +"""Allow numbered vignette scripts to import shared support from the repository. + +Python imports ``sitecustomize`` from the script directory during startup. By +adding the repository root here, each vignette remains directly executable as +``python examples/vignettes/NN_topic.py`` without packaging the examples. +""" +from __future__ import annotations + +import sys +from pathlib import Path + +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) From 20e34d3219a254bb9ce81f9298a203dfd3ec69e7 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:59:02 -0400 Subject: [PATCH 10/17] Port classification vignette section by section --- examples/vignettes/08_classification.py | 250 ++++++++++++++++++++---- 1 file changed, 214 insertions(+), 36 deletions(-) diff --git a/examples/vignettes/08_classification.py b/examples/vignettes/08_classification.py index 3d1f3e0..6c0a4c5 100644 --- a/examples/vignettes/08_classification.py +++ b/examples/vignettes/08_classification.py @@ -1,72 +1,250 @@ -"""Canonical vignette 08: Classification. +"""08. Getting Started with NNS: Classification. -Source of truth: +Instructional Python port of: 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. +The examples use the exact R ``iris`` observations and preserve the vignette's +order: classification geometry, partitions, boosting, stacking, and parameter +interpretation. Figures are written to +``examples/vignettes/output/08_classification``. """ from __future__ import annotations +# %% [markdown] +# # Classification +# +# NNS classification uses the same partition-and-regression mechanism as the +# continuous model, with the response encoded as classes beginning at 1. The R +# vignette emphasizes that these are multidimensional partitions rather than a +# sequence of one-variable tree splits. + +import matplotlib.pyplot as plt import numpy as np -from nns import nns_boost, nns_reg, nns_stack +from examples._vignette_support import ( + fast_mode, + gap, + load_iris, + note, + output_dir, + save_figure, + section, + show, + subsection, + table, +) +from nns import nns_boost, nns_part, nns_reg, nns_stack + +OUT = output_dir(__file__) -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 accuracy(predicted: np.ndarray, actual: np.ndarray) -> float: + return float(np.mean(np.asarray(predicted, dtype=float) == np.asarray(actual, dtype=float))) def main() -> None: - x_train, y_train, x_test, y_test = _three_class_data() + iris_x, iris_y, levels = load_iris() + feature_names = ["Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"] + test_idx = np.arange(140, 150) + train_idx = np.arange(0, 140) + x_train, y_train = iris_x[train_idx], iris_y[train_idx] + x_test, y_test = iris_x[test_idx], iris_y[test_idx] + + section("Classification") + note( + "The response is encoded 1=setosa, 2=versicolor, 3=virginica. " + "This preserves the R contract that the base category is 1, not 0." + ) + table( + ["code", "class", "training rows", "test rows"], + [ + (code, label, int(np.sum(y_train == code)), int(np.sum(y_test == code))) + for code, label in enumerate(levels, start=1) + ], + ) + + # %% [markdown] + # ## Splits versus partitions + # + # A tree asks one split question at a time. NNS partitions the joint support + # through partial-moment quadrants, then uses the resulting regional centers + # as regression points. The figure below shows the first two iris predictors + # only so that the otherwise four-dimensional geometry can be inspected. + subsection("Splits versus partitions") + pair = nns_part(iris_x[:, 0], iris_x[:, 1], order=4, obs_req=0) + labels = np.asarray(pair["dt"]["quadrant"], dtype=str) + _, codes = np.unique(labels, return_inverse=True) + fig, ax = plt.subplots(figsize=(7, 5)) + ax.scatter(iris_x[:, 0], iris_x[:, 1], c=codes, cmap="tab20", s=28, alpha=0.75) + rp = pair["regression.points"] + ax.scatter(rp["x"], rp["y"], marker="x", s=70, linewidths=2, label="partition centers") + ax.set_xlabel(feature_names[0]) + ax.set_ylabel(feature_names[1]) + ax.set_title("NNS joint partitions on two iris dimensions") + ax.legend() + save_figure(fig, OUT, "01_iris_joint_partitions.png") + show("Two-dimensional partition result", pair) - reg = nns_reg(x_train, y_train, type="CLASS", point_est=x_test) + # %% [markdown] + # ## NNS partitions used by the classifier + # + # The R vignette prints ``NNS.reg(... )$rhs.partitions`` for all four + # predictors. Python currently returns the fitted reduction and equation but + # not the internal per-predictor partition list. We therefore expose the + # closest public diagnostic explicitly instead of pretending the object is + # present. + subsection("NNS partitions") + base_fit = nns_reg( + iris_x, + iris_y, + type="CLASS", + point_est=iris_x[:10], + residual_plot=False, + class_levels=levels, + ) + show("Classification regression result", base_fit) + gap( + "R exposes NNS.reg(... )$rhs.partitions for every predictor. The Python " + "public result does not yet expose that internal list. The two-dimensional " + "NNS.part result above is the supported public partition diagnostic." + ) + + # %% [markdown] + # ## NNS.boost() + # + # This is the exact R holdout: rows 141--150 of iris. The same 10 epochs, + # 10 learner trials, class balancing, and seed are used by default. + subsection("NNS.boost") boost = nns_boost( x_train, y_train, x_test, type="CLASS", - epochs=5, - learner_trials=10, - cv_size=0.25, + epochs=3 if fast_mode() else 10, + learner_trials=4 if fast_mode() else 10, balance=True, status=False, seed=123, + class_levels=levels, ) + boost_pred = np.asarray(boost["results"], dtype=float) + show("NNS.boost result", boost) + print(f"Boost holdout accuracy: {accuracy(boost_pred, y_test):.4f}") + + weights = boost["feature.weights"] + frequencies = boost["feature.frequency"] + weight_values = np.asarray(list(weights.values()), dtype=float) + frequency_values = np.asarray(list(frequencies.values()), dtype=float) + labels_for_plot = list(weights.keys()) + fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) + axes[0].bar(np.arange(len(weight_values)), weight_values) + axes[0].set_xticks(np.arange(len(weight_values)), labels_for_plot, rotation=35, ha="right") + axes[0].set_title("Boost feature weights") + axes[1].bar(np.arange(len(frequency_values)), frequency_values) + axes[1].set_xticks(np.arange(len(frequency_values)), list(frequencies.keys()), rotation=35, ha="right") + axes[1].set_title("Boost feature frequency") + save_figure(fig, OUT, "02_boost_feature_diagnostics.png") + + fig, ax = plt.subplots(figsize=(8, 4.5)) + positions = np.arange(y_test.size) + ax.plot(positions, y_test, marker="o", label="actual class") + ax.plot(positions, boost_pred, marker="x", linestyle="--", label="boost prediction") + ax.set_yticks([1, 2, 3], levels) + ax.set_xlabel("Held-out iris observation") + ax.set_title("NNS.boost: exact R vignette holdout") + ax.legend() + save_figure(fig, OUT, "03_boost_holdout_predictions.png") + note( + "Balanced boosting is stochastic. Python reproduces R's random stream on " + "deterministic and simple sampling paths, but the interleaved class-balance " + "draw order is a documented remaining bit-for-bit RNG gap. The statistical " + "procedure and returned diagnostics are the same." + ) + + # %% [markdown] + # ## Cross-validation classification using NNS.stack() + # + # The R vignette next combines the full and dimension-reduced classifiers by + # cross-validating ``n.best``, the reduction threshold, and the ensemble + # weights. The default here preserves its five-fold specification; CI fast + # mode uses one fold solely to reduce runtime. + subsection("Cross-validation classification using NNS.stack") stack = nns_stack( x_train, y_train, x_test, type="CLASS", balance=True, - folds=1, - cv_size=0.25, + folds=1 if fast_mode() else 5, + method=(1, 2), status=False, seed=123, ) + show("NNS.stack result", stack) + stack_pred = np.asarray(stack["stack"], dtype=float) + reg_pred = np.asarray(stack["reg"], dtype=float) + dim_pred = np.asarray(stack["dim.red"], dtype=float) + table( + ["component", "accuracy"], + [ + ("full regression", accuracy(reg_pred, y_test)), + ("dimension reduction", accuracy(dim_pred, y_test)), + ("stack", accuracy(stack_pred, y_test)), + ], + ) - 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"]) + fig, ax = plt.subplots(figsize=(9, 4.8)) + ax.plot(positions, y_test, marker="o", linewidth=2, label="actual") + ax.plot(positions, reg_pred, marker="x", linestyle="--", label="reg") + ax.plot(positions, dim_pred, marker="s", linestyle=":", label="dim.red") + ax.plot(positions, stack_pred, marker="d", linestyle="-.", label="stack") + ax.set_yticks([1, 2, 3], levels) + ax.set_xlabel("Held-out iris observation") + ax.set_title("Cross-validated NNS classification components") + ax.legend(ncol=4) + save_figure(fig, OUT, "04_stack_holdout_predictions.png") + + # %% [markdown] + # ## Depth, nearest-neighbor classification, and extremes + # + # The R vignette closes by explaining three controls: maximum partition depth, + # ``n.best = 1`` for a donor/nearest-neighbor classifier, and ``extreme`` in + # boosting. Each supported path is demonstrated directly. + subsection("Depth, n.best = 1, and extreme boosting") + nearest = nns_reg( + x_train, + y_train, + type="CLASS", + order="max", + n_best=1, + point_est=x_test, + class_levels=levels, + ) + nearest_pred = np.asarray(nearest["Point.est"], dtype=float) + extreme = nns_boost( + x_train, + y_train, + x_test, + type="CLASS", + depth="max", + learner_trials=3 if fast_mode() else 10, + epochs=2 if fast_mode() else 10, + balance=True, + extreme=True, + status=False, + seed=123, + class_levels=levels, + ) + extreme_pred = np.asarray(extreme["results"], dtype=float) + table( + ["model", "accuracy", "predictions"], + [ + ("order=max, n.best=1", accuracy(nearest_pred, y_test), nearest_pred.astype(int)), + ("boost extreme=True", accuracy(extreme_pred, y_test), extreme_pred.astype(int)), + ], + ) + show("Nearest-neighbor classification result", nearest) + show("Extreme boost result", extreme) if __name__ == "__main__": From f42255b9920a288dad7fa497911b0a6fe8d4244d Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 00:59:38 -0400 Subject: [PATCH 11/17] Port forecasting vignette section by section --- examples/vignettes/09_forecasting.py | 294 ++++++++++++++++++++++++++- 1 file changed, 290 insertions(+), 4 deletions(-) diff --git a/examples/vignettes/09_forecasting.py b/examples/vignettes/09_forecasting.py index c36c155..a3d8817 100644 --- a/examples/vignettes/09_forecasting.py +++ b/examples/vignettes/09_forecasting.py @@ -1,11 +1,297 @@ -"""Canonical vignette 09: Forecasting. +"""09. Getting Started with NNS: Forecasting. -Source of truth: +Instructional Python port of: NNS/vignettes/NNSvignette_09_Forecasting.Rmd -The maintained implementation lives in ``time_series_forecasting.py``. +The exact R ``AirPassengers`` series, training window, holdout, seasonal-factor +search, and optimizer sequence are preserved. Figures are written to +``examples/vignettes/output/09_forecasting``. """ -from time_series_forecasting import main +from __future__ import annotations + +# %% [markdown] +# # Forecasting +# +# NNS forecasting treats a time series as a nonlinear regression problem on +# lagged seasonal components. The R vignette moves from a fixed seasonal factor +# to validation-based selection, automatic seasonality, and finally +# ``NNS.ARMA.optim``. + +import matplotlib.pyplot as plt +import numpy as np + +from examples._vignette_support import ( + fast_mode, + gap, + load_air_passengers, + note, + output_dir, + save_figure, + section, + show, + subsection, + table, +) +from nns import nns_arma, nns_arma_optim, nns_seas + +OUT = output_dir(__file__) + + +def rmse(predicted: np.ndarray, actual: np.ndarray) -> float: + predicted = np.asarray(predicted, dtype=float) + actual = np.asarray(actual, dtype=float) + return float(np.sqrt(np.mean((predicted - actual) ** 2))) + + +def estimates(result: object) -> np.ndarray: + if isinstance(result, dict): + for key in ("Estimates", "results"): + if key in result: + return np.asarray(result[key], dtype=float) + raise KeyError("forecast dictionary has no Estimates/results field") + return np.asarray(result, dtype=float) + + +def forecast_figure( + history: np.ndarray, + actual: np.ndarray | None, + predicted: np.ndarray, + title: str, + filename: str, + lower: np.ndarray | None = None, + upper: np.ndarray | None = None, +) -> None: + fig, ax = plt.subplots(figsize=(10, 5)) + x_hist = np.arange(history.size) + x_fc = np.arange(history.size, history.size + predicted.size) + ax.plot(x_hist, history, label="history") + if actual is not None: + ax.plot(x_fc[: actual.size], actual, label="actual holdout") + ax.plot(x_fc, predicted, label="NNS forecast", linewidth=2) + if lower is not None and upper is not None: + ax.fill_between(x_fc, lower, upper, alpha=0.2, label="prediction interval") + ax.set_xlabel("Month index") + ax.set_ylabel("Passengers") + ax.set_title(title) + ax.legend() + save_figure(fig, OUT, filename) + + +def main() -> None: + air = load_air_passengers() + training_set = 100 + horizon = air.size - training_set + training = air[:training_set] + actual = air[training_set:] + + section("Forecasting") + print(f"AirPassengers observations: {air.size}") + print(f"Training observations: {training_set}; validation observations: {horizon}") + + # %% [markdown] + # ## Linear regression + # + # The R vignette begins with a 12-month seasonal factor and a 44-month + # validation horizon. ``method='lin'`` regresses each seasonal component + # linearly before recombining the forecasts. + subsection("Linear regression") + linear_raw = nns_arma( + air, + h=horizon, + training_set=training_set, + method="lin", + seasonal_factor=12, + plot=False, + seasonal_plot=False, + ) + linear = estimates(linear_raw) + linear_rmse = rmse(linear, actual) + print(f"Linear seasonal-factor-12 RMSE: {linear_rmse:.6f}") + show("Linear forecast", linear_raw) + forecast_figure(training, actual, linear, "Linear NNS.ARMA, seasonal factor 12", "01_linear.png") + + # %% [markdown] + # ## Nonlinear regression + # + # The same data and seasonal factor are now passed through the nonlinear NNS + # regression engine. The comparison isolates the regression form; nothing + # else about the experiment changes. + subsection("Nonlinear regression") + nonlinear_raw = nns_arma( + air, + h=horizon, + training_set=training_set, + method="nonlin", + seasonal_factor=12, + plot=False, + seasonal_plot=False, + ) + nonlinear = estimates(nonlinear_raw) + nonlinear_rmse = rmse(nonlinear, actual) + print(f"Nonlinear seasonal-factor-12 RMSE: {nonlinear_rmse:.6f}") + show("Nonlinear forecast", nonlinear_raw) + forecast_figure( + training, + actual, + nonlinear, + "Nonlinear NNS.ARMA, seasonal factor 12", + "02_nonlinear.png", + ) + table( + ["method", "seasonal factor", "RMSE"], + [("linear", 12, linear_rmse), ("nonlinear", 12, nonlinear_rmse)], + ) + + # %% [markdown] + # ## Cross-validation of the seasonal factor + # + # Following the R vignette, every period from 1 through 25 is evaluated on + # the same 44-observation holdout. CI fast mode uses a reduced grid only for + # runtime; normal execution preserves the complete search. + subsection("Cross-validation") + periods = np.arange(1, 26) if not fast_mode() else np.array([1, 6, 12, 18, 24]) + rows: list[tuple[int, float]] = [] + forecasts: dict[int, np.ndarray] = {} + for period in periods: + raw = nns_arma( + air, + h=horizon, + training_set=training_set, + method="nonlin", + seasonal_factor=int(period), + plot=False, + seasonal_plot=False, + ) + pred = estimates(raw) + forecasts[int(period)] = pred + rows.append((int(period), rmse(pred, actual))) + table(["Period", "RMSE"], rows) + best_period, best_rmse = min(rows, key=lambda item: item[1]) + print(f"Best validated period: {best_period}; RMSE={best_rmse:.6f}") + + fig, ax = plt.subplots(figsize=(9, 4.8)) + ax.plot([row[0] for row in rows], [row[1] for row in rows], marker="o") + ax.axvline(best_period, linestyle="--", label=f"best={best_period}") + ax.set_xlabel("Seasonal factor") + ax.set_ylabel("Validation RMSE") + ax.set_title("Cross-validation of NNS.ARMA seasonal factor") + ax.legend() + save_figure(fig, OUT, "03_period_cross_validation.png") + forecast_figure( + training, + actual, + forecasts[best_period], + f"Best nonlinear validation forecast, period {best_period}", + "04_best_period_forecast.png", + ) + + # %% [markdown] + # ## Automatic seasonality detection + # + # ``NNS.seas`` tests candidate component periods against the variability of + # the original series. The R vignette constrains candidates around a modulo + # of 12 and displays the ranked period table. + subsection("NNS.seas") + seasonality = nns_seas(air, modulo=12, mod_only=False, plot=False) + show("NNS.seas result", seasonality) + all_periods = seasonality["all.periods"] + fig, ax = plt.subplots(figsize=(9, 4.8)) + ax.scatter(all_periods["Period"], all_periods["Coefficient.of.Variation"], label="component CV") + ax.axhline( + float(all_periods["Variable.Coefficient.of.Variation"][0]), + linestyle="--", + label="series CV", + ) + ax.set_xlabel("Period") + ax.set_ylabel("Coefficient of variation") + ax.set_title("NNS seasonality diagnostics") + ax.legend() + save_figure(fig, OUT, "05_seasonality.png") + + # %% [markdown] + # ## Cross-validating combinations with NNS.ARMA.optim + # + # The optimizer compares linear, nonlinear, and combined forecasts over + # candidate seasonal-factor combinations. The exact R candidate grid is + # ``12, 18, ..., 60``; fast mode reduces it solely for CI. + subsection("NNS.ARMA.optim") + candidate_periods = list(range(12, 61, 6)) if not fast_mode() else [12, 18, 24] + optim = nns_arma_optim( + air, + training_set=training_set, + seasonal_factor=candidate_periods, + obj_fn=rmse, + objective="min", + pred_int=0.95, + print_trace=not fast_mode(), + plot=False, + ) + show("NNS.ARMA.optim validation result", optim) + optim_pred = np.asarray(optim["results"], dtype=float) + print(f"Optimizer holdout RMSE: {rmse(optim_pred, actual):.6f}") + forecast_figure( + training, + actual, + optim_pred, + "NNS.ARMA.optim validation forecast", + "06_arma_optim_validation.png", + np.asarray(optim["lower.pred.int"], dtype=float), + np.asarray(optim["upper.pred.int"], dtype=float), + ) + + # %% [markdown] + # ## Extension estimates + # + # Once the validation specification is selected, the R vignette refits on the + # complete series and extends it 50 months. Prediction intervals are retained + # and displayed rather than reduced to a printed forecast vector. + subsection("Extension estimates") + extension_h = 12 if fast_mode() else 50 + extension = nns_arma_optim( + air, + h=extension_h, + seasonal_factor=candidate_periods, + obj_fn=rmse, + objective="min", + pred_int=0.95, + print_trace=not fast_mode(), + plot=False, + ) + show("NNS.ARMA.optim extension result", extension) + forecast_figure( + air, + None, + np.asarray(extension["results"], dtype=float), + f"AirPassengers extension: {extension_h} months", + "07_extension_forecast.png", + np.asarray(extension["lower.pred.int"], dtype=float), + np.asarray(extension["upper.pred.int"], dtype=float), + ) + + # %% [markdown] + # ## Parameter interpretation and multivariate forecasting + # + # The R vignette closes with parameter notes and points readers to NNS.VAR for + # multivariate systems. Python exposes ``nns_var`` but the canonical R + # forecasting vignette contains no executable VAR example to translate here. + subsection("Parameter interpretation") + note( + "seasonal_factor selects lag components; method chooses linear, nonlinear, " + "or combined regression; weights combine multiple periods; dynamic=True " + "re-estimates periods recursively; shrink averages linear estimates toward " + "component means; pred_int requests partial-moment prediction intervals." + ) + gap( + "R's NNS.ARMA seasonal.plot side effect is not separately exposed in the " + "Python return object. The translated seasonality figure above is generated " + "from NNS.seas diagnostics instead." + ) + note( + "NNS.VAR is available in Python, but the R vignette only references the " + "multivariate method and paper; it does not provide an executable example. " + "No synthetic VAR demonstration is inserted into the canonical port." + ) + if __name__ == "__main__": main() From 795678a50661eb2a02f03c3193851a6ca221eede Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 01:00:03 -0400 Subject: [PATCH 12/17] Validate instructional vignette execution and figures --- tests/docs/test_vignette_examples.py | 35 ++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/tests/docs/test_vignette_examples.py b/tests/docs/test_vignette_examples.py index f67af5c..f1aeda9 100644 --- a/tests/docs/test_vignette_examples.py +++ b/tests/docs/test_vignette_examples.py @@ -1,18 +1,43 @@ from __future__ import annotations +import os import subprocess import sys from pathlib import Path import pytest -EXAMPLE_DIR = Path(__file__).resolve().parents[2] / "examples" / "vignettes" +ROOT = Path(__file__).resolve().parents[2] +EXAMPLE_DIR = ROOT / "examples" / "vignettes" +CANONICAL = sorted(EXAMPLE_DIR.glob("[0-9][0-9]_*.py")) -@pytest.mark.parametrize("script", sorted(EXAMPLE_DIR.glob("*.py")), ids=lambda p: p.stem) -def test_vignette_example(script: Path) -> None: - subprocess.run( +@pytest.mark.parametrize("script", CANONICAL, ids=lambda p: p.stem) +def test_canonical_vignette_executes_and_generates_figures(script: Path) -> None: + output = EXAMPLE_DIR / "output" / script.stem + if output.exists(): + for figure in output.glob("*.png"): + figure.unlink() + + env = os.environ.copy() + env.update( + { + "NNS_VIGNETTE_FAST": "1", + "MPLBACKEND": "Agg", + "PYTHONPATH": str(ROOT), + } + ) + completed = subprocess.run( [sys.executable, str(script)], check=True, - cwd=Path(__file__).resolve().parents[2], + cwd=ROOT, + env=env, + capture_output=True, + text=True, + timeout=240, ) + + assert "Saved figure:" in completed.stdout + figures = sorted(output.glob("*.png")) + assert figures, f"{script.name} produced no figures" + assert all(path.stat().st_size > 0 for path in figures) From 6c4606aa07a43e13985081065634f50c71ea78a0 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 01:00:15 -0400 Subject: [PATCH 13/17] Enforce instructional vignette parity contract --- .../docs/test_canonical_vignette_manifest.py | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/tests/docs/test_canonical_vignette_manifest.py b/tests/docs/test_canonical_vignette_manifest.py index 3359bf0..d44e716 100644 --- a/tests/docs/test_canonical_vignette_manifest.py +++ b/tests/docs/test_canonical_vignette_manifest.py @@ -5,6 +5,7 @@ ROOT = Path(__file__).resolve().parents[2] EXAMPLE_DIR = ROOT / "examples" / "vignettes" +DATA_DIR = ROOT / "examples" / "data" EXPECTED = { "01": ("NNSvignette_01_Overview.Rmd", "01_overview.py"), "02": ("NNSvignette_02_Partial_Moments.Rmd", "02_partial_moments.py"), @@ -23,13 +24,42 @@ def test_canonical_vignette_set_is_complete_and_ordered() -> None: assert numbered == [python for _, python in EXPECTED.values()] -def test_each_python_entrypoint_names_its_r_source() -> None: +def test_each_python_vignette_is_instructional_not_a_wrapper() -> None: for r_source, python_name in EXPECTED.values(): text = (EXAMPLE_DIR / python_name).read_text(encoding="utf-8") assert r_source in text + assert "# %% [markdown]" in text + assert "save_figure(" in text + assert "def main()" in text + assert len(text.splitlines()) >= 90 + assert "from overview import main" not in text + assert "from time_series_forecasting import main" not 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()) + for criterion in ( + "section_order", + "same_data", + "figures", + "returned_structures", + "explicit_gaps", + ): + assert criterion in text + + +def test_exact_r_datasets_are_committed_locally() -> None: + expected_rows = {"iris.csv": 151, "mtcars.csv": 33, "AirPassengers.csv": 145} + for filename, line_count in expected_rows.items(): + path = DATA_DIR / filename + assert path.exists() + assert len(path.read_text(encoding="utf-8").splitlines()) == line_count + + +def test_parity_ledger_covers_all_nine_vignettes() -> None: + text = (EXAMPLE_DIR / "PARITY.md").read_text(encoding="utf-8") + for vignette_id, (_, python_name) in EXPECTED.items(): + assert f"## {vignette_id}" in text + assert python_name in text From 2cf343682949a34b19df76f04cb1f30082c276c6 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 01:00:26 -0400 Subject: [PATCH 14/17] Define full vignette parity acceptance criteria --- examples/vignettes/manifest.yml | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/examples/vignettes/manifest.yml b/examples/vignettes/manifest.yml index d948f50..e016028 100644 --- a/examples/vignettes/manifest.yml +++ b/examples/vignettes/manifest.yml @@ -2,7 +2,26 @@ canonical_source: repository: OVVO-Financial/NNS branch: NNS-Beta-Version directory: vignettes - policy: R defines the statistical narrative, datasets, and section order. + policy: R defines the statistical narrative, datasets, section order, and interpretation. + +acceptance_criteria: + section_order: Preserve the R vignette section order and statistical narrative. + demonstrations: Reproduce every executable R demonstration supported by Python. + same_data: Use the same R datasets or committed faithful equivalents. + figures: Generate equivalent figures instead of replacing graphics with console summaries. + returned_structures: Display important returned structures and numerical results. + explicit_gaps: Label genuine R/Python API gaps at the point where the R section uses them. + instructional_form: Numbered Python files are Jupytext-style instructional scripts, not smoke tests. + +runtime: + default: Full R-equivalent demonstrations and parameter grids. + fast_mode: Set NNS_VIGNETTE_FAST=1 only for CI runtime reduction; section coverage is unchanged. + figures: examples/vignettes/output// + +datasets: + - examples/data/iris.csv + - examples/data/mtcars.csv + - examples/data/AirPassengers.csv examples: - id: '01' From 05a9876263a8e29bf5ddf58405d0f7f54e27f4e9 Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 01:00:46 -0400 Subject: [PATCH 15/17] Add section-by-section vignette parity ledger --- examples/vignettes/PARITY.md | 101 +++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 examples/vignettes/PARITY.md diff --git a/examples/vignettes/PARITY.md b/examples/vignettes/PARITY.md new file mode 100644 index 0000000..60b28af --- /dev/null +++ b/examples/vignettes/PARITY.md @@ -0,0 +1,101 @@ +# Canonical R-to-Python vignette parity + +The R package is the source of truth. This ledger records what each Python +instructional script reproduces and names remaining API differences explicitly. +A missing Python facility is never replaced silently with an unrelated example. + +## 01 — `01_overview.py` + +Preserves the overview curriculum: partial-moment foundations, descriptive and +distributional tools, nonlinear dependence, normalization/rescaling, ANOVA and +stochastic superiority, regression/boosting/stacking/causality, forecasting, +maximum-entropy simulation, and stochastic dominance. Figures accompany the +major statistical demonstrations. + +**Known gaps:** interactive R plotting devices are rendered as saved Matplotlib +figures. This changes presentation, not the statistical calculation. + +## 02 — `02_partial_moments.py` + +Reproduces the mean, variance, standard deviation, higher-moment, mode, +covariance, covariance-matrix, correlation, empirical/continuous CDF, copula, +numerical-integration, and Bayes demonstrations in R order. Co-partial matrices +and reconstructed covariance are printed, and CDF/copula figures are saved. + +**Known gaps:** R's interactive 3-D `rgl` views are translated to static +Matplotlib 3-D figures. + +## 03 — `03_correlation_and_dependence.py` + +Reproduces linear, quadratic, cubic, sinusoidal, circular, and discrete examples; +contrasts Pearson correlation with NNS correlation/dependence; evaluates +asymmetry, copula dependence, and causation; and saves relationship and matrix +figures. + +**Known gaps:** Python's public dependence API returns pairwise scalar results; +matrix displays are assembled transparently from those public pair calls. + +## 04 — `04_normalization_and_rescaling.py` + +Preserves linear and nonlinear normalization, affine min/max rescaling, +risk-neutral terminal and discounted transformations, and the R interpretation +of shape-preserving transformations. Input/output distributions are plotted. + +**Known gaps:** R chart-type side effects are replaced with explicit saved +figures generated from returned arrays. + +## 05 — `05_sampling_and_simulation.py` + +Preserves empirical-CDF inversion across partial-moment degrees, +maximum-entropy bootstrap sampling, dependence-targeted Monte Carlo sampling, +and ensemble diagnostics. Distribution and path figures are generated. + +**Known gaps:** stochastic draws use NumPy's local generator, so simulations are +distributionally equivalent but not generally bit-for-bit equal to R draws. + +## 06 — `06_comparing_distributions.py` + +Preserves two-sample and multi-group NNS ANOVA, stochastic superiority, +first/second/third-order stochastic dominance, efficient-set extraction, and +stochastic-dominance clustering. CDF, superiority, and cluster figures are +saved, while the full returned structures are printed. + +**Known gaps:** R's ANOVA plotting side effect is reconstructed from the public +returned statistics and empirical CDFs. + +## 07 — `07_clustering_and_regression.py` + +Preserves partition orders 1–4, X-only partitions, clusters used in regression, +univariate and full-grid multivariate regression, interpolation/extrapolation, +dimension reduction, thresholding, classification, `NNS.stack`, duplicated +predictor dimensions, smoothing, and univariate/multivariate imputation. The +exact iris data are committed locally. Six figure groups and important model +objects are produced. + +**Known gaps:** Python does not yet expose R's Voronoi tessellation side effect or +`NNS.reg(... )$rhs.partitions`; both are labeled in the relevant sections and +supported public partition diagnostics are shown instead. + +## 08 — `08_classification.py` + +Preserves the splits-versus-partitions narrative, exact iris holdout rows +141–150, `NNS.boost` with balancing and the R parameterization, +`NNS.stack` cross-validation, and the depth/nearest-neighbor/extreme controls. +Feature diagnostics, holdout paths, and returned model structures are displayed. + +**Known gaps:** Python does not yet expose `rhs.partitions`. Balanced boosting +has a documented bit-for-bit RNG ordering gap from R's interleaved class draws; +the procedure and diagnostics remain equivalent. + +## 09 — `09_forecasting.py` + +Preserves the exact AirPassengers series, 100-observation training set, +44-observation validation set, linear/nonlinear seasonal-factor-12 comparison, +periods 1–25 validation, `NNS.seas`, `NNS.ARMA.optim` over periods 12–60 by 6, +and the 50-period extension with prediction intervals. Every forecast stage is +plotted and its returned structure displayed. + +**Known gaps:** Python does not separately expose R's `seasonal.plot` side +effect, so the equivalent diagnostic is generated from `NNS.seas`. The R +vignette mentions `NNS.VAR` but contains no executable VAR demonstration; the +Python canonical port therefore does not invent one. From 6fe73da77927d508c604fa6919be425097397a4d Mon Sep 17 00:00:00 2001 From: OVVO-Financial Date: Thu, 16 Jul 2026 01:01:05 -0400 Subject: [PATCH 16/17] Document full instructional vignette contract --- examples/vignettes/README.md | 44 ++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/examples/vignettes/README.md b/examples/vignettes/README.md index f5f0a1a..e77e48a 100644 --- a/examples/vignettes/README.md +++ b/examples/vignettes/README.md @@ -1,11 +1,11 @@ -# Canonical NNS examples +# Canonical NNS instructional vignettes -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. +The R package is the source of truth for the NNS curriculum. Its nine numbered +vignettes define the statistical narrative, datasets, section order, examples, +and interpretation. The numbered Python files are section-by-section +instructional ports—not wrappers, API samplers, or smoke tests. -| # | Canonical topic | R source | Python entry point | +| # | Canonical topic | R source | Python vignette | |---|---|---|---| | 01 | Overview | `NNSvignette_01_Overview.Rmd` | `01_overview.py` | | 02 | Partial Moments | `NNSvignette_02_Partial_Moments.Rmd` | `02_partial_moments.py` | @@ -17,18 +17,38 @@ language-appropriate containers and syntax. | 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. +Each script: + +- preserves the R section order and explanatory progression; +- uses the same datasets, with exact local copies of `iris`, `mtcars`, and + `AirPassengers` committed under `examples/data`; +- reproduces every executable demonstration supported by the Python API; +- saves equivalent figures under `examples/vignettes/output/`; +- prints the returned structures and numerical diagnostics needed to understand + the example; and +- labels genuine Python API gaps in the section where the R feature appears. + +The complete section and gap audit is in [`PARITY.md`](PARITY.md). The +machine-readable acceptance contract is in [`manifest.yml`](manifest.yml). -Run one example: +Run one full vignette: ```bash -uv run python examples/vignettes/02_partial_moments.py +uv run python examples/vignettes/07_clustering_and_regression.py ``` -Run the full example suite: +Run all nine full vignettes: ```bash uv run python examples/run_all_vignettes.py ``` + +CI uses a reduced-runtime mode without removing sections: + +```bash +NNS_VIGNETTE_FAST=1 uv run python examples/run_all_vignettes.py +``` + +The unnumbered scripts remain available as focused implementation examples and +backward-compatible entry points. New canonical material begins in R and is then +ported into the corresponding numbered Python vignette. From e014162e1e55535ab6e08a7ca913152e41fe0a09 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 16 Jul 2026 18:07:19 +0000 Subject: [PATCH 17/17] Fix and complete the canonical vignette ports Repair the failures in the full R vignette port suite: - 01_overview: replace the thin delegate shim with a full nine-section instructional port of NNSvignette_01_Overview.Rmd, including figures for the partial-moment decomposition, mode estimation, copula surfaces, nonlinear regression, ARMA forecasting, and Monte Carlo replicates. - 04_normalization: pass a list instead of a generator to np.column_stack, which numpy rejects at runtime. - 05_sampling: replace the assertion-style smoke script with a full section-by-section port of NNSvignette_05_Sampling.Rmd with CDF, inverse-CDF, Monte Carlo, and dependence-transfer figures. - 06_comparing_distributions: replace the smoke script with a full port of NNSvignette_06_Comparing_Distributions.Rmd, reconstructing the R ANOVA plot side effects from the returned statistics and adding superiority, dominance, and SD-cluster figures. - 07_clustering_and_regression: restore the file body that was committed truncated after its module docstring, porting every R section including partition orders, X-only partitioning, multivariate regression over the full grid, dimension reduction, stacking, smoothing, and imputation. - 08_classification and _vignette_support: stop rounding non-numeric arrays in the compact result renderer, and display rhs.partitions now that the multivariate result exposes it (correcting the stale gap note). - PARITY.md: update the 07/08 ledger entries to match the ports. - Lint compliance for ruff check across the examples tree, with a scoped E402 ignore for the sys.path bootstrap in instructional scripts, and ignore the generated figure output directory. All nine canonical scripts pass tests/docs/test_vignette_examples.py in fast mode and run clean in full mode; ruff and mypy pass. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_014q6V4fJ7zz6KLAgu3ddhDp --- .gitignore | 3 + examples/_vignette_support.py | 23 +- examples/vignettes/01_overview.py | 484 +++++++++++++++++- examples/vignettes/02_partial_moments.py | 9 +- .../03_correlation_and_dependence.py | 18 +- .../04_normalization_and_rescaling.py | 9 +- .../vignettes/05_sampling_and_simulation.py | 329 +++++++++++- .../vignettes/06_comparing_distributions.py | 336 ++++++++++-- .../vignettes/07_clustering_and_regression.py | 442 +++++++++++++++- examples/vignettes/08_classification.py | 18 +- examples/vignettes/09_forecasting.py | 5 +- examples/vignettes/PARITY.md | 16 +- pyproject.toml | 3 + 13 files changed, 1596 insertions(+), 99 deletions(-) diff --git a/.gitignore b/.gitignore index 10bf304..c83cce6 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,6 @@ Rplots.pdf # MkDocs build output /site/ + +# Generated vignette figures +/examples/vignettes/output/ diff --git a/examples/_vignette_support.py b/examples/_vignette_support.py index 0f532ca..5139e0a 100644 --- a/examples/_vignette_support.py +++ b/examples/_vignette_support.py @@ -9,9 +9,10 @@ import csv import os +from collections.abc import Iterable, Sequence from pathlib import Path from pprint import pformat -from typing import Any, Iterable, Sequence +from typing import Any import matplotlib @@ -68,12 +69,17 @@ def show(label: str, value: Any, *, precision: int = 6) -> None: def _compact(value: Any) -> Any: if isinstance(value, np.ndarray): + numeric = np.issubdtype(value.dtype, np.number) if value.size <= 20: - return np.round(value, 6).tolist() + return np.round(value, 6).tolist() if numeric else value.tolist() + flat = value.reshape(-1) + head, tail = flat[:5], flat[-5:] + if numeric: + head, tail = np.round(head, 6), np.round(tail, 6) return { "shape": list(value.shape), - "head": np.round(value.reshape(-1)[:5], 6).tolist(), - "tail": np.round(value.reshape(-1)[-5:], 6).tolist(), + "head": head.tolist(), + "tail": tail.tolist(), } if isinstance(value, dict): return {key: _compact(item) for key, item in value.items()} @@ -86,12 +92,17 @@ def table(headers: Sequence[str], rows: Iterable[Sequence[Any]], *, precision: i rendered = [] for row in rows: rendered.append( - [f"{item:.{precision}f}" if isinstance(item, (float, np.floating)) else str(item) for item in row] + [ + f"{item:.{precision}f}" if isinstance(item, (float, np.floating)) else str(item) + for item in row + ] ) widths = [len(str(header)) for header in headers] for row in rendered: widths = [max(width, len(cell)) for width, cell in zip(widths, row, strict=True)] - print(" ".join(str(header).ljust(width) for header, width in zip(headers, widths, strict=True))) + print( + " ".join(str(header).ljust(width) for header, width in zip(headers, widths, strict=True)) + ) print(" ".join("-" * width for width in widths)) for row in rendered: print(" ".join(cell.ljust(width) for cell, width in zip(row, widths, strict=True))) diff --git a/examples/vignettes/01_overview.py b/examples/vignettes/01_overview.py index 55de7a7..7fe8217 100644 --- a/examples/vignettes/01_overview.py +++ b/examples/vignettes/01_overview.py @@ -1,12 +1,486 @@ -"""Canonical vignette 01: Overview. +"""01. Getting Started with NNS: Overview. -Source of truth: +Instructional Python port of: 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. +A complete hands-on curriculum for Nonlinear Nonparametric Statistics using +partial moments, preserving the R vignette's nine-section sequence. Figures are +written to ``examples/vignettes/output/01_overview``. """ -from overview import main +from __future__ import annotations + +import itertools +import sys +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) + +from examples._vignette_support import ( + fast_mode, + gap, + load_iris, + load_mtcars, + note, + output_dir, + save_figure, + section, + show, + subsection, + table, +) +from nns import ( + co_lpm, + fsd_uni, + lpm, + lpm_ratio, + lpm_var, + nns_anova, + nns_arma_optim, + nns_boost, + nns_causation, + nns_copula, + nns_dep, + nns_mc, + nns_meboot, + nns_mode, + nns_moments, + nns_norm, + nns_reg, + nns_rescale, + nns_sd_cluster, + nns_seas, + nns_ss, + nns_stack, + pm_matrix, + sd_efficient_set, + ssd_uni, + tsd_uni, + upm, +) + +OUT = output_dir(__file__) + + +def main() -> None: + section("Orientation") + print( + "Goal: a complete curriculum for Nonlinear Nonparametric Statistics (NNS) " + "using partial moments. LPM(n, t, X) integrates (t - x)^n below the target " + "t; UPM(n, t, X) integrates (x - t)^n above it. The empirical estimators " + "replace the population CDF with the empirical one." + ) + + # %% [markdown] + # # 1. Foundations — Partial Moments & Variance Decomposition + # + # Classical variance treats upside and downside symmetrically. Partial + # moments separate them around a global target, so at t = mean(x): + # Var(X) = UPM(2, mu, X) + LPM(2, mu, X) (with Bessel adjustment). + section("1. Foundations — Partial Moments & Variance Decomposition") + rng = np.random.default_rng(42) + y = rng.normal(size=3000) + mu = float(np.mean(y)) + n = y.size + lower2 = float(lpm(2, mu, y)) + upper2 = float(upm(2, mu, y)) + decomposed = (lower2 + upper2) * (n / (n - 1)) + table( + ["quantity", "value"], + [ + ("LPM(2, mean, y)", lower2), + ("UPM(2, mean, y)", upper2), + ("(LPM2 + UPM2) * n/(n-1)", decomposed), + ("var(y) sample", float(np.var(y, ddof=1))), + ], + ) + + print("\nEmpirical CDF via LPM.ratio(0, t, y):") + table( + ["target", "LPM.ratio", "empirical"], + [(t, float(lpm_ratio(0, t, y)), float(np.mean(y <= t))) for t in (-1.0, 0.0, 1.0)], + ) + + z = rng.exponential(size=3000) - 1.0 + mu_z = float(np.mean(z)) + print("\nAsymmetry on a skewed distribution (expect imbalance):") + table( + ["quantity", "value"], + [("LPM(2, mean, z)", float(lpm(2, mu_z, z))), ("UPM(2, mean, z)", float(upm(2, mu_z, z)))], + ) + + fig, axes = plt.subplots(1, 2, figsize=(12, 4.6)) + for ax, (values, target, name) in zip( + axes, + ( + (y, mu, "N(0,1): symmetric partial moments"), + (z, mu_z, "exp(1) - 1: asymmetric partial moments"), + ), + strict=True, + ): + _, edges, patches = ax.hist(values, bins=40, alpha=0.8) + for edge, patch in zip(edges[:-1], patches, strict=False): + patch.set_facecolor("steelblue" if edge < target else "darkorange") + ax.axvline(target, color="black", linewidth=1.6, label="target = mean") + ax.set_title(name) + ax.legend() + save_figure(fig, OUT, "01_partial_moment_decomposition.png") + print( + "The equality holds because deviations are measured against the global " + "mean; LPM.ratio(0, t, x) constructs an empirical CDF from partial-moment counts." + ) + + # %% [markdown] + # # 2. Descriptive & Distributional Tools + section("2. Descriptive & Distributional Tools") + subsection("Higher moments from partial moments: NNS.moments") + show("NNS.moments(y)", nns_moments(y)) + + subsection("Mode estimation: NNS.mode") + rng_mode = np.random.default_rng(23) + multimodal = np.concatenate( + (rng_mode.normal(-2.0, 0.5, 1500), rng_mode.normal(2.0, 0.5, 1500)) + ) + modes = np.atleast_1d(np.asarray(nns_mode(multimodal, multi=True), dtype=float)) + show("NNS.mode(multimodal, multi = TRUE)", modes) + fig, ax = plt.subplots(figsize=(8, 4.6)) + ax.hist(multimodal, bins=60, alpha=0.75) + for mode_value in modes: + ax.axvline(mode_value, color="red", linewidth=2) + ax.set_title("Multimodal sample with NNS.mode(multi = TRUE) estimates") + save_figure(fig, OUT, "02_mode_estimation.png") + + subsection("CDF tables via LPM ratios") + quantile_probs = np.arange(0.05, 0.951, 0.1) + qgrid = np.asarray([lpm_var(p, 0, z) for p in quantile_probs], dtype=float) + table( + ["threshold (LPM.VaR)", "CDF (LPM.ratio)"], + [(float(q), float(lpm_ratio(0, q, z))) for q in qgrid], + ) + + # %% [markdown] + # # 3. Dependence & Nonlinear Association + # + # Pearson r captures linear monotone relationships. U-shapes, saturation and + # asymmetric tails can produce near-zero r despite strong dependence. + section("3. Dependence & Nonlinear Association") + rng_dep = np.random.default_rng(1) + x_dep = rng_dep.uniform(-1.0, 1.0, 2000) + y_dep = x_dep**2 + rng_dep.normal(scale=0.05, size=2000) + pearson = float(np.corrcoef(x_dep, y_dep)[0, 1]) + dependence = float(nns_dep(x_dep, y_dep)["Dependence"]) + table(["measure", "value"], [("Pearson r", pearson), ("NNS.dep", dependence)]) + + frame = np.column_stack((x_dep, y_dep, x_dep * y_dep + rng_dep.normal(scale=0.05, size=2000))) + pm = pm_matrix(1, 1, "means", frame, True, names=["a", "b", "c"]) + show("PM.matrix(1, 1, target = 'means')", pm) + show("NNS.copula(X, continuous = TRUE)", float(nns_copula(frame, continuous=True))) + + subsection("Copula visualization") + gap( + "R renders the Co.LPM copula surfaces interactively with rgl::plot3d. The " + "same surfaces are saved here as static Matplotlib 3-D figures." + ) + rng_cop = np.random.default_rng(123) + grid_n = 40 if fast_mode() else 60 + x_cop = rng_cop.normal(size=100) + y_cop = rng_cop.normal(size=100) + grid = np.asarray(list(itertools.product(x_cop[:grid_n], y_cop[:grid_n])), dtype=float) + raw_surface = np.asarray(co_lpm(0, x_cop, y_cop, grid[:, 0], grid[:, 1]), dtype=float) + u_x = np.asarray(lpm_ratio(0, x_cop, x_cop), dtype=float) + u_y = np.asarray(lpm_ratio(0, y_cop, y_cop), dtype=float) + u_grid = np.asarray(list(itertools.product(u_x[:grid_n], u_y[:grid_n])), dtype=float) + uniform_surface = np.asarray(co_lpm(0, u_x, u_y, u_grid[:, 0], u_grid[:, 1]), dtype=float) + fig = plt.figure(figsize=(12, 5)) + ax = fig.add_subplot(1, 2, 1, projection="3d") + ax.scatter(grid[:, 0], grid[:, 1], raw_surface, c="red", s=6, alpha=0.6) + ax.set_title("Co.LPM surface on raw margins") + ax = fig.add_subplot(1, 2, 2, projection="3d") + ax.scatter(u_grid[:, 0], u_grid[:, 1], uniform_surface, c="blue", s=6, alpha=0.6) + ax.set_title("Co.LPM surface on uniform margins") + save_figure(fig, OUT, "03_copula_surfaces.png") + + # %% [markdown] + # # 4. Normalization and Rescaling + section("4. Normalization and Rescaling") + subsection("NNS.norm") + rng_norm = np.random.default_rng(42) + panel = np.column_stack( + ( + rng_norm.normal(0.0, 1.0, 100), + rng_norm.normal(0.0, 5.0, 100), + rng_norm.normal(10.0, 1.0, 100), + rng_norm.normal(10.0, 10.0, 100), + ) + ) + linear_normalized = np.asarray(nns_norm(panel, linear=True), dtype=float) + show("First rows of NNS.norm(X, linear = TRUE)", linear_normalized[:6]) + print( + "Linear mode equalizes means; nonlinear mode weights each variable by its " + "dependence with the others." + ) + + subsection("Risk-neutral rescale (pricing context)") + prices = 100.0 + np.cumsum(rng_norm.normal(0.0, 1.0, 260)) + risk_neutral = nns_rescale( + prices, a=100.0, b=0.03, method="riskneutral", time_to_maturity=1.0, type="Terminal" + ) + table( + ["quantity", "value"], + [ + ("target 100 * exp(0.03)", 100.0 * float(np.exp(0.03))), + ("mean of rescaled series", float(np.mean(risk_neutral))), + ], + ) + gap("Python spells the R argument T as time_to_maturity; the operation is identical.") + + # %% [markdown] + # # 5. Hypothesis Testing, ANOVA & Stochastic Superiority + # + # Instead of distributional assumptions, groups are compared through + # LPM-based CDFs; the output is a degree of certainty, not a p-value. + section("5. Hypothesis Testing, ANOVA & Stochastic Superiority") + subsection("Two-sample and multi-group NNS.ANOVA") + rng_test = np.random.default_rng(42) + control = rng_test.normal(0.0, 1.0, 200) + treatment = rng_test.normal(0.35, 1.2, 180) + show("NNS.ANOVA(control, treatment)", nns_anova(control, treatment, random_seed=123)) + + groups = [ + rng_test.normal(0.0, 1.1, 150), + rng_test.normal(0.2, 1.0, 150), + rng_test.normal(-0.1, 0.9, 150), + ] + show("Multi-group NNS.ANOVA(control = list(g1, g2, g3), means.only = TRUE)", + nns_anova(groups, means_only=True, random_seed=123)) + + subsection("Stochastic superiority NNS.SS") + rng_ss = np.random.default_rng(123) + x_ss = rng_ss.normal(0.0, 1.0, 1000) + y_ss = rng_ss.normal(1.0, 1.0, 1000) + show("NNS.SS(x, y)", nns_ss(x_ss, y_ss)) + print( + "y has the higher mean, so P* is below 0.5: a draw from x is less likely " + "to exceed a draw from y. P* = P(X > Y) + 0.5 * P(X = Y)." + ) + with_interval = nns_ss( + x_ss, + y_ss, + confidence_interval=True, + reps=99 if fast_mode() else 999, + ci=0.95, + random_seed=123, + ) + show( + "NNS.SS(x, y, confidence.interval = TRUE)", + { + key: with_interval[key] + for key in ("p_gt", "p_tie", "p_star", "lower", "upper") + if key in with_interval + }, + ) + x_discrete = rng_ss.integers(1, 6, 100).astype(float) + y_discrete = rng_ss.integers(1, 6, 100).astype(float) + show("NNS.SS on discrete samples (explicit tie adjustment)", nns_ss(x_discrete, y_discrete)) + + # %% [markdown] + # # 6. Regression, Boosting, Stacking & Causality + # + # NNS.reg learns partitioned relationships using partial-moment weights — + # linear where appropriate, nonlinear where needed. + section("6. Regression, Boosting, Stacking & Causality") + subsection("Nonlinear regression") + rng_reg = np.random.default_rng(123) + x_train = rng_reg.uniform(-2.0, 2.0, 1000) + y_train = np.sin(np.pi * x_train) + rng_reg.normal(scale=0.2, size=1000) + x_test = np.linspace(-2.0, 2.0, 100) + regression = nns_reg(x_train, y_train, point_est=x_test) + show("NNS.reg returned structure", regression) + + fig, ax = plt.subplots(figsize=(8.5, 5)) + ax.scatter(x_train, y_train, s=8, alpha=0.35, label="training data") + ax.plot( + x_test, + np.asarray(regression["Point.est"], dtype=float), + color="red", + linewidth=2, + label="NNS point estimates", + ) + ax.plot( + x_test, + np.sin(np.pi * x_test), + color="black", + linestyle="--", + linewidth=1.4, + label="true sin(pi x)", + ) + ax.set_title("NNS.reg on sin(pi x) with noise") + ax.legend() + save_figure(fig, OUT, "04_nonlinear_regression.png") + + subsection("Classification via boosting and stacking") + iris_x, iris_y, levels = load_iris() + test_idx = np.arange(140, 150) + train_idx = np.arange(0, 140) + boost = nns_boost( + iris_x[train_idx], + iris_y[train_idx], + iris_x[test_idx], + type="CLASS", + epochs=3 if fast_mode() else 10, + learner_trials=4 if fast_mode() else 10, + status=False, + balance=True, + seed=123, + class_levels=levels, + ) + boost_accuracy = float(np.mean(np.asarray(boost["results"], dtype=float) == iris_y[test_idx])) + show("NNS.boost feature weights", boost["feature.weights"]) + show("NNS.boost feature frequency", boost["feature.frequency"]) + + stacked = nns_stack( + iris_x[train_idx], + iris_y[train_idx], + iris_x[test_idx], + type="CLASS", + balance=True, + folds=1 if fast_mode() else 2, + status=False, + seed=123, + class_levels=levels, + ) + stack_accuracy = float(np.mean(np.asarray(stacked["stack"], dtype=float) == iris_y[test_idx])) + table( + ["ensemble", "holdout accuracy (iris rows 141-150)"], + [("NNS.boost", boost_accuracy), ("NNS.stack", stack_accuracy)], + ) + + subsection("Directional causality") + mtcars = load_mtcars() + hp_to_mpg = nns_causation(mtcars["hp"], mtcars["mpg"]) + mpg_to_hp = nns_causation(mtcars["mpg"], mtcars["hp"]) + show("NNS.caus(mtcars$hp, mtcars$mpg)", hp_to_mpg) + show("NNS.caus(mtcars$mpg, mtcars$hp)", mpg_to_hp) + print("Asymmetry in the conditional-dependence scores indicates direction.") + + # %% [markdown] + # # 7. Time Series & Forecasting + # + # NNS seasonality uses the coefficient of variation instead of ACF/PACF, and + # NNS.ARMA blends multiple seasonal periods into linear or nonlinear + # regression forecasts. + section("7. Time Series & Forecasting") + rng_ts = np.random.default_rng(42) + length = 240 if fast_mode() else 480 + horizon = 24 if fast_mode() else 48 + raw_series = np.sin(np.arange(1, length + 1) / 8.0) + rng_ts.normal(scale=0.35, size=length) + series = (raw_series - raw_series.mean()) / raw_series.std(ddof=1) + + seasonal = nns_seas(series, plot=False) + show( + "NNS.seas head of all.periods", + np.asarray(seasonal["all.periods"]["Period"], dtype=float)[:6], + ) + periods = np.asarray(seasonal["periods"], dtype=int).ravel() + train_n = series.size - horizon + usable = [int(p) for p in periods if p < train_n / 4][: 2 if fast_mode() else 3] + print("Seasonal periods validated by NNS.ARMA.optim:", usable) + + optimized = nns_arma_optim( + series, + h=horizon, + seasonal_factor=usable, + print_trace=False, + plot=False, + ) + show("NNS.ARMA.optim returned structure", optimized) + forecast = np.asarray(optimized["results"], dtype=float) + fig, ax = plt.subplots(figsize=(10, 5)) + ax.plot(np.arange(series.size), series, linewidth=1.2, label="observed (scaled)") + horizon_index = np.arange(series.size, series.size + forecast.size) + ax.plot(horizon_index, forecast, color="red", linewidth=2, label="NNS.ARMA.optim forecast") + lower = optimized.get("lower.pred.int") + upper = optimized.get("upper.pred.int") + if lower is not None and upper is not None: + ax.fill_between( + horizon_index, + np.asarray(lower, dtype=float), + np.asarray(upper, dtype=float), + alpha=0.25, + label="prediction interval", + ) + ax.set_title(f"NNS.ARMA.optim forecast, h = {horizon}") + ax.legend() + save_figure(fig, OUT, "05_arma_optim_forecast.png") + note( + "R passes every detected period straight to NNS.ARMA.optim; the Python " + "optimizer validates seasonal factors against the training span, so the " + "detected periods are filtered to the supported range first." + ) + + # %% [markdown] + # # 8. Simulation & Bootstrap & Risk-Neutral Rescaling + section("8. Simulation & Bootstrap") + subsection("Maximum entropy bootstrap NNS.meboot") + x_ts = np.cumsum(rng_ts.normal(scale=0.7, size=350)) + meboot = nns_meboot(x_ts, reps=5, rho=1.0, random_seed=123) + replicates = np.asarray(meboot["replicates"], dtype=float) + show("dim(NNS.meboot replicates)", np.asarray(replicates.shape)) + + subsection("Monte Carlo over the full correlation space NNS.MC") + mc = nns_mc(x_ts, reps=5, lower_rho=-1.0, upper_rho=1.0, by=0.5, exp=1.0, random_seed=123) + print("length(mc$ensemble):", np.asarray(mc["ensemble"]).size) + print("names(mc$replicates):", list(mc["replicates"])) + show("head(mc$replicates$'rho = 0')", np.asarray(mc["replicates"]["rho = 0"], dtype=float)[:6]) + + fig, ax = plt.subplots(figsize=(9.5, 5)) + colors = plt.get_cmap("rainbow")(np.linspace(0.0, 1.0, len(mc["replicates"]))) + for (name, reps), color in zip(mc["replicates"].items(), colors, strict=True): + ax.plot(np.asarray(reps, dtype=float)[:, 0], color=color, linewidth=1.0, label=name) + ax.plot(x_ts, color="black", linewidth=2.5, label="original") + ax.set_title("NNS.MC replicates across the correlation space") + ax.legend(fontsize=8) + save_figure(fig, OUT, "06_monte_carlo_replicates.png") + + # %% [markdown] + # # 9. Portfolio & Stochastic Dominance + # + # Stochastic dominance orders uncertain prospects for broad classes of + # risk-averse utilities; partial moments supply nonparametric estimators. + section("9. Portfolio & Stochastic Dominance") + rng_sd = np.random.default_rng(42) + ra = rng_sd.normal(0.005, 0.03, 240) + rb = rng_sd.normal(0.003, 0.02, 240) + rc = rng_sd.normal(0.006, 0.04, 240) + table( + ["test", "result (1 = dominance)"], + [ + ("NNS.FSD.uni(RA, RB)", fsd_uni(ra, rb)), + ("NNS.SSD.uni(RA, RB)", ssd_uni(ra, rb)), + ("NNS.TSD.uni(RA, RB)", tsd_uni(ra, rb)), + ], + ) + returns_matrix = np.column_stack((ra, rb, rc)) + names = ["A", "B", "C"] + clusters = nns_sd_cluster(returns_matrix, degree=1, names=names) + show("NNS.SD.cluster(Rmat, degree = 1)", clusters["Clusters"]) + efficient = sd_efficient_set(returns_matrix, degree=1) + show( + "NNS.SD.efficient.set(Rmat, degree = 1)", + [names[i - 1] for i in np.asarray(efficient, dtype=int)], + ) + + section("Appendix") + print( + "The measure-theoretic sketch and the grouped quick reference in the R " + "vignette are narrative; each linked topic is executable in the numbered " + "vignettes 02 through 09 of this catalog." + ) + if __name__ == "__main__": main() diff --git a/examples/vignettes/02_partial_moments.py b/examples/vignettes/02_partial_moments.py index c8aec00..11a47b1 100644 --- a/examples/vignettes/02_partial_moments.py +++ b/examples/vignettes/02_partial_moments.py @@ -19,7 +19,14 @@ if str(ROOT) not in sys.path: sys.path.insert(0, str(ROOT)) -from examples._vignette_support import empirical_cdf, output_dir, save_figure, section, show, subsection +from examples._vignette_support import ( + empirical_cdf, + output_dir, + save_figure, + section, + show, + subsection, +) from nns import ( co_lpm, co_upm, diff --git a/examples/vignettes/03_correlation_and_dependence.py b/examples/vignettes/03_correlation_and_dependence.py index cb276bc..2700ec9 100644 --- a/examples/vignettes/03_correlation_and_dependence.py +++ b/examples/vignettes/03_correlation_and_dependence.py @@ -16,7 +16,15 @@ if str(ROOT) not in sys.path: sys.path.insert(0, str(ROOT)) -from examples._vignette_support import gap, output_dir, partition_scatter, save_figure, section, show, subsection +from examples._vignette_support import ( + gap, + output_dir, + partition_scatter, + save_figure, + section, + show, + subsection, +) from nns import nns_copula, nns_dep, nns_part OUT = output_dir(__file__) @@ -140,7 +148,13 @@ def main() -> None: pairs = [(0, 1), (0, 2), (1, 2)] for ax, (left, right) in zip(axes, pairs, strict=True): ax.scatter(independent[:, left], independent[:, right], s=8, alpha=0.35, label="sample") - ax.scatter(reference[:, left], reference[:, right], s=8, alpha=0.15, label="independence overlay") + ax.scatter( + reference[:, left], + reference[:, right], + s=8, + alpha=0.15, + label="independence overlay", + ) ax.set_title(f"X{left + 1} vs X{right + 1}") axes[0].legend() fig.suptitle(f"Multivariate NNS copula dependence = {dependence:.4f}") diff --git a/examples/vignettes/04_normalization_and_rescaling.py b/examples/vignettes/04_normalization_and_rescaling.py index 9d405c2..663ac90 100644 --- a/examples/vignettes/04_normalization_and_rescaling.py +++ b/examples/vignettes/04_normalization_and_rescaling.py @@ -17,6 +17,7 @@ sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from _vignette_support import gap, output_dir, save_figure, section, show, subsection, table + from nns import nns_norm, nns_rescale OUT = output_dir(__file__) @@ -24,7 +25,11 @@ def describe_columns(values: np.ndarray) -> list[tuple[str, float, float]]: return [ - (f"column_{index + 1}", float(np.mean(values[:, index])), float(np.std(values[:, index], ddof=1))) + ( + f"column_{index + 1}", + float(np.mean(values[:, index])), + float(np.std(values[:, index], ddof=1)), + ) for index in range(values.shape[1]) ] @@ -152,7 +157,7 @@ def main() -> None: pair_lin = np.asarray(nns_norm(pair, linear=True)) pair_nonlin = np.asarray(nns_norm(pair, linear=False)) pair_minmax = np.column_stack( - ((pair[:, i] - pair[:, i].min()) / np.ptp(pair[:, i]) for i in range(pair.shape[1])) + [(pair[:, i] - pair[:, i].min()) / np.ptp(pair[:, i]) for i in range(pair.shape[1])] ) fig, axes = plt.subplots(2, 2, figsize=(11, 8)) diff --git a/examples/vignettes/05_sampling_and_simulation.py b/examples/vignettes/05_sampling_and_simulation.py index 019ade5..dbf02bf 100644 --- a/examples/vignettes/05_sampling_and_simulation.py +++ b/examples/vignettes/05_sampling_and_simulation.py @@ -1,46 +1,329 @@ -"""Canonical vignette 05: Sampling and simulation. +"""05. Getting Started with NNS: Sampling and Simulation. -Source of truth: +Instructional Python port of: NNS/vignettes/NNSvignette_05_Sampling.Rmd -Combines partial-moment CDF/inverse-CDF sampling with maximum-entropy bootstrap -and dependence-targeted Monte Carlo examples. +Run directly to print the important returned structures and write the figures to +``examples/vignettes/output/05_sampling_and_simulation``. """ from __future__ import annotations +import math +import sys +from pathlib import Path + +import matplotlib.pyplot as plt import numpy as np -from nns import lpm_ratio, lpm_var, nns_mc, nns_meboot +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) + +from examples._vignette_support import ( + empirical_cdf, + fast_mode, + note, + output_dir, + save_figure, + section, + show, + subsection, + table, +) +from nns import lpm_ratio, lpm_var, nns_copula, nns_mc, nns_meboot + +OUT = output_dir(__file__) + +DEGREES = (0.0, 0.25, 0.5, 1.0, 2.0) + + +def spearman(a: np.ndarray, b: np.ndarray) -> float: + ranks_a = np.argsort(np.argsort(a)).astype(float) + ranks_b = np.argsort(np.argsort(b)).astype(float) + return float(np.corrcoef(ranks_a, ranks_b)[0, 1]) def main() -> None: + section("Sampling") + print( + "NNS offers novel sampling methods from any distribution, as well as " + "simulation that maintains the dependence of the original variables." + ) + + # %% [markdown] + # ## CDFs + # + # ### Empirical CDF + # + # The empirical CDF is the base construct: F(x) = P(X <= x). R uses + # ``ecdf(x)`` and evaluates the resulting step function at 0 and 1. + subsection("Empirical CDF") rng = np.random.default_rng(123) x = rng.normal(size=100) + p0 = float(np.mean(x <= 0.0)) + p1 = float(np.mean(x <= 1.0)) + table(["target", "empirical CDF"], [(0.0, p0), (1.0, p1)]) + + # %% [markdown] + # ### Lower Partial Moment CDF (LPM.ratio) + # + # With degree 0, LPM.ratio is identical to the empirical CDF: + # LPM(0, t, X) counts the share of observations at or below t. + subsection("LPM.ratio degree 0 equals the empirical CDF") + table( + ["target", "ecdf", "LPM.ratio(0, target, x)"], + [ + (0.0, p0, float(lpm_ratio(0, 0.0, x))), + (1.0, p1, float(lpm_ratio(0, 1.0, x))), + ], + ) targets = np.sort(x) + lpm_cdf = np.asarray(lpm_ratio(0, targets, x), dtype=float) 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) + max_gap = float(np.max(np.abs(lpm_cdf - empirical))) + print("Maximum |LPM.ratio - ecdf| over every observation:", max_gap) + + steps, probs = empirical_cdf(x) + fig, ax = plt.subplots(figsize=(7, 6)) + ax.step(steps, probs, where="post", color="black", label="ecdf") + ax.scatter(targets, lpm_cdf, color="red", s=18, zorder=3, label="LPM.ratio") + ax.set_xlabel("x") + ax.set_ylabel("F(x)") + ax.set_title("Empirical CDF and LPM.ratio(degree = 0)") + ax.legend(loc="upper left") + save_figure(fig, OUT, "01_ecdf_vs_lpm_ratio.png") + + # %% [markdown] + # ### LPM.ratio degree > 0 + # + # Increasing the degree parameter to any positive real number generates + # different CDFs of the initial distribution. + subsection("LPM.ratio for degrees > 0") + fig, ax = plt.subplots(figsize=(8, 6)) + ax.step(steps, probs, where="post", color="black", linewidth=3, label="ecdf") + colors = plt.get_cmap("rainbow")(np.linspace(0.0, 1.0, len(DEGREES))) + for degree, color in zip(DEGREES, colors, strict=True): + curve = np.asarray(lpm_ratio(degree, targets, x), dtype=float) + ax.plot(targets, curve, color=color, linewidth=2, label=f"LPM.ratio(degree = {degree:g})") + normal_grid = np.linspace(-3.0, 3.0, 200) + normal_cdf = 0.5 * (1.0 + np.asarray([math.erf(v / math.sqrt(2.0)) for v in normal_grid])) + ax.plot( + normal_grid, + normal_cdf, + color="black", + linestyle=":", + linewidth=2, + label="N(0,1) approximation", + ) + ax.set_xlabel("x") + ax.set_ylabel("F(x)") + ax.set_title("CDF shapes from the LPM degree") + ax.legend(loc="upper left", fontsize=8) + save_figure(fig, OUT, "02_lpm_ratio_degrees.png") - percentiles = np.linspace(0.01, 0.99, 99) + # %% [markdown] + # ### Generating PDFs with LPM.VaR + # + # LPM.VaR provides inverse-CDF estimates: sampling percentile grids under a + # degree produces new samples of the underlying distribution. + subsection("Inverse CDF sampling via LPM.VaR") + percentiles = np.linspace(0.0, 1.0, 100) 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) + for degree in DEGREES + } + + fig, axes = plt.subplots(2, 3, figsize=(14, 8)) + axes = axes.ravel() + axes[0].step(steps, probs, where="post", color="black", linewidth=2) + axes[0].set_title("eCDF via LPM.ratio()") + reference_counts, reference_edges = np.histogram(samples[0.0], bins=15) + for ax, (degree, values), color in zip(axes[1:], samples.items(), colors, strict=True): + ax.stairs(reference_counts, reference_edges, color="black", linewidth=2) + ax.hist(values, bins=15, color=color, alpha=0.5) + ax.set_title(f"Inverse CDF via LPM.VaR(degree {degree:g})") + ax.set_xlabel("x") + ax.set_ylabel("freq") + save_figure(fig, OUT, "03_lpm_var_inverse_cdf.png") + + print("First 10 samples from each degree compared with the original x:") + sorted_x = np.sort(x) + table( + ["original x"] + [f"degree {degree:g}" for degree in DEGREES], + [ + tuple([float(sorted_x[i])] + [float(samples[degree][i]) for degree in DEGREES]) + for i in range(10) + ], + ) + print( + "Degree 0 reproduces the empirical quantiles exactly; higher degrees " + "concentrate samples toward the distribution center." + ) + + section("Simulation") + + # %% [markdown] + # ## Bootstrapping (NNS.meboot) + # + # NNS.meboot is based on the maximum entropy bootstrap, designed for time + # series and avoiding the IID assumption. Sampling from specified + # correlations ensures the full spectrum of future paths is sampled; + # typical Monte Carlo samples are restricted to [-0.3, 0.3] correlation + # with the original data. NNS.MC is a streamlined wrapper generating one + # replicate for each rho in a sequence. + subsection("NNS.MC replicates across the correlation space") + mc = nns_mc(x, reps=1, lower_rho=-1.0, upper_rho=1.0, by=0.5, random_seed=123) + replicates = { + name: np.asarray(mc["replicates"][name], dtype=float).ravel() for name in mc["replicates"] } - 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"])) + show("NNS.MC returned structure", mc) + + fig, ax = plt.subplots(figsize=(9, 6)) + rho_colors = plt.get_cmap("rainbow")(np.linspace(0.0, 1.0, len(replicates))) + for (name, path), color in zip(replicates.items(), rho_colors, strict=True): + ax.plot(path, color=color, linewidth=1.2, label=name) + ax.plot(x, color="black", linewidth=3, label="original x") + ax.set_xlabel("observation") + ax.set_title("NNS.MC replicates versus the original series") + ax.legend(fontsize=8) + save_figure(fig, OUT, "04_mc_replicates.png") + + print("Replicate Spearman correlations with the original x:") + table( + ["replicate", "Spearman correlation"], + [(name, spearman(path, x)) for name, path in replicates.items()], + ) + note( + "The rho target lives on each replicate individually, and it is measured " + "in the metric that `type` targets. The ensemble is the per-observation " + "mean of replicates, so correlations taken on the ensemble read higher " + "than the per-replicate rho. Calibrate rho on the replicates." + ) + + # %% [markdown] + # ### target_drift Specification + # + # A target drift can be requested for the replicates with target_drift. + subsection("target_drift specification") + drifted = nns_mc( + x, + reps=1, + lower_rho=-1.0, + upper_rho=1.0, + by=0.5, + target_drift=0.05, + random_seed=123, + ) + drift_replicates = { + name: np.asarray(drifted["replicates"][name], dtype=float).ravel() + for name in drifted["replicates"] + } + fig, ax = plt.subplots(figsize=(9, 6)) + ax.plot(x, color="black", linewidth=3, label="original x") + for (name, path), color in zip(drift_replicates.items(), rho_colors, strict=True): + ax.plot(path, color=color, linewidth=1.2, label=name) + ax.set_xlabel("observation") + ax.set_title("NNS.MC replicates with target_drift = 0.05") + ax.legend(fontsize=8) + save_figure(fig, OUT, "05_mc_replicates_target_drift.png") + + # %% [markdown] + # ## Simulating a Multivariate Dependence Structure + # + # Analogous to an empirical copula transformation: (1) determine the + # dependence structure with LPM.ratio(1, x, x); (2) generate or supply new + # data of any distribution or length; (3) map the new data through LPM.VaR + # at the stored positions. + subsection("Multivariate dependence transfer") + rng = np.random.default_rng(123) + n = 250 if fast_mode() else 1000 + x_sim = rng.normal(size=n) + y_sim = rng.normal(size=n) + z_sim = rng.normal(size=n) + original = np.column_stack((x_sim, y_sim, z_sim, x_sim)) + print("x is duplicated in the panel to avoid total independence (example only).") + + dependence_structure = np.column_stack( + [ + np.asarray(lpm_ratio(1, original[:, j], original[:, j]), dtype=float) + for j in range(original.shape[1]) + ] + ) + new_data = rng.normal(10.0, 20.0, size=(original.shape[0] * 2, original.shape[1])) + new_dep_data = np.column_stack( + [ + np.asarray( + [lpm_var(p, 1, new_data[:, j]) for p in dependence_structure[:, j]], + dtype=float, + ) + for j in range(original.shape[1]) + ] + ) + + show("head(original.data)", original[:6]) + show("head(new.dep.data)", new_dep_data[:6]) + table( + ["panel", "NNS.copula"], + [ + ("original.data", float(nns_copula(original))), + ("new.dep.data", float(nns_copula(new_dep_data))), + ], + ) + print( + "Similar multivariate dependence with radically different values, since " + "N(10, 20) draws replaced the original N(0, 1) observations." + ) + + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + axes[0].scatter(original[:, 0], original[:, 1], s=10, alpha=0.5) + axes[0].set_title("original.data: x vs y") + axes[1].scatter(new_dep_data[:, 0], new_dep_data[:, 1], s=10, alpha=0.5, color="darkorange") + axes[1].set_title("new.dep.data: same dependence, new scale") + for ax in axes: + ax.set_xlabel("first variable") + ax.set_ylabel("second variable") + save_figure(fig, OUT, "06_dependence_transfer.png") + + # %% [markdown] + # ## Alternative Using NNS.meboot + # + # To keep simulated values close to the original data, apply NNS.meboot to + # each variable with rho = 0.95 and use the replicate ensembles. + subsection("NNS.meboot alternative") + reps = 10 if fast_mode() else 100 + boots = [ + nns_meboot(original[:, j], reps=reps, rho=0.95, random_seed=123) + for j in range(original.shape[1]) + ] + boot_matrix = np.column_stack([np.asarray(b["ensemble"], dtype=float).ravel() for b in boots]) + + print("Ensemble Spearman correlations with original.data:") + table( + ["variable", "Spearman correlation"], + [ + (name, spearman(boot_matrix[:, j], original[:, j])) + for j, name in enumerate(["x", "y", "z", "x (duplicate)"]) + ], + ) + show("head(new.boot.dep.matrix)", boot_matrix[:6]) + table( + ["panel", "NNS.copula"], + [ + ("original.data", float(nns_copula(original))), + ("new.boot.dep.matrix", float(nns_copula(boot_matrix))), + ], + ) + print("Similar dependence with similar values.") + note( + "The ensemble is the mean of the replicates, so its correlation with the " + "original series reads above the per-replicate rho of 0.95; the rho " + "alignment is calibrated on the individual replicates." + ) + note( + "Stochastic draws use NumPy generators, so simulations are " + "distributionally equivalent to R but not bit-for-bit identical draws." + ) if __name__ == "__main__": diff --git a/examples/vignettes/06_comparing_distributions.py b/examples/vignettes/06_comparing_distributions.py index 0d864a8..06bf8d8 100644 --- a/examples/vignettes/06_comparing_distributions.py +++ b/examples/vignettes/06_comparing_distributions.py @@ -1,61 +1,321 @@ -"""Canonical vignette 06: Comparing distributions. +"""06. Getting Started with NNS: Comparing Distributions. -Source of truth: +Instructional Python port of: NNS/vignettes/NNSvignette_06_Comparing_Distributions.Rmd -Covers NNS ANOVA certainty, stochastic superiority, pairwise stochastic -dominance, the stochastic-dominance efficient set, and SD clustering. +Run directly to print the important returned structures and write the figures to +``examples/vignettes/output/06_comparing_distributions``. """ from __future__ import annotations +import sys +from pathlib import Path + +import matplotlib.pyplot as plt import numpy as np +from scipy import stats + +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) -from nns import ( - fsd_uni, - nns_anova, - nns_sd_cluster, - nns_ss, - sd_efficient_set, - ssd_uni, - tsd_uni, +from examples._vignette_support import ( + empirical_cdf, + fast_mode, + gap, + load_mtcars, + note, + output_dir, + save_figure, + section, + show, + subsection, + table, ) +from nns import fsd, nns_anova, nns_sd_cluster, nns_ss, sd_efficient_set + +OUT = output_dir(__file__) + + +def anova_figure( + control: np.ndarray, + treatment: np.ndarray, + result: dict[str, object], + title: str, + filename: str, +) -> None: + """Reconstruct R's NNS.ANOVA plot side effect from the returned statistics.""" + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + control_stat = float(np.asarray(result["Control"], dtype=float)) + treatment_stat = float(np.asarray(result["Treatment"], dtype=float)) + grand = float(np.asarray(result["Grand_Statistic"], dtype=float)) + axes[0].bar( + ["control", "treatment", "grand"], + [control_stat, treatment_stat, grand], + color=["steelblue", "darkorange", "gray"], + ) + lower = result.get("Lower Bound Robust Certainty") + upper = result.get("Upper Bound Robust Certainty") + robust = result.get("Robust Certainty Estimate") + label = f"Certainty = {float(np.asarray(result['Certainty'], dtype=float)):.4f}" + if robust is not None: + robust_value = float(np.asarray(robust, dtype=float)) + bounds = (float(np.asarray(lower, dtype=float)), float(np.asarray(upper, dtype=float))) + label += f"\nrobust = {robust_value:.4f} [{bounds[0]:.4f}, {bounds[1]:.4f}]" + axes[0].set_title(label) + axes[0].set_ylabel("group statistic") + + for values, name, color in ( + (control, "control", "steelblue"), + (treatment, "treatment", "darkorange"), + ): + ordered, probs = empirical_cdf(values) + axes[1].step(ordered, probs, where="post", label=name, color=color) + axes[1].axvline(grand, linestyle="--", color="gray", label="grand statistic") + axes[1].set_xlabel("value") + axes[1].set_ylabel("F(x)") + axes[1].set_title(title) + axes[1].legend() + save_figure(fig, OUT, filename) def main() -> None: + section("Comparing Distributions") + print( + "NNS tests whether distributions came from the same population, or share " + "the same mean or median, through NNS.ANOVA. The output is a Certainty " + "statistic on [0, 1], with 1 representing identical distributions." + ) + gap( + "R's NNS.ANOVA(..., plot = TRUE) draws its diagnostics as a side effect. " + "The Python API returns the statistics without plotting, so each figure " + "here is reconstructed from the returned values and empirical CDFs." + ) + + # %% [markdown] + # ## Test if Same Population + # + # Do automatic and manual transmissions have significantly different mpg + # distributions per the mtcars dataset? + subsection("Test if same population: mtcars mpg by transmission") + mtcars = load_mtcars() + mpg_auto = mtcars["mpg"][mtcars["am"] == 1] + mpg_manual = mtcars["mpg"][mtcars["am"] == 0] + same_population = nns_anova( + mpg_manual, + mpg_auto, + robust=True, + n_boot=100 if fast_mode() else 1000, + random_seed=123, + ) + show( + "NNS.ANOVA(control = manual mpg, treatment = automatic mpg, robust = TRUE)", + same_population, + ) + anova_figure( + mpg_manual, + mpg_auto, + same_population, + "mtcars mpg by transmission", + "01_anova_mtcars.png", + ) + print( + "The Certainty shows these two distributions clearly do not come from the " + "same population. The nonparametric Mann-Whitney-Wilcoxon test agrees:" + ) + wilcox = stats.mannwhitneyu(mpg_manual, mpg_auto, alternative="two-sided") + table(["statistic", "p-value"], [(float(wilcox.statistic), float(wilcox.pvalue))]) + + # %% [markdown] + # ## Test if means are Equal + # + # Two Normal samples share a mean of zero with different standard + # deviations; both NNS.ANOVA and the t-test should be comfortable that the + # means are equal. + subsection("Test if means are equal") 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 + x = rng.normal(0.0, 1.0, 1000) + y = rng.normal(0.0, 2.0, 1000) + equal_means = nns_anova( + x, + y, + means_only=True, + robust=True, + n_boot=100 if fast_mode() else 1000, + random_seed=123, + ) + show("NNS.ANOVA(x, y, means.only = TRUE, robust = TRUE)", equal_means) + anova_figure(x, y, equal_means, "Equal means, unequal variance", "02_anova_equal_means.png") + ttest = stats.ttest_ind(x, y, equal_var=False) + table(["t statistic", "p-value"], [(float(ttest.statistic), float(ttest.pvalue))]) + + # %% [markdown] + # ## Test if means are Unequal + # + # Shifting y's mean to 1 shows the sensitivity of both methods, which firmly + # reject equal means. The effect size interval reported by NNS.ANOVA should + # contain the specified shift of 1. + subsection("Test if means are unequal") + y_shifted = rng.normal(1.0, 1.0, 1000) + unequal_means = nns_anova( + x, + y_shifted, + means_only=True, + robust=True, + n_boot=100 if fast_mode() else 1000, + random_seed=123, + ) + show("NNS.ANOVA(x, y, means.only = TRUE, robust = TRUE)", unequal_means) + anova_figure(x, y_shifted, unequal_means, "Unequal means", "03_anova_unequal_means.png") + ttest = stats.ttest_ind(x, y_shifted, equal_var=False) + table(["t statistic", "p-value"], [(float(ttest.statistic), float(ttest.pvalue))]) + print( + "Effect size bounds " + f"[{float(np.asarray(unequal_means['Effect_Size_LB'], dtype=float)):.4f}, " + f"{float(np.asarray(unequal_means['Effect_Size_UB'], dtype=float)):.4f}] " + "contain the specified shift of 1." + ) - 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) + # %% [markdown] + # ## Medians + # + # Setting both means.only = TRUE and medians = TRUE tests medians instead. + subsection("Medians") + unequal_medians = nns_anova( + x, + y_shifted, + means_only=True, + medians=True, + robust=True, + n_boot=100 if fast_mode() else 1000, + random_seed=123, + ) + show("NNS.ANOVA(x, y, means.only = TRUE, medians = TRUE, robust = TRUE)", unequal_medians) - 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 + # %% [markdown] + # # Stochastic Superiority + # + # Rather than testing equality, stochastic superiority measures the + # probability that a random draw from one distribution exceeds a random + # draw from another: P* = P(X > Y) + 0.5 * P(X = Y). A value of 0.5 means + # no directional advantage. + section("Stochastic Superiority") + superiority = nns_ss(x, y_shifted) + show("NNS.SS(x, y)", superiority) + print( + "y was generated with the higher mean, so P* for x versus y is below " + "0.5: a draw from x is less likely to exceed a draw from y." + ) - 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] + subsection("Bootstrap confidence interval") + with_interval = nns_ss( + x, + y_shifted, + confidence_interval=True, + reps=99 if fast_mode() else 999, + ci=0.95, + random_seed=123, + ) + interval_keys = ["p_gt", "p_tie", "p_star", "lower", "upper"] + show( + "NNS.SS(x, y, confidence.interval = TRUE)", + {key: with_interval[key] for key in interval_keys if key in with_interval}, ) - 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) + subsection("Discrete variables and ties") + x_discrete = rng.integers(1, 6, 100).astype(float) + y_discrete = rng.integers(1, 6, 100).astype(float) + discrete_superiority = nns_ss(x_discrete, y_discrete) + show("NNS.SS on discrete 1..5 samples", discrete_superiority) print( - "FSD/SSD/TSD shifted over base:", - fsd_uni(shifted, x), - ssd_uni(shifted, x), - tsd_uni(shifted, x), + "Ties occur with positive probability, so p_tie and p_star reflect the " + "adjustment explicitly." + ) + + fig, ax = plt.subplots(figsize=(8, 5.5)) + for values, name, color in ((x, "x", "steelblue"), (y_shifted, "y", "darkorange")): + ordered, probs = empirical_cdf(values) + ax.step(ordered, probs, where="post", label=name, color=color) + p_star = float(np.asarray(superiority["p_star"], dtype=float)) + ax.set_title(f"Stochastic superiority: P*(x > y) = {p_star:.4f}") + ax.set_xlabel("value") + ax.set_ylabel("F(x)") + ax.legend() + save_figure(fig, OUT, "04_stochastic_superiority.png") + + # %% [markdown] + # # Stochastic Dominance + # + # First, second, and third degree dominance tests are available via + # NNS.FSD, NNS.SSD and NNS.TSD. NNS.FSD correctly identifies the shift in + # the y variable specified in the unequal means example. + section("Stochastic Dominance") + fsd_result = fsd(y_shifted, x) + print("NNS.FSD(y, x):", fsd_result, "(1 indicates y first-degree dominates x)") + fig, ax = plt.subplots(figsize=(8, 5.5)) + for values, name, color in ((x, "x", "steelblue"), (y_shifted, "y", "darkorange")): + ordered, probs = empirical_cdf(values) + ax.step(ordered, probs, where="post", label=name, color=color) + ax.set_title("First-degree dominance: the y CDF sits at or right of x everywhere") + ax.set_xlabel("value") + ax.set_ylabel("F(x)") + ax.legend() + save_figure(fig, OUT, "05_stochastic_dominance.png") + gap( + "R's NNS.FSD plots both CDFs as a side effect and returns the degree of " + "dominance. Python returns the indicator, so the CDF comparison figure is " + "generated from the same inputs." + ) + + # %% [markdown] + # ## Stochastic Dominant Efficient Sets + # + # x2, x4, x6, x8 dominate their preceding distributions yet do not dominate + # one another, forming the first-degree efficient set. + subsection("Stochastic dominant efficient sets") + rng = np.random.default_rng(123) + base = [rng.normal(size=1000) for _ in range(4)] + columns: list[np.ndarray] = [] + for series in base: + columns.extend((series, series + 1.0)) + panel = np.column_stack(columns) + names = [f"x{i + 1}" for i in range(panel.shape[1])] + efficient = sd_efficient_set(panel, degree=1) + efficient_names = [names[i - 1] for i in np.asarray(efficient, dtype=int)] + show("NNS.SD.efficient.set(cbind(x1..x8), degree = 1)", efficient_names) + + # %% [markdown] + # ## Stochastic Dominant Clusters + # + # Clusters are assigned to non-dominated constituents; R also renders a + # dendrogram from the returned linkage. + subsection("Stochastic dominant clusters") + clusters = nns_sd_cluster(panel, degree=1, dendrogram=True, names=names) + show("NNS.SD.cluster(cbind(x1..x8), degree = 1, dendrogram = TRUE)", clusters["Clusters"]) + + fig, ax = plt.subplots(figsize=(9, 5.5)) + cluster_map = clusters["Clusters"] + palette = plt.get_cmap("tab10") + for cluster_index, (cluster_name, members) in enumerate(dict(cluster_map).items()): + for member in np.atleast_1d(np.asarray(members, dtype=str)): + column = names.index(str(member)) + ordered, probs = empirical_cdf(panel[:, column]) + ax.step( + ordered, + probs, + where="post", + color=palette(cluster_index), + label=f"{member} ({cluster_name})", + ) + ax.set_xlabel("value") + ax.set_ylabel("F(x)") + ax.set_title("First-degree stochastic dominance clusters") + ax.legend(fontsize=8) + save_figure(fig, OUT, "06_sd_clusters.png") + note( + "The returned Dendrogram entry carries the linkage used by R's plot; the " + "cluster CDFs above color each member by its assigned cluster." ) - print("SD efficient set:", efficient) - print("SD clusters:", clusters["Clusters"]) if __name__ == "__main__": diff --git a/examples/vignettes/07_clustering_and_regression.py b/examples/vignettes/07_clustering_and_regression.py index 6c3b724..6221fce 100644 --- a/examples/vignettes/07_clustering_and_regression.py +++ b/examples/vignettes/07_clustering_and_regression.py @@ -8,4 +8,444 @@ """ from __future__ import annotations -import \ No newline at end of file +import itertools +import sys +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +ROOT = Path(__file__).resolve().parents[2] +if str(ROOT) not in sys.path: + sys.path.insert(0, str(ROOT)) + +from examples._vignette_support import ( + fast_mode, + gap, + load_iris, + note, + output_dir, + partition_scatter, + regression_scatter, + save_figure, + section, + show, + subsection, + table, +) +from nns import nns_part, nns_reg, nns_stack + +OUT = output_dir(__file__) + + +def main() -> None: + section("Clustering and Regression") + print( + "NNS.part is both a partitional and hierarchical clustering method: it " + "iteratively splits the joint distribution into partial-moment quadrants and " + "assigns a quadrant identification at each partition. The quadrant means are " + "the regression points reused by NNS.reg." + ) + + # %% [markdown] + # ## NNS Partitioning NNS.part() + # + # R draws each order's Voronoi tessellation as a side effect of + # ``NNS.part(..., Voronoi = TRUE)``. Python has no Voronoi rendering flag, so + # each returned partition is drawn directly from its quadrant memberships and + # regression points. + subsection("NNS.part(): orders 1 through 4") + x = np.round(np.arange(-5.0, 5.0001, 0.05), 10) + y = x**3 + + fig, axes = plt.subplots(2, 2, figsize=(11, 9)) + for order, ax in zip(range(1, 5), axes.ravel(), strict=True): + part = nns_part(x, y, order=order, obs_req=0) + partition_scatter(ax, part, f"NNS.part order = {order}") + save_figure(fig, OUT, "01_partition_orders.png") + + final = nns_part(x, y, order=4, obs_req=0) + show("NNS.part(x, y, order = 4, obs.req = 0) returned structure", final) + gap( + "R's NNS.part(..., Voronoi = TRUE) draws the tessellation itself. The Python " + "partition result contains the same observation quadrant assignments and " + "regression points, which these figures color directly." + ) + + # %% [markdown] + # ### X-only Partitioning + # + # ``type = "XONLY"`` partitions on the x values alone, using the entire + # bandwidth for the regression point derivation. Identifications are limited + # to 1's and 2's (left and right of each partition) rather than the four + # quadrant values of joint partitioning. + subsection("X-only partitioning") + fig, axes = plt.subplots(2, 2, figsize=(11, 9)) + for order, ax in zip(range(1, 5), axes.ravel(), strict=True): + part = nns_part(x, y, order=order, type="XONLY") + partition_scatter(ax, part, f'NNS.part type = "XONLY", order = {order}') + save_figure(fig, OUT, "02_partition_orders_xonly.png") + + xonly = nns_part(x, y, order=4, type="XONLY") + show('NNS.part(x, y, order = 4, type = "XONLY") returned structure', xonly) + ids = np.unique(np.asarray(xonly["dt"]["quadrant"], dtype=str)) + print("Distinct X-only partition identifications use only 1s and 2s per level:") + print(sorted(ids)[:8], "..." if ids.size > 8 else "") + + # %% [markdown] + # ## Clusters Used in Regression + # + # The left column shows the partitions for orders 1 through 3; the right + # column shows the corresponding NNS regression built from those clusters. + subsection("Clusters used in regression") + fig, axes = plt.subplots(3, 2, figsize=(11, 13)) + for order in range(1, 4): + part = nns_part(x, y, order=order, obs_req=0, type="XONLY") + partition_scatter(axes[order - 1][0], part, f"NNS.part order = {order}") + reg = nns_reg(x, y, order=order) + regression_scatter(axes[order - 1][1], x, y, reg, f"NNS.reg order = {order}") + save_figure(fig, OUT, "03_clusters_used_in_regression.png") + + # %% [markdown] + # # NNS Regression NNS.reg() + # + # ## Univariate + subsection("Univariate NNS.reg") + uni = nns_reg(x, y) + fig, ax = plt.subplots(figsize=(8, 5)) + regression_scatter(ax, x, y, uni, "NNS.reg(x, y)") + save_figure(fig, OUT, "04_univariate_regression.png") + show("Univariate NNS.reg returned structure", uni) + + # %% [markdown] + # ## Multivariate + # + # The R vignette fits f(x, y) = x^3 + 3y - y^3 - 3x over the full + # expand.grid of x with itself using ``order = "max"``. The full 201-point + # grid (40,401 rows) is preserved outside fast mode; CI uses a coarser grid + # solely to reduce runtime. + subsection("Multivariate NNS.reg") + grid_axis = np.round(np.arange(-5.0, 5.0001, 0.25 if fast_mode() else 0.05), 10) + grid = np.asarray(list(itertools.product(grid_axis, grid_axis)), dtype=float) + g = grid[:, 0] ** 3 + 3.0 * grid[:, 1] - grid[:, 1] ** 3 - 3.0 * grid[:, 0] + multi = nns_reg(grid, g, order="max") + show("Multivariate NNS.reg returned structure", multi) + show("Per-regressor partitions (rhs.partitions)", multi["rhs.partitions"]) + show("Regression point matrix (RPM)", multi["RPM"]) + + fitted = np.asarray(multi["Fitted.xy"]["y.hat"], dtype=float) + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + axes[0].scatter(g, fitted, s=6, alpha=0.4) + axes[0].plot([g.min(), g.max()], [g.min(), g.max()], linestyle="--", linewidth=1.5) + axes[0].set_xlabel("observed g") + axes[0].set_ylabel("fitted g") + axes[0].set_title(f"Multivariate fit; R2 = {float(multi['R2']):.4f}") + surface = axes[1].tricontourf(grid[:, 0], grid[:, 1], fitted, levels=20, cmap="viridis") + fig.colorbar(surface, ax=axes[1]) + axes[1].set_xlabel("x") + axes[1].set_ylabel("y") + axes[1].set_title("Fitted surface for f(x, y) = x^3 + 3y - y^3 - 3x") + save_figure(fig, OUT, "05_multivariate_regression.png") + + # %% [markdown] + # ## Inter/Extrapolation + # + # ``point.est`` accepts any data of the regressors' dimension; estimates are + # returned in ``$Point.est``. + subsection("Inter/extrapolation with point.est") + points = np.asarray([[-5.5, -5.5], [0.25, 0.25], [5.5, 5.5]], dtype=float) + extrap = nns_reg(grid, g, order="max", point_est=points) + f_true = points[:, 0] ** 3 + 3.0 * points[:, 1] - points[:, 1] ** 3 - 3.0 * points[:, 0] + table( + ["x", "y", "true f(x, y)", "NNS point estimate"], + [ + (float(p[0]), float(p[1]), float(t), float(e)) + for p, t, e in zip( + points, f_true, np.asarray(extrap["Point.est"], dtype=float), strict=True + ) + ], + ) + note( + "The first row lies outside the training support, demonstrating " + "extrapolation; the second interpolates between grid nodes." + ) + + # %% [markdown] + # ## NNS Dimension Reduction Regression + # + # ``dim.red.method = "cor"`` reduces all regressors to a single dimension + # using the returned equation. + subsection("Dimension reduction regression") + iris_x, iris_y, levels = load_iris() + feature_names = ["Sepal.Length", "Sepal.Width", "Petal.Length", "Petal.Width"] + dim_red = nns_reg(iris_x, iris_y, dim_red_method="cor", class_levels=levels) + equation = dim_red["equation"] + coefficients = np.asarray(equation["Coefficient"], dtype=float) + show("Dimension reduction equation", equation) + terms = " + ".join( + f"{coefficients[i]:.3f}*{feature_names[i]}" for i in range(len(feature_names)) + ) + print(f"\nSpecies = ({terms}) / {coefficients[-1]:g}") + + # %% [markdown] + # ### Threshold + # + # ``threshold = 0.75`` drops regressors whose absolute correlation falls + # below the requested level, reducing the denominator accordingly. + subsection("Dimension reduction with threshold") + reduced = nns_reg(iris_x, iris_y, dim_red_method="cor", threshold=0.75, class_levels=levels) + equation = reduced["equation"] + coefficients = np.asarray(equation["Coefficient"], dtype=float) + show("Thresholded dimension reduction equation", equation) + terms = " + ".join( + f"{coefficients[i]:.3f}*{feature_names[i]}" for i in range(len(feature_names)) + ) + print(f"\nSpecies = ({terms}) / {coefficients[-1]:g}") + + point_est = nns_reg( + iris_x, + iris_y, + dim_red_method="cor", + threshold=0.75, + point_est=iris_x[:10], + class_levels=levels, + ) + show("Point estimates for iris rows 1 through 10", np.asarray(point_est["Point.est"])) + + # %% [markdown] + # # Classification + # + # ``type = "CLASS"`` rounds estimates to class values. The base category of + # the response must be 1, not 0. + subsection("Classification") + classified = nns_reg(iris_x, iris_y, type="CLASS", point_est=iris_x[:10], class_levels=levels) + show( + "Classification point estimates for iris rows 1 through 10", + np.asarray(classified["Point.est"]), + ) + note( + "The response is encoded 1=setosa, 2=versicolor, 3=virginica, preserving " + "the R base-1 contract." + ) + + # %% [markdown] + # # Cross-Validation NNS.stack() + # + # NNS.stack cross-validates n.best for the multivariate regression and the + # threshold for the dimension-reduced regression, then ensembles both. The + # objective function is the classification rate, matching the R call + # ``obj.fn = expression(mean(round(predicted) == actual))``. + subsection("NNS.stack cross-validation") + stack = nns_stack( + iris_x, + iris_y, + iris_x[:10], + type="CLASS", + dim_red_method="cor", + obj_fn=lambda predicted, actual: float(np.mean(np.round(predicted) == actual)), + objective="max", + folds=1, + status=False, + seed=123, + class_levels=levels, + ) + show("NNS.stack returned structure", stack) + + # %% [markdown] + # # Increasing Dimensions + # + # Multicollinearity is not an issue for nonparametric regression, so an + # ill-fit univariate model can duplicate its regressor, cross-validate + # n.best, and regress in the doubled space. + subsection("Increasing dimensions") + rng = np.random.default_rng(123) + x_rand = rng.normal(size=100) + y_rand = rng.normal(size=100) + doubled = np.column_stack((x_rand, x_rand)) + params = nns_stack(doubled, y_rand, method=1, folds=1, status=False, seed=123) + n_best = params["NNS.reg.n.best"] + show("Cross-validated parameters from NNS.stack(method = 1)", params) + + doubled_fit = nns_reg( + doubled, + y_rand, + n_best=n_best, + point_est=doubled, + confidence_interval=0.95, + ) + estimates = np.asarray(doubled_fit["Point.est"], dtype=float) + residuals = y_rand - estimates + pred_int = doubled_fit["pred.int"] + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + order_idx = np.argsort(x_rand) + axes[0].scatter(x_rand, y_rand, s=16, alpha=0.6, label="observed") + axes[0].plot( + x_rand[order_idx], estimates[order_idx], linewidth=1.4, label="NNS fit (cbind(x, x))" + ) + if isinstance(pred_int, dict) and pred_int: + keys = list(pred_int) + lower = np.asarray(pred_int[keys[0]], dtype=float) + upper = np.asarray(pred_int[keys[-1]], dtype=float) + axes[0].fill_between( + x_rand[order_idx], + lower[order_idx], + upper[order_idx], + alpha=0.2, + label="95% interval", + ) + axes[0].set_title(f"Duplicated regressors, n.best = {n_best:g}") + axes[0].legend() + axes[1].scatter(estimates, residuals, s=16, alpha=0.6) + axes[1].axhline(0.0, linewidth=1.0, linestyle="--") + axes[1].set_xlabel("fitted") + axes[1].set_ylabel("residual") + axes[1].set_title("Residual diagnostic (R residual.plot equivalent)") + save_figure(fig, OUT, "06_increasing_dimensions.png") + + # %% [markdown] + # # Smoothing Option + # + # Smoothness is not required, but ``smooth = TRUE`` applies a smoothing + # spline to the internally generated regression points. + subsection("Smoothing option") + rough = nns_reg(x_rand, y_rand) + smooth = nns_reg(x_rand, y_rand, smooth=True) + fig, axes = plt.subplots(1, 2, figsize=(12, 5)) + regression_scatter(axes[0], x_rand, y_rand, rough, "NNS.reg(x, y)") + regression_scatter(axes[1], x_rand, y_rand, smooth, "NNS.reg(x, y, smooth = TRUE)") + save_figure(fig, OUT, "07_smoothing_option.png") + show("Smoothed regression returned structure", smooth) + + # %% [markdown] + # # Imputation + # + # Imputation is a direct application of nearest-neighbor regression: the + # observed (X, y) pairs are the training set and the predictors of the + # missing rows are point.est. With ``order = "max", n.best = 1`` each + # missing y is filled by its closest donor, so imputations remain strictly + # within the support of the observed data. + + # %% [markdown] + # ## Univariate Imputation + # + # The increasing dimensions trick duplicates x into cbind(x, x) so the + # distance function operates in a 2-D space, sharpening donor selection. + subsection("Univariate imputation") + rng = np.random.default_rng(123) + n = 200 if fast_mode() else 400 + x_uni = np.sort(rng.uniform(-3.0, 3.0, n)) + y_uni = np.sin(x_uni) + 0.2 * x_uni**2 + rng.normal(0.0, 0.25, n) + missing = rng.binomial(1, 0.25, n) == 1 + + x2_train = np.column_stack((x_uni[~missing], x_uni[~missing])) + x2_missing = np.column_stack((x_uni[missing], x_uni[missing])) + imputed_uni = np.asarray( + nns_reg(x2_train, y_uni[~missing], point_est=x2_missing, order="max", n_best=1)[ + "Point.est" + ], + dtype=float, + ) + completed = y_uni.copy() + completed[missing] = imputed_uni + table( + ["quantity", "value"], + [ + ("missing rows", int(np.sum(missing))), + ( + "imputation RMSE vs truth", + float(np.sqrt(np.mean((imputed_uni - y_uni[missing]) ** 2))), + ), + ], + ) + fig, ax = plt.subplots(figsize=(8, 6)) + ax.scatter( + x_uni, + y_uni, + s=26, + facecolors="none", + edgecolors="steelblue", + linewidths=1.4, + label="Observed", + ) + ax.scatter( + x_uni[missing], imputed_uni, s=30, marker="s", color="red", label="Imputed (NNS 1-NN)" + ) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_title("NNS 1-NN Imputation") + ax.legend(loc="upper left", frameon=False) + save_figure(fig, OUT, "08_univariate_imputation.png") + + # %% [markdown] + # ## Multivariate Imputation + # + # The same form applies directly with the full predictor set: observed rows + # train the model and incomplete rows are point.est. + subsection("Multivariate imputation") + rng = np.random.default_rng(123) + n = 300 if fast_mode() else 800 + predictors = np.column_stack( + (rng.normal(size=n), rng.uniform(-2.0, 2.0, n), rng.normal(0.0, 1.0, n)) + ) + + def truth(p: np.ndarray) -> np.ndarray: + return ( + 1.1 * p[:, 0] + - 0.8 * p[:, 1] + + 0.5 * p[:, 2] + + 0.6 * p[:, 0] * p[:, 1] + - 0.4 * p[:, 1] * p[:, 2] + + 0.3 * np.sin(1.3 * p[:, 0]) + ) + + y_multi = truth(predictors) + rng.normal(0.0, 0.4, n) + missing = rng.binomial(1, 0.30, n) == 1 + imputed_multi = np.asarray( + nns_reg( + predictors[~missing], + y_multi[~missing], + point_est=predictors[missing], + order="max", + n_best=1, + )["Point.est"], + dtype=float, + ) + table( + ["quantity", "value"], + [ + ("missing rows", int(np.sum(missing))), + ( + "imputation RMSE vs truth", + float(np.sqrt(np.mean((imputed_multi - y_multi[missing]) ** 2))), + ), + ], + ) + fig, ax = plt.subplots(figsize=(9, 5.5)) + index = np.arange(n) + ax.scatter( + index, + y_multi, + s=24, + facecolors="none", + edgecolors="steelblue", + linewidths=1.2, + label="Observed", + ) + ax.scatter( + index[missing], imputed_multi, s=26, marker="s", color="red", label="Imputed (NNS 1-NN)" + ) + ax.set_xlabel("Observation index") + ax.set_ylabel("y") + ax.set_title("NNS 1-NN Multivariate Imputation") + ax.legend(loc="upper left", frameon=False) + save_figure(fig, OUT, "09_multivariate_imputation.png") + note( + "The R vignette closes with bootstrap multiple imputation pooled via " + "Rubin's rules; that section is narrative in R (no executable chunk) and " + "links to the NNS_MI_vs_MICE comparison, so it is summarized here rather " + "than reinvented." + ) + + +if __name__ == "__main__": + main() diff --git a/examples/vignettes/08_classification.py b/examples/vignettes/08_classification.py index 6c0a4c5..8e46f5c 100644 --- a/examples/vignettes/08_classification.py +++ b/examples/vignettes/08_classification.py @@ -17,13 +17,11 @@ # continuous model, with the response encoded as classes beginning at 1. The R # vignette emphasizes that these are multidimensional partitions rather than a # sequence of one-variable tree splits. - import matplotlib.pyplot as plt import numpy as np from examples._vignette_support import ( fast_mode, - gap, load_iris, note, output_dir, @@ -89,10 +87,8 @@ def main() -> None: # ## NNS partitions used by the classifier # # The R vignette prints ``NNS.reg(... )$rhs.partitions`` for all four - # predictors. Python currently returns the fitted reduction and equation but - # not the internal per-predictor partition list. We therefore expose the - # closest public diagnostic explicitly instead of pretending the object is - # present. + # predictors; the Python multivariate result returns the same per-predictor + # partition list. subsection("NNS partitions") base_fit = nns_reg( iris_x, @@ -103,11 +99,7 @@ def main() -> None: class_levels=levels, ) show("Classification regression result", base_fit) - gap( - "R exposes NNS.reg(... )$rhs.partitions for every predictor. The Python " - "public result does not yet expose that internal list. The two-dimensional " - "NNS.part result above is the supported public partition diagnostic." - ) + show("Per-predictor partitions (rhs.partitions)", base_fit["rhs.partitions"]) # %% [markdown] # ## NNS.boost() @@ -141,7 +133,9 @@ def main() -> None: axes[0].set_xticks(np.arange(len(weight_values)), labels_for_plot, rotation=35, ha="right") axes[0].set_title("Boost feature weights") axes[1].bar(np.arange(len(frequency_values)), frequency_values) - axes[1].set_xticks(np.arange(len(frequency_values)), list(frequencies.keys()), rotation=35, ha="right") + axes[1].set_xticks( + np.arange(len(frequency_values)), list(frequencies.keys()), rotation=35, ha="right" + ) axes[1].set_title("Boost feature frequency") save_figure(fig, OUT, "02_boost_feature_diagnostics.png") diff --git a/examples/vignettes/09_forecasting.py b/examples/vignettes/09_forecasting.py index a3d8817..8ba29f0 100644 --- a/examples/vignettes/09_forecasting.py +++ b/examples/vignettes/09_forecasting.py @@ -16,7 +16,6 @@ # lagged seasonal components. The R vignette moves from a fixed seasonal factor # to validation-based selection, automatic seasonality, and finally # ``NNS.ARMA.optim``. - import matplotlib.pyplot as plt import numpy as np @@ -108,7 +107,9 @@ def main() -> None: linear_rmse = rmse(linear, actual) print(f"Linear seasonal-factor-12 RMSE: {linear_rmse:.6f}") show("Linear forecast", linear_raw) - forecast_figure(training, actual, linear, "Linear NNS.ARMA, seasonal factor 12", "01_linear.png") + forecast_figure( + training, actual, linear, "Linear NNS.ARMA, seasonal factor 12", "01_linear.png" + ) # %% [markdown] # ## Nonlinear regression diff --git a/examples/vignettes/PARITY.md b/examples/vignettes/PARITY.md index 60b28af..f2fa5ba 100644 --- a/examples/vignettes/PARITY.md +++ b/examples/vignettes/PARITY.md @@ -69,12 +69,14 @@ Preserves partition orders 1–4, X-only partitions, clusters used in regression univariate and full-grid multivariate regression, interpolation/extrapolation, dimension reduction, thresholding, classification, `NNS.stack`, duplicated predictor dimensions, smoothing, and univariate/multivariate imputation. The -exact iris data are committed locally. Six figure groups and important model +exact iris data are committed locally. Nine figure groups and important model objects are produced. -**Known gaps:** Python does not yet expose R's Voronoi tessellation side effect or -`NNS.reg(... )$rhs.partitions`; both are labeled in the relevant sections and -supported public partition diagnostics are shown instead. +**Known gaps:** Python does not yet expose R's Voronoi tessellation side effect; +the partitions are drawn directly from the returned quadrant memberships instead. +The multivariate `rhs.partitions` and `RPM` structures are exposed and displayed. +R's closing bootstrap multiple-imputation section is narrative (no executable +chunk) and is summarized rather than reinvented. ## 08 — `08_classification.py` @@ -83,9 +85,9 @@ Preserves the splits-versus-partitions narrative, exact iris holdout rows `NNS.stack` cross-validation, and the depth/nearest-neighbor/extreme controls. Feature diagnostics, holdout paths, and returned model structures are displayed. -**Known gaps:** Python does not yet expose `rhs.partitions`. Balanced boosting -has a documented bit-for-bit RNG ordering gap from R's interleaved class draws; -the procedure and diagnostics remain equivalent. +**Known gaps:** Balanced boosting has a documented bit-for-bit RNG ordering gap +from R's interleaved class draws; the procedure and diagnostics remain +equivalent. The multivariate `rhs.partitions` structure is exposed and displayed. ## 09 — `09_forecasting.py` diff --git a/pyproject.toml b/pyproject.toml index b7f1cab..3b674e1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -110,6 +110,9 @@ select = ["E", "F", "I", "B", "UP", "N", "RUF", "TID"] [tool.ruff.lint.per-file-ignores] "tests/**" = ["TID251"] +# Instructional vignette scripts insert the repository root on sys.path before +# importing shared support, so module-level imports cannot all sit at the top. +"examples/**" = ["E402"] "scripts/generate_api_reference.py" = ["E501"] "scripts/regenerate_r_cache.py" = ["TID251"] "scripts/install_local_r_nns.py" = ["TID251"]