From 97dd09a619139f520017f2ace533947bc9a159e1 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 17 Jul 2026 17:45:44 +0000 Subject: [PATCH 1/2] diff: add dy_d_best, the original v0.5.7 dy.d_ finite-difference estimator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds nns.dy_d_best, a faithful Python port of the original NNS 0.5.7 dy.d_ used to produce the results in Vinod & Viole (2020), SSRN 3681436. It differs from the shipped dy_d (which sources f(x +/- h) from nns_stack(method=(1, 2))): it uses nns_reg(dim_red_method='equal', point_only=True) directly (no smoothing), a quantile-spacing step, a degree-0 quantile grid, plain-mean aggregation, distance = 2*h_step, and a plain nanmean blend across the v0.5.7 bandwidths h_s = 1/log(size(x), [2, 10]) extended by 10*h_s and doubled when nns_dep(x[:, wrt], y) < 0.5. The wrapper is byte-for-byte the same algorithm as the restored R NNS.dy.d_ (OVVO-Financial/NNS), so R and Python agree given the same NNS.reg engine; against the current engine both differ from the 0.5.7 paper values because the regression engine itself changed since 0.5.7. Adds invariant tests (export, shapes, single-wrt 1-D output, mixed, determinism, pinned golden values). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.8 --- src/nns/__init__.py | 1 + src/nns/diff.py | 276 ++++++++++++++++++++++++++++++++++ tests/invariants/test_diff.py | 77 ++++++++++ 3 files changed, 354 insertions(+) diff --git a/src/nns/__init__.py b/src/nns/__init__.py index ba5b401..1d8f8ce 100644 --- a/src/nns/__init__.py +++ b/src/nns/__init__.py @@ -26,6 +26,7 @@ "d_lpm": ("nns.co_moments", "d_lpm"), "dpm_nd": ("nns.dependence", "dpm_nd"), "dy_d": ("nns.diff", "dy_d"), + "dy_d_best": ("nns.diff", "dy_d_best"), "dy_dx": ("nns.diff", "dy_dx"), "d_upm": ("nns.co_moments", "d_upm"), "ecdf_pm": ("nns.classical", "ecdf_pm"), diff --git a/src/nns/diff.py b/src/nns/diff.py index fdf397f..d8c2353 100644 --- a/src/nns/diff.py +++ b/src/nns/diff.py @@ -262,6 +262,282 @@ def _combine_dy_d_outputs( } +def dy_d_best( + x: NDArray[Any], + y: NDArray[Any], + wrt: int | NDArray[np.int64], + eval_points: str | float | NDArray[np.float64] = "obs", + *, + mixed: bool = False, + messages: bool = False, + factor_levels: Sequence[Sequence[Any] | None] | None = None, +) -> dict[str, NDArray[np.float64]]: + """Original v0.5.7 ``dy.d_`` finite-difference estimator (base-R faithful). + + This restores the exact finite-difference design used to produce the results + in Viole (2020), "Comparing Old and New Partial Derivative Estimates from + Nonlinear Nonparametric Regressions" (SSRN 3681436). It is a straight port of + the original ``data.table`` implementation to NumPy, so it carries none of the + later modifications (no ``gravity`` aggregation, no smoothing spline, no + descending-band weighting). Specifically it uses: + + * bandwidths ``h_s = 1/log(size(x), [2, 10])`` extended by ``10 * h_s`` and + doubled when ``nns_dep(x[:, wrt], y) < 0.5``; + * a quantile-spacing step ``h_step = |mean(diff(VaR(seq(.01, 1, h), 0, x)))|``; + * a degree-0 quantile evaluation grid ``VaR(seq(0, 1, .05), 0, .)``; + * ``nns_reg(dim_red_method="equal", smooth=False, point_only=True)`` estimates; + * **plain mean** aggregation over the grid, ``distance = 2 * h_step`` so that + ``First = (upper - lower) / (2 h_step)`` and + ``Second = (upper - 2 f(x) + lower) / (2 h_step) ** 2``; + * a cumulative +/- evaluation window across bandwidths; and + * a plain ``rowMeans`` (``nanmean``) blend across bandwidths. + + Values from the very largest bandwidths (where ``seq(.01, 1, h)`` collapses to + a single probability) evaluate to NaN and are dropped by the ``nanmean`` blend, + matching the R ``rowMeans(..., na.rm = TRUE)`` behaviour. + + ``wrt`` uses R-style 1-based indexing into the factor-expanded predictor + matrix. A scalar/1-D ``eval_points`` evaluates only the selected regressor; a + 2-D array evaluates complete predictor tuples. + """ + raw_x = np.asarray(x) + if raw_x.ndim != 2: + raise ValueError("Please ensure (x) is a matrix or data.frame type object.") + if raw_x.shape[1] < 2: + raise ValueError("Please use NNS::dy.dx(...) for univariate partial derivatives.") + + raw_y = np.asarray(y).reshape(-1) + if raw_y.size != raw_x.shape[0]: + raise ValueError("x and y must have compatible row counts.") + if _dy_d_has_missing(raw_x) or _dy_d_has_missing(raw_y): + raise ValueError("You have some missing values, please address.") + + x_values = _dy_d_expand_predictors(x, factor_levels=factor_levels) + y_values = np.asarray(raw_y, dtype=np.float64) + if _dy_d_has_missing(x_values) or _dy_d_has_missing(y_values): + raise ValueError("You have some missing values, please address.") + + wrt_values = _dy_d_validate_wrt(wrt, x_values.shape[1]) + outputs = [ + _dy_d_best_scalar( + x_values, + y_values, + int(wrt_value) - 1, + eval_points, + mixed=bool(mixed), + messages=bool(messages), + wrt_label=int(wrt_value), + ) + for wrt_value in wrt_values + ] + if len(outputs) == 1: + return outputs[0] + return _combine_dy_d_outputs(outputs) + + +def _dy_d_best_reg_estimates( + x: NDArray[np.float64], + y: NDArray[np.float64], + test_points: NDArray[np.float64], +) -> NDArray[np.float64]: + """f(x +/- h) via NNS.reg(dim.red.method='equal', point.only=TRUE) (no smoothing).""" + from nns.regression import nns_reg + + result = nns_reg( + x, + y, + point_est=np.asarray(test_points, dtype=np.float64), + dim_red_method="equal", + threshold=0.0, + order=None, + point_only=True, + smooth=False, + ) + return np.asarray(result["Point.est"], dtype=np.float64).reshape(-1) + + +def _r_seq(start: float, stop: float, by: float) -> NDArray[np.float64]: + """Reproduce R's seq(start, stop, by) (empty diff when a single point results).""" + if not np.isfinite(by) or by == 0.0: + return np.asarray([start], dtype=np.float64) + count = int(np.floor((stop - start) / by + 1e-9)) + return start + by * np.arange(0, count + 1, dtype=np.float64) + + +def _v057_bandwidths(x_size: int, wrt_dependence: float) -> NDArray[np.float64]: + """h_s = 1/log(length(x), c(2, 10)); c(h_s, 10*h_s); doubled if dependence < .5.""" + base = np.log(np.asarray([2.0, 10.0])) / np.log(float(x_size)) + h_s = np.concatenate([base, 10.0 * base]) + if wrt_dependence < 0.5: + h_s = 2.0 * h_s + return h_s + + +def _dy_d_best_scalar( + x_values: NDArray[np.float64], + y_values: NDArray[np.float64], + wrt_index: int, + eval_points: str | float | NDArray[np.float64], + *, + mixed: bool, + messages: bool, + wrt_label: int, +) -> dict[str, NDArray[np.float64]]: + from nns.dependence import nns_dep + from nns.var import lpm_var + + _n_rows, n_predictors = x_values.shape + if wrt_index < 0 or wrt_index >= n_predictors: + raise ValueError("`wrt` must select exactly one column of the expanded predictor matrix.") + if n_predictors != 2: + mixed = False + + if messages: + print( + "Currently generating NNS.reg finite difference estimates...Regressor " + f"{wrt_label}\r" + ) + + eval_values, vector_branch = _dy_d_eval_points(x_values, wrt_index, eval_points) + + dependence = float(nns_dep(x_values[:, wrt_index], y_values)["Dependence"]) + h_s = _v057_bandwidths(x_values.size, dependence) + + def _quantile_step(column: int, band: float) -> float: + probs = _r_seq(0.01, 1.0, band) + quantiles = np.asarray( + [lpm_var(float(p), 0.0, x_values[:, column]) for p in probs], dtype=np.float64 + ) + if quantiles.size < 2: + return float("nan") + return float(abs(np.mean(np.diff(quantiles)))) + + if vector_branch: + eval_vec = np.asarray(eval_values, dtype=np.float64).reshape(-1) + window_min = eval_vec.copy() + window_max = eval_vec.copy() + grid = np.column_stack( + [ + np.asarray( + [lpm_var(float(p), 0.0, x_values[:, col]) for p in _r_seq(0.0, 1.0, 0.05)], + dtype=np.float64, + ) + for col in range(n_predictors) + ] + ) + sample_size = grid.shape[0] + else: + eval_mat = _as_eval_matrix(eval_values, n_predictors) + window_min = eval_mat.copy() + window_max = eval_mat.copy() + + firsts: list[NDArray[np.float64]] = [] + seconds: list[NDArray[np.float64]] = [] + mixeds: list[NDArray[np.float64]] = [] + + for index, band in enumerate(h_s, start=1): + h_step = _quantile_step(wrt_index, float(band)) + finite = np.isfinite(h_step) and h_step != 0.0 + distance = 2.0 * h_step + + if vector_branch: + k = eval_vec.size + if finite: + window_min = window_min - h_step + window_max = window_max + h_step + deriv_points = np.vstack([grid] * (3 * k)) + wrt_column = np.repeat( + np.ravel(np.vstack((window_min, eval_vec, window_max)).T), sample_size + )[: deriv_points.shape[0]] + deriv_points[:, wrt_index] = wrt_column + position = np.tile( + np.repeat(np.asarray(["l", "m", "u"], dtype=object), sample_size), k + )[: deriv_points.shape[0]] + ids = np.repeat(np.arange(k), 3 * sample_size)[: deriv_points.shape[0]] + if messages: + print(f"Currently evaluating the {deriv_points.shape[0]} required points\r") + estimates = _dy_d_best_reg_estimates(x_values, y_values, deriv_points) + lower = np.asarray( + [np.mean(estimates[(position == "l") & (ids == g)]) for g in range(k)] + ) + two_fx = 2.0 * np.asarray( + [np.mean(estimates[(position == "m") & (ids == g)]) for g in range(k)] + ) + upper = np.asarray( + [np.mean(estimates[(position == "u") & (ids == g)]) for g in range(k)] + ) + firsts.append((upper - lower) / distance) + seconds.append((upper - two_fx + lower) / distance**2) + else: + firsts.append(np.full(k, np.nan)) + seconds.append(np.full(k, np.nan)) + mixed_eval_points: NDArray[np.float64] | None = None + else: + n_eval = eval_mat.shape[0] + if finite: + window_min[:, wrt_index] = window_min[:, wrt_index] - h_step + window_max[:, wrt_index] = window_max[:, wrt_index] + h_step + deriv_points = np.vstack((window_min, eval_mat, window_max)) + if messages: + print( + "Currently generating NNS.reg finite difference estimates...bandwidth " + f"{index} of {h_s.size}\r" + ) + estimates = _dy_d_best_reg_estimates(x_values, y_values, deriv_points) + lower = estimates[:n_eval] + two_fx = 2.0 * estimates[n_eval : 2 * n_eval] + upper = estimates[2 * n_eval :] + firsts.append((upper - lower) / distance) + seconds.append((upper - two_fx + lower) / distance**2) + else: + firsts.append(np.full(n_eval, np.nan)) + seconds.append(np.full(n_eval, np.nan)) + mixed_eval_points = eval_mat + + if mixed: + if vector_branch: + tuple_eval = np.asarray(eval_values, dtype=np.float64).reshape(-1) + if tuple_eval.size != 2: + raise ValueError("Mixed Derivatives are only for 2 IV") + mixed_eval_points = tuple_eval.reshape(1, 2) + assert mixed_eval_points is not None + if mixed_eval_points.shape[1] != 2: + raise ValueError("Mixed Derivatives are only for 2 IV") + step_1 = _quantile_step(0, float(band)) + step_2 = _quantile_step(1, float(band)) + if np.isfinite(step_1) and np.isfinite(step_2) and step_1 != 0.0 and step_2 != 0.0: + corners = np.vstack( + [ + np.array( + [ + [p[0] + step_1, p[1] + step_2], + [p[0] - step_1, p[1] + step_2], + [p[0] + step_1, p[1] - step_2], + [p[0] - step_1, p[1] - step_2], + ] + ) + for p in mixed_eval_points + ] + ) + mixed_estimates = _dy_d_best_reg_estimates(x_values, y_values, corners) + z = mixed_estimates.reshape(-1, 4) + mixeds.append((z[:, 0] + z[:, 3] - z[:, 1] - z[:, 2]) / (4.0 * step_1 * step_2)) + else: + mixeds.append(np.full(mixed_eval_points.shape[0], np.nan)) + + def _row_nanmean(bands: list[NDArray[np.float64]]) -> NDArray[np.float64]: + matrix = np.column_stack([np.asarray(b, dtype=np.float64).reshape(-1) for b in bands]) + with np.errstate(invalid="ignore"): + return np.asarray(np.nanmean(matrix, axis=1), dtype=np.float64) + + output = {"First": _row_nanmean(firsts), "Second": _row_nanmean(seconds)} + if mixed: + output["Mixed"] = _row_nanmean(mixeds) + if messages: + print("\r") + return output + + def _dy_d_scalar( x_values: NDArray[np.float64], y_values: NDArray[np.float64], diff --git a/tests/invariants/test_diff.py b/tests/invariants/test_diff.py index 1c74520..de3faea 100644 --- a/tests/invariants/test_diff.py +++ b/tests/invariants/test_diff.py @@ -115,3 +115,80 @@ def test_dy_d_point_modes_preserve_linear_slope_direction() -> None: for eval_points in ("mean", "median", "last"): assert dy_d(x, positive, wrt=1, eval_points=eval_points)["First"][0] > 0.0 assert dy_d(x, negative, wrt=1, eval_points=eval_points)["First"][0] < 0.0 + + +def test_dy_d_best_is_exported() -> None: + import nns + + assert "dy_d_best" in nns.__all__ + assert callable(nns.dy_d_best) + + +def test_dy_d_best_vectorized_wrt_tuple_returns_first_second() -> None: + from nns import dy_d_best + + x = np.random.RandomState(0).randn(60, 3) + y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2] + + result = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3))) + + assert result.keys() == {"First", "Second"} + assert result["First"].shape == (1, 3) + assert result["Second"].shape == (1, 3) + assert np.all(np.isfinite(result["First"])) + + +def test_dy_d_best_single_wrt_returns_one_dimensional_first() -> None: + from nns import dy_d_best + + x = np.random.RandomState(0).randn(60, 3) + y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2] + + # A single `wrt` returns 1-D arrays (one value per evaluation point). + result = dy_d_best(x, y, wrt=1, eval_points=np.array([[0.0, 0.0, 0.0]])) + + assert result["First"].shape == (1,) + assert np.all(np.isfinite(result["First"])) + + +def test_dy_d_best_two_column_mixed_returns_mixed() -> None: + from nns import dy_d_best + + x = np.random.RandomState(1).randn(60, 2) + y = x[:, 0] * x[:, 1] + + result = dy_d_best(x, y, wrt=np.array([1, 2]), eval_points=np.array([[0.0, 0.0]]), mixed=True) + + assert result.keys() == {"First", "Second", "Mixed"} + assert result["First"].shape == (1, 2) + assert result["Mixed"].shape == (1, 2) + + +def test_dy_d_best_is_deterministic() -> None: + from nns import dy_d_best + + x = np.random.RandomState(0).randn(60, 3) + y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2] + + a = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3))) + b = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3))) + + assert np.allclose(a["First"], b["First"]) + assert np.allclose(a["Second"], b["Second"]) + + +def test_dy_d_best_matches_pinned_values() -> None: + # Golden values guard the v0.5.7 finite-difference design against regressions. + from nns import dy_d_best + + x = np.random.RandomState(0).randn(60, 3) + y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2] + + result = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3))) + + np.testing.assert_allclose( + np.ravel(result["First"]), + np.array([0.974410, 1.085425, 0.978869]), + rtol=0.0, + atol=1e-4, + ) From 93b529497922c2a0ba581f303e6dda1b32a89c64 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 17 Jul 2026 19:07:21 +0000 Subject: [PATCH 2/2] dy_d_best: reconcile v0.5.7 dy.d_ with the current engine (dy.dx step + NNS.stack) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Reworks dy_d_best from the raw-engine v0.5.7 port into the reconciled estimator: * h_step now shares the dy.dx() logic - a locally-adaptive step centred on the evaluation point's percentile, h_step = VaR(p+H,1,x) - VaR(p-H,1,x) with p = lpm_ratio(1, eval, x) - and drops the cumulative window. This is what makes the estimates uniform across identically-distributed regressors (previously 0.29-0.62 scatter on the 5-way product benchmark). * estimates now come from nns_stack on the equal-weight synthetic regressor X* via the increased-dimension trick cbind(X*, X*), with method=(1, 2), dim_red_method='equal', order='max', folds=5. The cross-validated n.best regularises the sharper current engine. On Example 2 of Vinod & Viole (2020, SSRN 3681436) this yields a tight, uniform First = [0.94, 0.95, 0.93, 0.92, 0.93] at (1,1,1,1,1) (vs the paper's ~0.65 on the retired 0.5.7 engine), and the mixed derivative for x1^2*x2^2 at (.5,.5) recovers ~1.0. Bandwidths, First=(u-l)/(2h), Second=(u-2f+l)/h^2 blended with nanmean across bandwidths. Golden test values updated. Byte-for-byte the same algorithm as the reconciled R dy.d_ (OVVO-Financial/NNS). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.8 --- src/nns/diff.py | 303 ++++++++++++++++++---------------- tests/invariants/test_diff.py | 2 +- 2 files changed, 161 insertions(+), 144 deletions(-) diff --git a/src/nns/diff.py b/src/nns/diff.py index d8c2353..781612f 100644 --- a/src/nns/diff.py +++ b/src/nns/diff.py @@ -272,33 +272,36 @@ def dy_d_best( messages: bool = False, factor_levels: Sequence[Sequence[Any] | None] | None = None, ) -> dict[str, NDArray[np.float64]]: - """Original v0.5.7 ``dy.d_`` finite-difference estimator (base-R faithful). - - This restores the exact finite-difference design used to produce the results - in Viole (2020), "Comparing Old and New Partial Derivative Estimates from - Nonlinear Nonparametric Regressions" (SSRN 3681436). It is a straight port of - the original ``data.table`` implementation to NumPy, so it carries none of the - later modifications (no ``gravity`` aggregation, no smoothing spline, no - descending-band weighting). Specifically it uses: - - * bandwidths ``h_s = 1/log(size(x), [2, 10])`` extended by ``10 * h_s`` and - doubled when ``nns_dep(x[:, wrt], y) < 0.5``; - * a quantile-spacing step ``h_step = |mean(diff(VaR(seq(.01, 1, h), 0, x)))|``; - * a degree-0 quantile evaluation grid ``VaR(seq(0, 1, .05), 0, .)``; - * ``nns_reg(dim_red_method="equal", smooth=False, point_only=True)`` estimates; - * **plain mean** aggregation over the grid, ``distance = 2 * h_step`` so that - ``First = (upper - lower) / (2 h_step)`` and - ``Second = (upper - 2 f(x) + lower) / (2 h_step) ** 2``; - * a cumulative +/- evaluation window across bandwidths; and - * a plain ``rowMeans`` (``nanmean``) blend across bandwidths. - - Values from the very largest bandwidths (where ``seq(.01, 1, h)`` collapses to - a single probability) evaluate to NaN and are dropped by the ``nanmean`` blend, - matching the R ``rowMeans(..., na.rm = TRUE)`` behaviour. + """Reconciled ``dy.d_`` partial-derivative estimator. + + This is the reconciliation of the original NNS 0.5.7 ``dy.d_`` (Vinod & Viole + 2020, SSRN 3681436) with the current regression engine. It keeps the v0.5.7 + finite-difference scaffolding but makes two changes that make the estimates + uniform across identically-distributed regressors and independent of the old + data.table machinery: + + * **Step (``h_step``) shares the ``dy.dx`` logic** - a locally-adaptive step + centred on the evaluation point's percentile, + ``h_step = VaR(p + H, 1, x) - VaR(p - H, 1, x)`` with + ``p = LPM.ratio(1, eval, x)`` - instead of a global quantile spacing with a + cumulative window. This is what removes the cross-regressor scatter. + * **Estimates come from** ``nns_stack`` on the equal-weight synthetic regressor + ``X*`` via the increased-dimension trick ``cbind(X*, X*)``, with + ``method=(1, 2)``, ``dim_red_method="equal"``, ``order="max"`` and + ``folds=5``. The stack's cross-validated ``n.best`` regularises the (sharper) + current engine back toward the paper's regime. + + Bandwidths follow v0.5.7: ``h_s = 1/log(size(x), [2, 10])`` extended by + ``10 * h_s`` and doubled when ``nns_dep(x[:, wrt], y) < 0.5``. First and second + derivatives use the central-difference forms + ``First = (upper - lower) / (2 h_step)`` and + ``Second = (upper - 2 f(x) + lower) / h_step ** 2`` (matching ``dy.dx``), + blended across bandwidths with a plain ``nanmean``. ``wrt`` uses R-style 1-based indexing into the factor-expanded predictor - matrix. A scalar/1-D ``eval_points`` evaluates only the selected regressor; a - 2-D array evaluates complete predictor tuples. + matrix. A scalar/1-D ``eval_points`` evaluates only the selected regressor + (averaged over the distribution of the other regressors); a 2-D array + evaluates complete predictor tuples. """ raw_x = np.asarray(x) if raw_x.ndim != 2: @@ -340,24 +343,35 @@ def _dy_d_best_reg_estimates( y: NDArray[np.float64], test_points: NDArray[np.float64], ) -> NDArray[np.float64]: - """f(x +/- h) via NNS.reg(dim.red.method='equal', point.only=TRUE) (no smoothing).""" - from nns.regression import nns_reg + """f(x +/- h) via NNS.stack on the equal-weight synthetic regressor X*. - result = nns_reg( - x, - y, - point_est=np.asarray(test_points, dtype=np.float64), + Reduces both the training design and the test points to the equal-weight + synthetic regressor X* = rowMeans(.), then estimates with the NNS + increased-dimension trick cbind(X*, X*) under + ``method=(1, 2), dim_red_method="equal", order="max", folds=5``. The + cross-validated ``n.best`` is what regularises the current engine. + """ + from nns.stack import nns_stack + + x_star = np.asarray(x, dtype=np.float64).mean(axis=1) + test_star = np.asarray(test_points, dtype=np.float64).mean(axis=1) + result = nns_stack( + ivs_train=np.column_stack([x_star, x_star]), + dv_train=y, + ivs_test=np.column_stack([test_star, test_star]), + method=(1, 2), dim_red_method="equal", - threshold=0.0, - order=None, - point_only=True, - smooth=False, + status=False, + order="max", + folds=5, + ncores=1, + dist=None, ) - return np.asarray(result["Point.est"], dtype=np.float64).reshape(-1) + return np.asarray(result["stack"], dtype=np.float64).reshape(-1) def _r_seq(start: float, stop: float, by: float) -> NDArray[np.float64]: - """Reproduce R's seq(start, stop, by) (empty diff when a single point results).""" + """Reproduce R's seq(start, stop, by).""" if not np.isfinite(by) or by == 0.0: return np.asarray([start], dtype=np.float64) count = int(np.floor((stop - start) / by + 1e-9)) @@ -373,6 +387,17 @@ def _v057_bandwidths(x_size: int, wrt_dependence: float) -> NDArray[np.float64]: return h_s +def _dydx_step(column: NDArray[np.float64], eval_value: float, band: float) -> float: + """Locally-adaptive dy.dx step: VaR(p + H, 1, x) - VaR(p - H, 1, x), p = CDF(eval).""" + from nns.core import lpm_ratio + from nns.var import lpm_var + + p = float(lpm_ratio(1.0, float(eval_value), column)) + upper = lpm_var(min(1.0, p + band), 1.0, column) + lower = lpm_var(max(0.0, p - band), 1.0, column) + return float(upper - lower) + + def _dy_d_best_scalar( x_values: NDArray[np.float64], y_values: NDArray[np.float64], @@ -400,22 +425,16 @@ def _dy_d_best_scalar( eval_values, vector_branch = _dy_d_eval_points(x_values, wrt_index, eval_points) - dependence = float(nns_dep(x_values[:, wrt_index], y_values)["Dependence"]) + column = x_values[:, wrt_index] + dependence = float(nns_dep(column, y_values)["Dependence"]) h_s = _v057_bandwidths(x_values.size, dependence) - def _quantile_step(column: int, band: float) -> float: - probs = _r_seq(0.01, 1.0, band) - quantiles = np.asarray( - [lpm_var(float(p), 0.0, x_values[:, column]) for p in probs], dtype=np.float64 - ) - if quantiles.size < 2: - return float("nan") - return float(abs(np.mean(np.diff(quantiles)))) + firsts: list[NDArray[np.float64]] = [] + seconds: list[NDArray[np.float64]] = [] + mixeds: list[NDArray[np.float64]] = [] if vector_branch: eval_vec = np.asarray(eval_values, dtype=np.float64).reshape(-1) - window_min = eval_vec.copy() - window_max = eval_vec.copy() grid = np.column_stack( [ np.asarray( @@ -426,104 +445,71 @@ def _quantile_step(column: int, band: float) -> float: ] ) sample_size = grid.shape[0] + k = eval_vec.size + for band in h_s: + steps = np.array([_dydx_step(column, ev, float(band)) for ev in eval_vec]) + blocks: list[NDArray[np.float64]] = [] + position: list[str] = [] + ids: list[int] = [] + for g in range(k): + lower = grid.copy() + middle = grid.copy() + upper = grid.copy() + lower[:, wrt_index] = eval_vec[g] - steps[g] + middle[:, wrt_index] = eval_vec[g] + upper[:, wrt_index] = eval_vec[g] + steps[g] + blocks.extend((lower, middle, upper)) + position.extend(["l"] * sample_size + ["m"] * sample_size + ["u"] * sample_size) + ids.extend([g] * (3 * sample_size)) + estimates = _dy_d_best_reg_estimates(x_values, y_values, np.vstack(blocks)) + pos = np.asarray(position, dtype=object) + idx = np.asarray(ids) + band_first = np.empty(k) + band_second = np.empty(k) + for g in range(k): + lo = np.mean(estimates[(pos == "l") & (idx == g)]) + mid = np.mean(estimates[(pos == "m") & (idx == g)]) + up = np.mean(estimates[(pos == "u") & (idx == g)]) + h = steps[g] + if np.isfinite(h) and h != 0.0: + band_first[g] = (up - lo) / (2.0 * h) + band_second[g] = (up - 2.0 * mid + lo) / (h**2) + else: + band_first[g] = np.nan + band_second[g] = np.nan + firsts.append(band_first) + seconds.append(band_second) + mixed_eval: NDArray[np.float64] | None = None else: eval_mat = _as_eval_matrix(eval_values, n_predictors) - window_min = eval_mat.copy() - window_max = eval_mat.copy() - - firsts: list[NDArray[np.float64]] = [] - seconds: list[NDArray[np.float64]] = [] - mixeds: list[NDArray[np.float64]] = [] - - for index, band in enumerate(h_s, start=1): - h_step = _quantile_step(wrt_index, float(band)) - finite = np.isfinite(h_step) and h_step != 0.0 - distance = 2.0 * h_step - - if vector_branch: - k = eval_vec.size - if finite: - window_min = window_min - h_step - window_max = window_max + h_step - deriv_points = np.vstack([grid] * (3 * k)) - wrt_column = np.repeat( - np.ravel(np.vstack((window_min, eval_vec, window_max)).T), sample_size - )[: deriv_points.shape[0]] - deriv_points[:, wrt_index] = wrt_column - position = np.tile( - np.repeat(np.asarray(["l", "m", "u"], dtype=object), sample_size), k - )[: deriv_points.shape[0]] - ids = np.repeat(np.arange(k), 3 * sample_size)[: deriv_points.shape[0]] - if messages: - print(f"Currently evaluating the {deriv_points.shape[0]} required points\r") - estimates = _dy_d_best_reg_estimates(x_values, y_values, deriv_points) - lower = np.asarray( - [np.mean(estimates[(position == "l") & (ids == g)]) for g in range(k)] - ) - two_fx = 2.0 * np.asarray( - [np.mean(estimates[(position == "m") & (ids == g)]) for g in range(k)] - ) - upper = np.asarray( - [np.mean(estimates[(position == "u") & (ids == g)]) for g in range(k)] - ) - firsts.append((upper - lower) / distance) - seconds.append((upper - two_fx + lower) / distance**2) - else: - firsts.append(np.full(k, np.nan)) - seconds.append(np.full(k, np.nan)) - mixed_eval_points: NDArray[np.float64] | None = None - else: - n_eval = eval_mat.shape[0] - if finite: - window_min[:, wrt_index] = window_min[:, wrt_index] - h_step - window_max[:, wrt_index] = window_max[:, wrt_index] + h_step - deriv_points = np.vstack((window_min, eval_mat, window_max)) - if messages: - print( - "Currently generating NNS.reg finite difference estimates...bandwidth " - f"{index} of {h_s.size}\r" - ) - estimates = _dy_d_best_reg_estimates(x_values, y_values, deriv_points) - lower = estimates[:n_eval] - two_fx = 2.0 * estimates[n_eval : 2 * n_eval] - upper = estimates[2 * n_eval :] - firsts.append((upper - lower) / distance) - seconds.append((upper - two_fx + lower) / distance**2) - else: - firsts.append(np.full(n_eval, np.nan)) - seconds.append(np.full(n_eval, np.nan)) - mixed_eval_points = eval_mat - - if mixed: - if vector_branch: - tuple_eval = np.asarray(eval_values, dtype=np.float64).reshape(-1) - if tuple_eval.size != 2: - raise ValueError("Mixed Derivatives are only for 2 IV") - mixed_eval_points = tuple_eval.reshape(1, 2) - assert mixed_eval_points is not None - if mixed_eval_points.shape[1] != 2: - raise ValueError("Mixed Derivatives are only for 2 IV") - step_1 = _quantile_step(0, float(band)) - step_2 = _quantile_step(1, float(band)) - if np.isfinite(step_1) and np.isfinite(step_2) and step_1 != 0.0 and step_2 != 0.0: - corners = np.vstack( - [ - np.array( - [ - [p[0] + step_1, p[1] + step_2], - [p[0] - step_1, p[1] + step_2], - [p[0] + step_1, p[1] - step_2], - [p[0] - step_1, p[1] - step_2], - ] - ) - for p in mixed_eval_points - ] + n_eval = eval_mat.shape[0] + for band in h_s: + steps = np.array( + [_dydx_step(column, eval_mat[i, wrt_index], float(band)) for i in range(n_eval)] + ) + finite = np.isfinite(steps) & (steps != 0.0) + lower = eval_mat.copy() + upper = eval_mat.copy() + lower[:, wrt_index] = eval_mat[:, wrt_index] - steps + upper[:, wrt_index] = eval_mat[:, wrt_index] + steps + if messages: + print( + "Currently generating NNS.reg finite difference estimates...bandwidth\r" ) - mixed_estimates = _dy_d_best_reg_estimates(x_values, y_values, corners) - z = mixed_estimates.reshape(-1, 4) - mixeds.append((z[:, 0] + z[:, 3] - z[:, 1] - z[:, 2]) / (4.0 * step_1 * step_2)) - else: - mixeds.append(np.full(mixed_eval_points.shape[0], np.nan)) + estimates = _dy_d_best_reg_estimates( + x_values, y_values, np.vstack((lower, eval_mat, upper)) + ) + lo = estimates[:n_eval] + mid = estimates[n_eval : 2 * n_eval] + up = estimates[2 * n_eval :] + with np.errstate(invalid="ignore", divide="ignore"): + first = (up - lo) / (2.0 * steps) + second = (up - 2.0 * mid + lo) / (steps**2) + first[~finite] = np.nan + second[~finite] = np.nan + firsts.append(first) + seconds.append(second) + mixed_eval = eval_mat def _row_nanmean(bands: list[NDArray[np.float64]]) -> NDArray[np.float64]: matrix = np.column_stack([np.asarray(b, dtype=np.float64).reshape(-1) for b in bands]) @@ -531,8 +517,39 @@ def _row_nanmean(bands: list[NDArray[np.float64]]) -> NDArray[np.float64]: return np.asarray(np.nanmean(matrix, axis=1), dtype=np.float64) output = {"First": _row_nanmean(firsts), "Second": _row_nanmean(seconds)} + if mixed: + if vector_branch: + tuple_eval = np.asarray(eval_values, dtype=np.float64).reshape(-1) + if tuple_eval.size != 2: + raise ValueError("Mixed Derivatives are only for 2 IV") + mixed_points = tuple_eval.reshape(1, 2) + else: + assert mixed_eval is not None + if mixed_eval.shape[1] != 2: + raise ValueError("Mixed Derivatives are only for 2 IV") + mixed_points = mixed_eval + for band in h_s: + band_vals: list[float] = [] + for point in mixed_points: + s1 = _dydx_step(x_values[:, 0], point[0], float(band)) + s2 = _dydx_step(x_values[:, 1], point[1], float(band)) + if not (np.isfinite(s1) and np.isfinite(s2) and s1 != 0.0 and s2 != 0.0): + band_vals.append(np.nan) + continue + corners = np.array( + [ + [point[0] + s1, point[1] + s2], + [point[0] - s1, point[1] + s2], + [point[0] + s1, point[1] - s2], + [point[0] - s1, point[1] - s2], + ] + ) + z = _dy_d_best_reg_estimates(x_values, y_values, corners) + band_vals.append((z[0] + z[3] - z[1] - z[2]) / (4.0 * s1 * s2)) + mixeds.append(np.asarray(band_vals, dtype=np.float64)) output["Mixed"] = _row_nanmean(mixeds) + if messages: print("\r") return output diff --git a/tests/invariants/test_diff.py b/tests/invariants/test_diff.py index de3faea..81a589f 100644 --- a/tests/invariants/test_diff.py +++ b/tests/invariants/test_diff.py @@ -188,7 +188,7 @@ def test_dy_d_best_matches_pinned_values() -> None: np.testing.assert_allclose( np.ravel(result["First"]), - np.array([0.974410, 1.085425, 0.978869]), + np.array([0.754963, 0.792164, 0.699243]), rtol=0.0, atol=1e-4, )