diff --git a/docs/api_reference.md b/docs/api_reference.md index 7575285..e8d6bf4 100644 --- a/docs/api_reference.md +++ b/docs/api_reference.md @@ -315,7 +315,11 @@ NNS ANOVA-style comparison helper covering binary, multi-group, pairwise, and de Closest R API: `NNS.norm`. -NNS normalization helper for numeric matrix-style inputs. +NNS normalization helper. Accepts a 2-D matrix or, like R's list input, a +sequence of 1-D vectors (one per variable). Equal-length vectors are +column-stacked and normalized through the matrix path; unequal-length +vectors force `linear=True` exactly as R does and return a list of scaled +arrays. #### `nns_distance` @@ -361,13 +365,17 @@ Rescales inputs using NNS conventions. Closest R APIs: `NNS.FSD`, `NNS.SSD`, and `NNS.TSD`. -Compute first-, second-, and third-order stochastic dominance. +Compute first-, second-, and third-order stochastic dominance. The two samples +may differ in length; curves are evaluated on the merged threshold grid exactly +as the R functions do. #### `fsd_uni`, `ssd_uni`, `tsd_uni` Closest R APIs: `NNS.FSD.uni`, `NNS.SSD.uni`, and `NNS.TSD.uni`. -Univariate wrappers for stochastic dominance workflows. +Univariate wrappers for stochastic dominance workflows. Unlike R's C++ `.uni` +routines, unequal-length samples are supported with the same merged-grid +semantics as the pairwise tests. #### `nns_sd_cluster` diff --git a/docs/api_status.md b/docs/api_status.md index 24ab029..85ea464 100644 --- a/docs/api_status.md +++ b/docs/api_status.md @@ -57,7 +57,7 @@ invariant, and property coverage. | Bootstrap/Monte Carlo: `nns_meboot`, `nns_mc` | implemented | medium | Deterministic diagnostics are parity-tested; exact stochastic replicate parity with R is not expected. | | Stochastic dominance/superiority: `fsd`, `ssd`, `tsd`, `.uni` wrappers, `nns_ss`, `nns_sd_cluster`, `sd_efficient_set` | implemented | medium | Public structures and deterministic paths are covered. SD uses exact pure-NumPy prefix-pair kernels plus a degree-1 discrete order-statistic matrix path; R's C++ core remains faster on full finance fixtures. Stochastic intervals use NNS Python RNG. | | ANOVA: `nns_anova` | implemented | high | Binary, multi-group, pairwise, and degenerate `NaN` conventions are covered. | -| Normalization: `nns_norm` | implemented | high | Numeric matrix path is implemented. | +| Normalization: `nns_norm` | implemented | high | Numeric matrix path and R's list-of-vectors path are implemented; unequal-length vectors force linear scaling as in R. | | Categorical helpers: `encode_factor_codes`, `factor_2_dummy`, `factor_2_dummy_fr`, `prepare_factor_predictors` | implemented | high | Explicit `levels=` / `factor_levels=` should be used to reproduce R factor ordering. `prepare_factor_predictors(...)` exposes the regression-ready full-rank design matrix path. | | Scalar differentiation: `nns_diff`, `dy_dx` | implemented | high | `dy_dx(..., eval_point="overall")` and numeric evaluation points are covered. | | Multivariate differentiation: `dy_d` | partial | medium-high | Scalar and vectorized point/distribution modes are covered on focused fixtures. Mixed derivatives are supported for two-regressor inputs where defined; multi-row matrix mixed derivatives use pointwise Python semantics rather than R's order-dependent list-matrix packing quirk. | diff --git a/docs/conventions.md b/docs/conventions.md index 4f15d64..8c2c838 100644 --- a/docs/conventions.md +++ b/docs/conventions.md @@ -73,7 +73,14 @@ NNS Python does not plot the dendrogram; it only returns the object data. The stochastic-dominance implementation is deliberately pure NumPy. It mirrors R's C++ SD core mathematically by sorting each column once, storing prefix sums, and evaluating dominance on each pair's merged threshold grid rather than on one -global all-column grid. The full prefix-pair dominance matrix remains available +global all-column grid. Pairwise tests (`fsd`, `ssd`, `tsd`, and the `*_uni` +wrappers) accept samples of unequal length: each sample's curve is evaluated on +the merged grid exactly as R's `NNS.FSD`/`NNS.SSD`/`NNS.TSD` compute +`LPM(degree, sort(c(x, y)), sample)` per sample. R's C++ `.uni` walkers assume +equal-length inputs, so for unequal lengths the Python `*_uni` wrappers are an +intentional extension carrying the same merged-grid semantics. The +matrix-based efficient-set and cluster routines operate on data columns and +therefore remain equal-length by construction. The full prefix-pair dominance matrix remains available internally for verification and fallback. Large degree-1 discrete calls use an exact order-statistic dominance matrix: with equal-length empirical samples, one sample first-order stochastically dominates another exactly when every diff --git a/src/nns/norm.py b/src/nns/norm.py index 07dc3f7..dfe77b0 100644 --- a/src/nns/norm.py +++ b/src/nns/norm.py @@ -1,6 +1,7 @@ from __future__ import annotations -from typing import cast +from collections.abc import Sequence +from typing import cast, overload import numpy as np from numpy.typing import NDArray @@ -8,9 +9,44 @@ from nns.dependence import nns_dep -def nns_norm(x: NDArray[np.float64], linear: bool = False) -> NDArray[np.float64]: - """Normalize a numeric matrix following R's NNS.norm scaling.""" - values = _as_matrix(x) +@overload +def nns_norm(x: NDArray[np.float64], linear: bool = ...) -> NDArray[np.float64]: ... + + +@overload +def nns_norm( + x: Sequence[NDArray[np.float64]], + linear: bool = ..., +) -> NDArray[np.float64] | list[NDArray[np.float64]]: ... + + +def nns_norm( + x: NDArray[np.float64] | Sequence[NDArray[np.float64]], + linear: bool = False, +) -> NDArray[np.float64] | list[NDArray[np.float64]]: + """Normalize variables following R's NNS.norm scaling. + + Two input conventions are supported, matching R's ``NNS.norm(X, ...)``: + + * A 2-D array whose columns are variables. Returns the scaled 2-D array. + * A list or tuple of 1-D arrays, one per variable (R's list input; the + elements are variables/columns, not observation rows). Equal-length + vectors are column-stacked and normalized through the matrix path, + returning a 2-D array — mirroring R, where ``mapply`` simplifies the + equal-length list result to a matrix. Unequal-length vectors force + ``linear=True`` exactly as R does (dependence-based scale factors need + aligned columns) and return a list of scaled arrays, one per input + vector. + """ + if isinstance(x, np.ndarray): + values = _as_matrix(x) + else: + series = [_as_vector(item, index) for index, item in enumerate(x)] + if not series: + raise ValueError("x must be non-empty.") + if len({item.size for item in series}) > 1: + return _norm_unequal_series(series) + values = _as_matrix(np.column_stack(series)) means = np.mean(values, axis=0) means = means.copy() means[means == 0.0] = 1e-10 @@ -30,6 +66,14 @@ def nns_norm(x: NDArray[np.float64], linear: bool = False) -> NDArray[np.float64 return scaled +def _norm_unequal_series(series: list[NDArray[np.float64]]) -> list[NDArray[np.float64]]: + means = np.array([float(np.mean(item)) for item in series]) + means[means == 0.0] = 1e-10 + ratio_grid = means[:, np.newaxis] * (1.0 / means[np.newaxis, :]) + scales = np.mean(ratio_grid, axis=0) + return [item * scale for item, scale in zip(series, scales, strict=True)] + + def _scale_factor(values: NDArray[np.float64]) -> NDArray[np.float64]: if values.shape[1] < 10: return cast(NDArray[np.float64], np.abs(np.corrcoef(values, rowvar=False))) @@ -44,6 +88,17 @@ def _scale_factor(values: NDArray[np.float64]) -> NDArray[np.float64]: return deps +def _as_vector(x: NDArray[np.float64], index: int) -> NDArray[np.float64]: + values = np.asarray(x, dtype=np.float64) + if values.ndim != 1: + raise ValueError(f"x[{index}] must be a 1D numeric vector.") + if values.size == 0: + raise ValueError(f"x[{index}] must be non-empty.") + if not np.all(np.isfinite(values)): + raise ValueError(f"x[{index}] must contain only finite values.") + return values + + def _as_matrix(x: NDArray[np.float64]) -> NDArray[np.float64]: values = np.asarray(x, dtype=np.float64) if values.ndim != 2: diff --git a/src/nns/stochastic_dominance.py b/src/nns/stochastic_dominance.py index fa0a836..5cc4537 100644 --- a/src/nns/stochastic_dominance.py +++ b/src/nns/stochastic_dominance.py @@ -4,6 +4,12 @@ curve equality guard: equal LPM/CDF curves are non-dominance even when samples differ below meaningful double precision. Efficient-set output follows the R C++ routine's LPM-at-global-maximum ordering and original-index tie break. + +Pairwise tests accept samples of unequal length: each sample's LPM/CDF curve is +evaluated on the merged threshold grid, exactly as R's NNS.FSD/NNS.SSD/NNS.TSD +compute ``LPM(degree, sort(c(x, y)), sample)`` per sample. (R's C++ ``.uni`` +walkers assume equal lengths; the Python ``*_uni`` wrappers extend the same +merged-grid semantics to unequal lengths.) """ from __future__ import annotations @@ -55,7 +61,7 @@ class _SDOrderStatPrecomputed: def fsd(x: NDArray[np.float64], y: NDArray[np.float64]) -> int: - """First-order stochastic dominance.""" + """First-order stochastic dominance; ``x`` and ``y`` may differ in length.""" x_values = _as_sd_values(x, "x") y_values = _as_sd_values(y, "y") return _sd_result(x_values, y_values, 1) @@ -70,7 +76,7 @@ def fsd_uni(x: NDArray[np.float64], y: NDArray[np.float64], type: str = "discret def ssd(x: NDArray[np.float64], y: NDArray[np.float64]) -> int: - """Second-order stochastic dominance.""" + """Second-order stochastic dominance; ``x`` and ``y`` may differ in length.""" x_values = _as_sd_values(x, "x") y_values = _as_sd_values(y, "y") return _sd_result(x_values, y_values, 2) @@ -84,7 +90,7 @@ def ssd_uni(x: NDArray[np.float64], y: NDArray[np.float64]) -> int: def tsd(x: NDArray[np.float64], y: NDArray[np.float64]) -> int: - """Third-order stochastic dominance.""" + """Third-order stochastic dominance; ``x`` and ``y`` may differ in length.""" x_values = _as_sd_values(x, "x") y_values = _as_sd_values(y, "y") return _sd_result(x_values, y_values, 3) @@ -754,8 +760,6 @@ def _dominates_uni( *, discrete: bool, ) -> bool: - if x.size != y.size: - raise ValueError("x and y must have the same length.") if np.array_equal(np.sort(x), np.sort(y)): return False if np.min(x) < np.min(y): diff --git a/tests/invariants/test_norm.py b/tests/invariants/test_norm.py index 7f4c001..61fa7d1 100644 --- a/tests/invariants/test_norm.py +++ b/tests/invariants/test_norm.py @@ -34,3 +34,60 @@ def test_nonlinear_nns_norm_preserves_shape_for_wide_matrix() -> None: assert result.shape == x.shape assert np.all(np.isfinite(result)) + + +def test_equal_length_sequence_matches_matrix_path() -> None: + x = np.linspace(1.0, 3.0, 50) + y = np.linspace(2.0, 8.0, 50) + z = np.linspace(10.0, 20.0, 50) + matrix = np.column_stack((x, y, z)) + + for linear in (False, True): + from_sequence = nns_norm([x, y, z], linear=linear) + from_matrix = nns_norm(matrix, linear=linear) + assert isinstance(from_sequence, np.ndarray) + np.testing.assert_allclose(from_sequence, from_matrix) + + +def test_unequal_length_sequence_forces_linear_scaling() -> None: + vec1 = np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]) + vec2 = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0]) + vec3 = np.array([0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3]) + + result = nns_norm([vec1, vec2, vec3]) + + assert isinstance(result, list) + assert [item.size for item in result] == [7, 6, 9] + + # Linear scaling equalizes every scaled mean at the grand mean of means, + # and the linear flag is forced regardless of its passed value (as in R). + grand_mean = np.mean([vec1.mean(), vec2.mean(), vec3.mean()]) + for item in result: + np.testing.assert_allclose(item.mean(), grand_mean) + forced = nns_norm([vec1, vec2, vec3], linear=True) + for got, expected in zip(result, forced, strict=True): + np.testing.assert_allclose(got, expected) + + +def test_zero_sum_length_differences_detected_as_unequal() -> None: + # Lengths (5, 6, 5) have pairwise diffs summing to zero; they must still + # take the unequal-length path rather than being treated as a matrix. + vec1 = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + vec2 = np.array([10.0, 20.0, 30.0, 40.0, 50.0, 60.0]) + vec3 = np.array([0.5, 0.6, 0.7, 0.8, 0.9]) + + result = nns_norm([vec1, vec2, vec3]) + + assert isinstance(result, list) + assert [item.size for item in result] == [5, 6, 5] + + +def test_sequence_input_validation() -> None: + import pytest + + with pytest.raises(ValueError, match="non-empty"): + nns_norm([]) + with pytest.raises(ValueError, match=r"x\[1\] must be a 1D"): + nns_norm([np.ones(3), np.ones((3, 2))]) + with pytest.raises(ValueError, match=r"x\[0\] must contain only finite"): + nns_norm([np.array([1.0, np.nan]), np.ones(3)]) diff --git a/tests/invariants/test_stochastic_dominance.py b/tests/invariants/test_stochastic_dominance.py index 3986f8d..8a738ba 100644 --- a/tests/invariants/test_stochastic_dominance.py +++ b/tests/invariants/test_stochastic_dominance.py @@ -1,8 +1,32 @@ from __future__ import annotations import numpy as np +import pytest +from numpy.typing import NDArray -from nns import fsd, ssd, tsd +from nns import fsd, fsd_uni, lpm, ssd, ssd_uni, tsd, tsd_uni + + +def _r_sd_reference(x: NDArray[np.float64], y: NDArray[np.float64], degree: int) -> int: + """R's NNS.FSD/NNS.SSD/NNS.TSD decision rule, transcribed literally.""" + grid = np.sort(np.concatenate((x, y))) + if degree == 1: + # LPM.ratio(0, grid, sample) == LPM(0, grid, sample) == ECDF + curve_x = np.asarray(lpm(0, grid, x), dtype=np.float64) + curve_y = np.asarray(lpm(0, grid, y), dtype=np.float64) + else: + curve_x = np.asarray(lpm(degree - 1, grid, x), dtype=np.float64) + curve_y = np.asarray(lpm(degree - 1, grid, y), dtype=np.float64) + + mean_ok_xy = degree == 1 or float(np.mean(x)) >= float(np.mean(y)) + mean_ok_yx = degree == 1 or float(np.mean(y)) >= float(np.mean(x)) + curves_identical = np.array_equal(curve_x, curve_y) + + if not np.any(curve_x > curve_y) and x.min() >= y.min() and mean_ok_xy and not curves_identical: + return 1 + if not np.any(curve_y > curve_x) and y.min() >= x.min() and mean_ok_yx and not curves_identical: + return -1 + return 0 def test_sd_antisymmetry() -> None: @@ -29,3 +53,60 @@ def test_self_does_not_dominate() -> None: assert fsd(x, x) == 0 assert ssd(x, x) == 0 assert tsd(x, x) == 0 + + +def test_unequal_length_shifted_samples_dominate() -> None: + x = np.array([2.0, 3.0, 4.0, 5.0, 6.0]) + y = np.array([1.0, 2.0, 3.0]) + + assert fsd(x, y) == 1 + assert ssd(x, y) == 1 + assert tsd(x, y) == 1 + assert fsd(y, x) == -1 + + assert fsd_uni(x, y) == 1 + assert fsd_uni(x, y, "continuous") == 1 + assert ssd_uni(x, y) == 1 + assert tsd_uni(x, y) == 1 + assert fsd_uni(y, x) == 0 + assert ssd_uni(y, x) == 0 + assert tsd_uni(y, x) == 0 + + +def test_unequal_length_antisymmetry() -> None: + rng = np.random.default_rng(7) + x = rng.normal(1.0, 1.0, 37) + y = rng.normal(0.0, 1.0, 61) + + assert fsd(x, y) == -fsd(y, x) + assert ssd(x, y) == -ssd(y, x) + assert tsd(x, y) == -tsd(y, x) + + +def test_mtcars_transmission_groups_match_r() -> None: + # mtcars mpg split by transmission: R's NNS.FSD returns "X FSD Y" for + # (manual, auto) despite the samples having different lengths (13 vs 19). + auto_mpg = np.array( + [21.4, 18.7, 18.1, 14.3, 24.4, 22.8, 19.2, 17.8, 16.4, 17.3, + 15.2, 10.4, 10.4, 14.7, 21.5, 15.5, 15.2, 13.3, 19.2] + ) + manual_mpg = np.array( + [21.0, 21.0, 22.8, 32.4, 30.4, 33.9, 27.3, 26.0, 30.4, 15.8, 19.7, 15.0, 21.4] + ) + + assert fsd(manual_mpg, auto_mpg) == 1 + assert fsd(auto_mpg, manual_mpg) == -1 + assert fsd_uni(manual_mpg, auto_mpg) == 1 + + +@pytest.mark.parametrize("degree", [1, 2, 3]) +@pytest.mark.parametrize("seed", range(8)) +def test_unequal_length_matches_r_decision_rule(degree: int, seed: int) -> None: + rng = np.random.default_rng(seed) + sizes = rng.integers(3, 60, size=2) + shift = rng.uniform(-0.5, 0.5) + x = rng.normal(shift, 1.0, int(sizes[0])) + y = rng.normal(0.0, rng.uniform(0.5, 1.5), int(sizes[1])) + + function = {1: fsd, 2: ssd, 3: tsd}[degree] + assert function(x, y) == _r_sd_reference(x, y, degree)