From f09d67b136a6eebd8c0dcd63a5a26d43313bf01c Mon Sep 17 00:00:00 2001 From: Thibault Vatter Date: Thu, 24 Sep 2026 15:01:10 +0200 Subject: [PATCH] fix: keep degree-1 and degree-2 fits finite and positive on tied, gappy data Three defects, found on one column of a density-estimation benchmark (1023 observations, 75 distinct values, a clump of ties at a cap, and a gap of about ten bandwidths below it). A point in the gap got a density of exactly zero, so any held-out observation there had log-likelihood -inf. Bandwidth selection. `ll_ibias2` averages the bias functional over all bins, weighted by their counts. At degree 1 and 2 the functional divides by the pilot density, which is zero up to FFT round-off in empty bins far from the data, so `0 * inf` or `0 * NaN` made the sum NaN on some builds and finite on others. `select_bandwidth` then took the normal-reference fallback (0.710) on the first and the plug-in value (0.2245) on the second. Only occupied bins contribute now, so every build selects the plug-in value. The degree-2 fit. The local quadratic of Hjort and Jones exists only where `1 + h^2 D > 0`, with `D` the second derivative of the pilot log-density. Next to a clump of ties the pilot log-density is as curved as the kernel's own and the condition fails. `sqrt` of a negative number then gave NaN, which `fit_lp` wrote as a density of zero. Where the local quadratic does not exist or is not finite, the local linear fit now stands in, with its own influence value for the degrees of freedom. The degree-1 fit. The local linear estimate is `f * exp(-h^2 b^2 / 2)` with `b = f' / f`, but `fit_lp` used `h` for `h^2`. The estimate was therefore not scale equivariant: on the same standard-normal sample rescaled by 0.1, 1 and 10, the density at 1.5 read 0.066, 0.148 and 0.155. It reads 0.1538 at every scale now. Every degree-1 fit changes, bounded ones included: those are fitted on the transformed scale, where the same `h` for `h^2` applies. Positivity. The binned kernel is truncated at four bandwidths, so a point in a wide enough gap still had density zero at any degree. The grid values between the smallest and largest observation are floored at machine epsilon times the largest value. The interpolating spline preserves positivity, so the density is positive everywhere in that range. Tests: scale equivariance of unbounded fits at every degree (the existing tests covered only bounded fits, where the transform hid the degree-1 error), positivity across a synthetic gap of twenty bandwidths, and the benchmark column as test data, pinning the plug-in bandwidth and positivity at every degree. All three fail before this change. Co-Authored-By: Claude Opus 5.5 (1M context) --- include/kde1d/dpik.hpp | 10 +- include/kde1d/kde1d.hpp | 61 +++-- test/data/tied_gap_sample.hpp | 358 +++++++++++++++++++++++++++++ test/test_numerical_invariants.cpp | 75 ++++++ 4 files changed, 485 insertions(+), 19 deletions(-) create mode 100644 test/data/tied_gap_sample.hpp diff --git a/include/kde1d/dpik.hpp b/include/kde1d/dpik.hpp index 8dba408..60b7b07 100644 --- a/include/kde1d/dpik.hpp +++ b/include/kde1d/dpik.hpp @@ -134,7 +134,15 @@ PluginBandwidthSelector::ll_ibias2(size_t degree) } else { throw std::invalid_argument("degree must be one of {0, 1, 2}."); } - return bin_counts_.cwiseProduct(arg).sum() / bin_counts_.sum(); + // The functional is an average over the observations, so only occupied bins + // contribute. An empty bin far from the data has a pilot density of zero up + // to FFT round-off, where `arg` is inf or NaN and `0 * arg` is NaN. + double ibias2 = 0.0; + for (Eigen::Index k = 0; k < arg.size(); ++k) { + if (bin_counts_(k) > 0.0) + ibias2 += bin_counts_(k) * arg(k); + } + return ibias2 / bin_counts_.sum(); } //! computes the integrated squared variance (without bandwidth and n terms). diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index 256e673..ef90fb5 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -288,7 +288,8 @@ class Kde1d const double& f2, const double& bandwidth, const double& s, - const double& weight); + const double& weight, + const size_t& degree); Eigen::VectorXd boundary_transform(const Eigen::VectorXd& x, bool inverse = false); Eigen::VectorXd boundary_correct(const Eigen::VectorXd& x, @@ -629,6 +630,8 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) boundary_offset_ = 1e-5 * boundary_scale_; } + const double data_min = xx.minCoeff(); + const double data_max = xx.maxCoeff(); xx = boundary_transform(xx); // bandwidth selection (from the request, so refitting re-selects) @@ -644,6 +647,16 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) // order grid points from left to right grid_points = finalize_grid(grid_points); + // Keep the estimate positive between the smallest and largest observation. + // The binned kernel is truncated, so a point in a wide enough gap between + // observations would otherwise get a density of exactly zero. + const double value_floor = + std::numeric_limits::epsilon() * values.maxCoeff(); + for (Eigen::Index k = 0; k < values.size(); ++k) { + if (grid_points(k) >= data_min && grid_points(k) <= data_max) + values(k) = std::max(values(k), value_floor); + } + Eigen::VectorXd influences = fitted.col(1); if (std::isnan(effective_xmin) && !std::isnan(effective_xmax)) influences.reverseInPlace(); @@ -965,27 +978,38 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, // degree > 0 f1 = kde_fft.kde_drv(1); - Eigen::VectorXd S = Eigen::VectorXd::Constant(f0.size(), bandwidth_); - Eigen::VectorXd b = f1.cwiseQuotient(f0); - if (degree_ == 2) { + if (degree_ == 2) f2 = kde_fft.kde_drv(2); - // D/R is notation from Hjort and Jones' AoS paper - Eigen::VectorXd D = f2.cwiseQuotient(f0) - b.cwiseProduct(b); - Eigen::VectorXd R = 1 / (1.0 + bandwidth_ * bandwidth_ * D.array()).sqrt(); - // this is our notation - S = (R / bandwidth_).array().pow(2); - b *= bandwidth_ * bandwidth_; - res.col(0) = bandwidth_ * S.cwiseSqrt().cwiseProduct(res.col(0)); - } - res.col(0) = res.col(0).array() * (-0.5 * b.array().pow(2) * S.array()).exp(); + const double B = bandwidth_ * bandwidth_; for (size_t k = 0; k < m; k++) { if (!std::isfinite(f0(k)) || f0(k) <= density_floor) { res.row(k).setZero(); continue; } - res(k, 1) = - calculate_infl(x.size(), f0(k), f1(k), f2(k), bandwidth_, S(k), wbin(k)); + const double b = f1(k) / f0(k); + // local linear fit (S = 1 / h^2) + double fit = f0(k) * std::exp(-0.5 * B * b * b); + double S = 1.0 / B; + size_t degree = 1; + if (degree_ == 2) { + // D/R is notation from Hjort and Jones' AoS paper: the local quadratic + // exists only where 1 + h^2 D > 0, which fails where the pilot's + // log-density is as curved as the kernel's own, as next to a clump of + // ties. There the local linear fit stands in. + const double denom = 1.0 + B * (f2(k) / f0(k) - b * b); + const double S2 = 1.0 / (B * denom); + const double fit2 = + f0(k) * std::sqrt(B * S2) * std::exp(-0.5 * B * B * b * b * S2); + if (denom > 0.0 && std::isfinite(fit2)) { + fit = fit2; + S = S2; + degree = 2; + } + } + res(k, 0) = fit; + res(k, 1) = calculate_infl( + x.size(), f0(k), f1(k), f2(k), bandwidth_, S, wbin(k), degree); if (!std::isfinite(res(k, 0)) || res(k, 0) < 0.0) res.row(k).setZero(); } @@ -1002,13 +1026,14 @@ Kde1d::calculate_infl(const size_t& n, const double& f2, const double& bandwidth, const double& s, - const double& weight) + const double& weight, + const size_t& degree) { double M_inverse00; double B = bandwidth * bandwidth; - if (degree_ == 0) { + if (degree == 0) { M_inverse00 = 1 / f0; - } else if (degree_ == 1) { + } else if (degree == 1) { Eigen::Matrix2d M; M(0, 0) = f0; M(0, 1) = B * f1; diff --git a/test/data/tied_gap_sample.hpp b/test/data/tied_gap_sample.hpp new file mode 100644 index 0000000..7f45045 --- /dev/null +++ b/test/data/tied_gap_sample.hpp @@ -0,0 +1,358 @@ +#pragma once + +// One column of a density-estimation fold: residual sugar from the UCI Wine +// Quality (red) data, standardized and capped at 6.696. It has 1023 +// observations but only 75 distinct values, a clump of ties at the cap, and no +// observation between 4.495 and 6.696. + +#include + +inline Eigen::VectorXd +tied_gap_sample() +{ + static const double values[] = { + -0.59307275376750412, -0.45368008558748252, 4.4947596348032866, + -0.45368008558748252, -0.3839837514974716, -0.52337641967749327, + -0.52337641967749327, -0.24459108331744986, 0.87055026212272313, + -0.3839837514974716, 6.6964668287065958, -0.3839837514974716, + -0.45368008558748252, -0.10519841513742838, 0.17358692122261488, + -0.66276908785751487, -0.91437285392245404, 0.24328325531262576, + 1.1493355984827667, -0.24459108331744986, -0.31428741740746075, + -0.17489474922743928, -0.24459108331744986, 0.17358692122261488, + 1.4629691018878155, -0.45368008558748252, 0.034194253042593402, + -0.52337641967749327, -0.3839837514974716, -0.10519841513742838, + -0.24459108331744986, -0.31428741740746075, -0.035502081047417484, + -0.31428741740746075, -0.52337641967749327, -0.66276908785751487, + -0.66276908785751487, 0.034194253042593402, -0.59307275376750412, + 0.87055026212272313, 0.31297958940263665, 0.17358692122261488, + -0.035502081047417484, -0.59307275376750412, 0.24328325531262576, + 0.17358692122261488, -0.45368008558748252, -0.10519841513742838, + -0.52337641967749327, -0.24459108331744986, -0.3839837514974716, + 0.1038905871326043, -0.52337641967749327, -0.52337641967749327, + 0.31297958940263665, 0.24328325531262576, -0.24459108331744986, + -0.3839837514974716, -0.3839837514974716, -0.66276908785751487, + -0.24459108331744986, -0.59307275376750412, -0.035502081047417484, + 2.0553879416529073, -0.59307275376750412, -0.00065391400241219698, + -0.17489474922743928, -0.45368008558748252, 0.24328325531262576, + -0.52337641967749327, -0.45368008558748252, -0.31428741740746075, + -0.59307275376750412, -0.66276908785751487, -0.10519841513742838, + -0.73246542194752573, 1.8114507723378697, -0.31428741740746075, + -0.3839837514974716, 1.009942930302745, -0.59307275376750412, + -0.10519841513742838, -0.31428741740746075, 0.1038905871326043, + -0.3839837514974716, -0.17489474922743928, 0.9402465962127341, + -0.24459108331744986, -0.59307275376750412, -0.17489474922743928, + 0.59176492576267992, -0.10519841513742838, -0.31428741740746075, + -0.10519841513742838, -0.45368008558748252, 0.24328325531262576, + -0.73246542194752573, -0.52337641967749327, 0.59176492576267992, + -0.52337641967749327, 0.31297958940263665, -0.31428741740746075, + -0.52337641967749327, 2.5432622802829834, 0.1038905871326043, + -0.31428741740746075, 1.1493355984827667, 0.034194253042593402, + 0.24328325531262576, 0.17358692122261488, -0.59307275376750412, + -0.035502081047417484, 0.1038905871326043, -0.035502081047417484, + -0.3839837514974716, -0.17489474922743928, -0.24459108331744986, + -0.59307275376750412, -0.66276908785751487, -0.035502081047417484, + 0.17358692122261488, -0.73246542194752573, -0.31428741740746075, + 0.17358692122261488, -0.3839837514974716, -0.17489474922743928, + 0.034194253042593402, -0.31428741740746075, 1.009942930302745, + -0.17489474922743928, -0.24459108331744986, -0.17489474922743928, + -0.59307275376750412, -0.8021617560375367, -0.52337641967749327, + -0.31428741740746075, 1.4281209348428097, 1.1493355984827667, + -0.035502081047417484, -0.31428741740746075, 0.17358692122261488, + -0.8021617560375367, -0.45368008558748252, -0.66276908785751487, + -0.24459108331744986, -0.24459108331744986, -0.17489474922743928, + 0.52206859167266906, 0.24328325531262576, -0.035502081047417484, + 0.52206859167266906, 2.8917439507330376, -0.17489474922743928, + -0.3839837514974716, 0.034194253042593402, 0.17358692122261488, + -0.73246542194752573, 0.17358692122261488, -0.17489474922743928, + -0.17489474922743928, -0.73246542194752573, -0.17489474922743928, + -0.3839837514974716, -0.55822458672249875, -0.3839837514974716, + -0.3839837514974716, -0.8021617560375367, 0.034194253042593402, + 2.2644769439229395, -0.52337641967749327, 1.009942930302745, + -0.52337641967749327, -0.45368008558748252, -0.31428741740746075, + -0.59307275376750412, -0.10519841513742838, -0.8021617560375367, + -0.3839837514974716, -0.035502081047417484, 0.59176492576267992, + -0.035502081047417484, -0.45368008558748252, 1.184183765527772, + -0.035502081047417484, -0.3839837514974716, -0.10519841513742838, + -0.035502081047417484, 0.45237225758265842, -0.3839837514974716, + 0.17358692122261488, -0.8021617560375367, 0.38267592349264756, + -0.31428741740746075, -0.45368008558748252, -0.73246542194752573, + -0.3839837514974716, -0.24459108331744986, -0.10519841513742838, + -0.24459108331744986, -0.8021617560375367, -0.52337641967749327, + -0.31428741740746075, 0.034194253042593402, -0.31428741740746075, + -0.45368008558748252, -0.24459108331744986, -0.17489474922743928, + -0.8021617560375367, -0.17489474922743928, 0.31297958940263665, + -0.31428741740746075, -0.31428741740746075, 0.034194253042593402, + 2.4735659461929722, 2.3341732780129507, -0.52337641967749327, + -0.31428741740746075, 0.73115759394270163, 0.17358692122261488, + 0.24328325531262576, -0.24459108331744986, -0.035502081047417484, + -0.24459108331744986, -0.3839837514974716, -0.52337641967749327, + 0.1038905871326043, -0.52337641967749327, 0.034194253042593402, + -0.52337641967749327, 1.9856916075628968, -0.91437285392245404, + -0.31428741740746075, 0.034194253042593402, 3.8674926279931885, + -0.27943925036245543, -0.87185809012754745, 0.034194253042593402, + -0.035502081047417484, -0.3839837514974716, -0.10519841513742838, + -0.45368008558748252, 0.034194253042593402, 0.17358692122261488, + 1.009942930302745, 1.2190319325727772, -0.45368008558748252, + -0.3839837514974716, 1.2190319325727772, -0.17489474922743928, + -0.59307275376750412, -0.24459108331744986, -0.52337641967749327, + 0.1038905871326043, -0.31428741740746075, 0.59176492576267992, + 0.034194253042593402, -0.035502081047417484, 4.4250633007132754, + -0.3839837514974716, -0.3839837514974716, -0.10519841513742838, + 0.034194253042593402, -0.17489474922743928, 2.125084275742918, + -0.17489474922743928, 0.034194253042593402, -0.17489474922743928, + -0.52337641967749327, -0.3839837514974716, -0.91437285392245404, + 4.0068852961732109, 2.125084275742918, 1.2887282666627884, + -0.035502081047417484, 2.0553879416529073, -0.035502081047417484, + -0.3839837514974716, -0.8021617560375367, -0.10519841513742838, + -0.31428741740746075, -0.10519841513742838, -0.17489474922743928, + -0.73246542194752573, -0.3839837514974716, 0.31297958940263665, + -0.45368008558748252, 0.45237225758265842, -0.3839837514974716, + -0.45368008558748252, -0.59307275376750412, -0.45368008558748252, + -0.035502081047417484, -0.66276908785751487, -0.24459108331744986, + -0.24459108331744986, -0.035502081047417484, -0.10519841513742838, + 0.24328325531262576, 1.2190319325727772, -0.3839837514974716, + 0.76600576098770701, -0.10519841513742838, -0.035502081047417484, + -0.17489474922743928, -0.24459108331744986, 0.1038905871326043, + -0.59307275376750412, -0.3839837514974716, -0.3839837514974716, + -0.24459108331744986, -0.17489474922743928, -0.66276908785751487, + 4.0068852961732109, -0.66276908785751487, 0.45237225758265842, + -0.52337641967749327, -0.17489474922743928, -0.10519841513742838, + -0.17489474922743928, -0.10519841513742838, 0.17358692122261488, + -0.91437285392245404, -0.45368008558748252, -0.52337641967749327, + -0.3839837514974716, 0.034194253042593402, -0.31428741740746075, + -0.87185809012754745, -0.24459108331744986, -0.52337641967749327, + -0.3839837514974716, -0.17489474922743928, -0.52337641967749327, + -0.17489474922743928, -0.17489474922743928, -0.66276908785751487, + -0.66276908785751487, -0.17489474922743928, -0.52337641967749327, + -0.10519841513742838, -0.31428741740746075, -0.17489474922743928, + 0.9402465962127341, -0.8021617560375367, 0.034194253042593402, + -0.52337641967749327, 1.7766026052928638, 0.45237225758265842, + -0.035502081047417484, -0.45368008558748252, -0.10519841513742838, + -0.24459108331744986, 0.1038905871326043, 3.7280999598131674, + -0.20974291627244457, 0.9402465962127341, -0.45368008558748252, + -0.59307275376750412, -0.31428741740746075, -0.3839837514974716, + -0.52337641967749327, -0.31428741740746075, -0.45368008558748252, + -0.24459108331744986, -0.31428741740746075, -0.24459108331744986, + -0.31428741740746075, -0.3839837514974716, 2.7871994495980208, + -0.31428741740746075, -0.24459108331744986, -0.3839837514974716, + -0.52337641967749327, -0.31428741740746075, -0.10519841513742838, + -0.59307275376750412, -0.035502081047417484, 0.17358692122261488, + 0.17358692122261488, -0.45368008558748252, -0.31428741740746075, + 0.034194253042593402, -0.3839837514974716, -0.17489474922743928, + -0.45368008558748252, -0.035502081047417484, -0.035502081047417484, + -0.31428741740746075, -0.24459108331744986, -0.45368008558748252, + -0.17489474922743928, -0.3839837514974716, -0.66276908785751487, + 0.034194253042593402, -0.59307275376750412, 0.034194253042593402, + -0.24459108331744986, -0.31428741740746075, 1.1493355984827667, + 0.83570209507771787, -0.31428741740746075, 0.52206859167266906, + -0.035502081047417484, 0.8008539280327126, 0.45237225758265842, + -0.52337641967749327, 1.2190319325727772, 0.17358692122261488, + -0.24459108331744986, -0.17489474922743928, -0.59307275376750412, + 6.6964668287065958, -0.17489474922743928, -0.45368008558748252, + -0.31428741740746075, -0.3839837514974716, -0.3839837514974716, + -0.035502081047417484, 0.034194253042593402, -0.45368008558748252, + -0.17489474922743928, -0.24459108331744986, -0.3839837514974716, + -0.52337641967749327, -0.10519841513742838, 0.034194253042593402, + -0.035502081047417484, -0.31428741740746075, -0.10519841513742838, + 1.3584246007527991, -0.66276908785751487, -0.24459108331744986, + -0.3839837514974716, -0.17489474922743928, 1.846298939382875, + -0.45368008558748252, 0.87055026212272313, 0.73115759394270163, + -0.17489474922743928, -0.31428741740746075, -0.24459108331744986, + -0.31428741740746075, -0.3839837514974716, -0.17489474922743928, + -0.52337641967749327, -0.17489474922743928, -0.3839837514974716, + 1.2190319325727772, 4.0068852961732109, -0.24459108331744986, + 0.1038905871326043, -0.73246542194752573, 0.38267592349264756, + -0.31428741740746075, -0.24459108331744986, -0.3839837514974716, + -0.10519841513742838, 0.034194253042593402, -0.24459108331744986, + 0.45237225758265842, 0.17358692122261488, -0.17489474922743928, + -0.10519841513742838, 0.034194253042593402, -0.59307275376750412, + 0.24328325531262576, -0.45368008558748252, -0.24459108331744986, + -0.45368008558748252, -0.45368008558748252, 0.73115759394270163, + 0.034194253042593402, 0.59176492576267992, -0.24459108331744986, + -0.59307275376750412, -0.24459108331744986, 1.2190319325727772, + -0.3839837514974716, -0.31428741740746075, -0.17489474922743928, + 0.31297958940263665, -0.73246542194752573, 2.125084275742918, + -0.52337641967749327, 0.1038905871326043, -0.3839837514974716, + -0.8021617560375367, -0.3839837514974716, -0.3839837514974716, + 1.009942930302745, -0.59307275376750412, -0.73246542194752573, + -0.14004658218243365, -0.10519841513742838, -0.66276908785751487, + -0.66276908785751487, -0.035502081047417484, 0.45237225758265842, + -0.66276908785751487, -0.73246542194752573, -0.3839837514974716, + -0.3839837514974716, -0.17489474922743928, -0.66276908785751487, + 0.034194253042593402, -0.17489474922743928, -0.8021617560375367, + -0.10519841513742838, 0.034194253042593402, -0.24459108331744986, + -0.17489474922743928, -0.24459108331744986, -0.035502081047417484, + -0.31428741740746075, 1.3584246007527991, 0.1038905871326043, + 0.034194253042593402, -0.3839837514974716, -0.3839837514974716, + -0.52337641967749327, 0.1038905871326043, 0.034194253042593402, + -0.17489474922743928, -0.59307275376750412, -0.45368008558748252, + -0.66276908785751487, -0.3839837514974716, -0.59307275376750412, + -0.24459108331744986, -0.45368008558748252, -0.52337641967749327, + -0.73246542194752573, -0.66276908785751487, 0.52206859167266906, + -0.17489474922743928, -0.035502081047417484, -0.52337641967749327, + -0.10519841513742838, -0.62792092081250961, -0.31428741740746075, + -0.91437285392245404, 1.009942930302745, -0.035502081047417484, + -0.3839837514974716, -0.59307275376750412, -0.10519841513742838, + 0.1038905871326043, -0.10519841513742838, 1.0796392643927555, + -0.31428741740746075, -0.10519841513742838, -0.3839837514974716, + -0.45368008558748252, -0.31428741740746075, -0.66276908785751487, + -0.10519841513742838, -0.66276908785751487, -0.035502081047417484, + 1.009942930302745, -0.45368008558748252, 0.1038905871326043, + -0.17489474922743928, -0.52337641967749327, -0.31428741740746075, + 0.1038905871326043, -0.52337641967749327, -0.3839837514974716, + 0.24328325531262576, 2.0553879416529073, -0.035502081047417484, + 1.5675136030228314, -0.035502081047417484, -0.10519841513742838, + -0.24459108331744986, 0.73115759394270163, -0.66276908785751487, + -0.59307275376750412, -0.24459108331744986, -0.3839837514974716, + -0.45368008558748252, 0.24328325531262576, 0.1038905871326043, + 0.31297958940263665, 0.24328325531262576, -0.035502081047417484, + -0.035502081047417484, -0.10519841513742838, -0.8021617560375367, + -0.31428741740746075, -0.59307275376750412, -0.45368008558748252, + 0.1038905871326043, -0.17489474922743928, -0.52337641967749327, + -0.24459108331744986, 0.17358692122261488, -0.3839837514974716, + -0.17489474922743928, -0.31428741740746075, 2.6826549484630049, + 0.17358692122261488, -0.035502081047417484, 0.38267592349264756, + -0.52337641967749327, 0.87055026212272313, -0.24459108331744986, + -0.45368008558748252, -0.17489474922743928, -0.10519841513742838, + -0.3839837514974716, -0.59307275376750412, 1.5675136030228314, + -0.3839837514974716, -0.17489474922743928, -0.31428741740746075, + -0.52337641967749327, -0.3839837514974716, -0.66276908785751487, + -0.035502081047417484, -0.91437285392245404, -0.59307275376750412, + 0.87055026212272313, -0.035502081047417484, -0.45368008558748252, + -0.035502081047417484, -0.31428741740746075, -0.17489474922743928, + 0.034194253042593402, -0.59307275376750412, 0.034194253042593402, + -0.8021617560375367, -0.24459108331744986, -0.31428741740746075, + 0.1038905871326043, 1.2887282666627884, 0.1038905871326043, + -0.31428741740746075, -0.31428741740746075, -0.24459108331744986, + -0.3839837514974716, -0.035502081047417484, 0.034194253042593402, + -0.52337641967749327, -0.45368008558748252, 0.17358692122261488, + -0.3839837514974716, 0.24328325531262576, -0.24459108331744986, + -0.45368008558748252, -0.3839837514974716, -0.10519841513742838, + -0.59307275376750412, -0.59307275376750412, 0.034194253042593402, + -0.59307275376750412, -0.59307275376750412, 0.24328325531262576, + 0.034194253042593402, -0.45368008558748252, 0.59176492576267992, + -0.10519841513742838, -0.3839837514974716, 0.17358692122261488, + 2.5432622802829834, -0.8021617560375367, 0.17358692122261488, + 4.3553669666232651, -0.45368008558748252, 0.9402465962127341, + -0.45368008558748252, 0.52206859167266906, 1.7766026052928638, + -0.24459108331744986, -0.3839837514974716, -0.52337641967749327, + -0.52337641967749327, 0.034194253042593402, 2.125084275742918, + -0.8021617560375367, 6.6964668287065958, -0.52337641967749327, + -0.24459108331744986, -0.24459108331744986, -0.45368008558748252, + -0.3839837514974716, 0.034194253042593402, 0.1038905871326043, + -0.59307275376750412, -0.17489474922743928, -0.52337641967749327, + -0.45368008558748252, -0.59307275376750412, -0.8021617560375367, + 0.24328325531262576, -0.59307275376750412, -0.55822458672249875, + 0.52206859167266906, -0.45368008558748252, -0.45368008558748252, + -0.24459108331744986, -0.45368008558748252, -0.3839837514974716, + -0.10519841513742838, -0.66276908785751487, -0.10519841513742838, + 0.45237225758265842, -0.035502081047417484, -0.035502081047417484, + -0.24459108331744986, -0.24459108331744986, -0.10519841513742838, + -0.3839837514974716, -0.24459108331744986, 0.034194253042593402, + -0.3839837514974716, -0.24459108331744986, 1.7766026052928638, + -0.035502081047417484, -0.24459108331744986, -0.10519841513742838, + -0.3839837514974716, 0.034194253042593402, -0.24459108331744986, + -0.24459108331744986, -0.24459108331744986, -0.31428741740746075, + -0.17489474922743928, -0.10519841513742838, -0.31428741740746075, + -0.3839837514974716, -0.45368008558748252, -0.10519841513742838, + -0.66276908785751487, -0.52337641967749327, -0.59307275376750412, + -0.24459108331744986, 0.45237225758265842, 0.17358692122261488, + -0.66276908785751487, -0.035502081047417484, -0.45368008558748252, + 0.17358692122261488, -0.45368008558748252, -0.45368008558748252, + 0.45237225758265842, -0.73246542194752573, 2.4735659461929722, + 1.0796392643927555, -0.73246542194752573, 0.1038905871326043, + -0.31428741740746075, -0.52337641967749327, -0.52337641967749327, + 0.45237225758265842, -0.31428741740746075, -0.17489474922743928, + -0.3839837514974716, -0.17489474922743928, 0.66146125985269077, + -0.45368008558748252, 0.45237225758265842, 3.6584036257231562, + -0.3839837514974716, -0.73246542194752573, -0.66276908785751487, + -0.035502081047417484, -0.10519841513742838, -0.31428741740746075, + -0.17489474922743928, -0.66276908785751487, -0.10519841513742838, + -0.24459108331744986, -0.45368008558748252, 0.17358692122261488, + -0.8021617560375367, -0.59307275376750412, -0.3839837514974716, + -0.52337641967749327, 0.31297958940263665, -0.17489474922743928, + -0.035502081047417484, 0.24328325531262576, -0.31428741740746075, + -0.52337641967749327, -0.45368008558748252, -0.24459108331744986, + -0.8021617560375367, 0.17358692122261488, 0.52206859167266906, + -0.52337641967749327, -0.52337641967749327, 0.034194253042593402, + -0.62792092081250961, 0.17358692122261488, -0.24459108331744986, + 0.31297958940263665, -0.035502081047417484, 6.6964668287065958, + 1.4281209348428097, -0.45368008558748252, -0.66276908785751487, + -0.45368008558748252, 0.034194253042593402, -0.3839837514974716, + -0.59307275376750412, -0.035502081047417484, -0.45368008558748252, + -0.73246542194752573, 0.034194253042593402, 0.034194253042593402, + -0.24459108331744986, 3.1008329530030698, -0.52337641967749327, + -0.17489474922743928, -0.10519841513742838, -0.10519841513742838, + 0.59176492576267992, 0.59176492576267992, -0.17489474922743928, + 1.4281209348428097, -0.87185809012754745, -0.45368008558748252, + -0.31428741740746075, -0.035502081047417484, -0.52337641967749327, + -0.73246542194752573, -0.52337641967749327, -0.45368008558748252, + 0.034194253042593402, 0.17358692122261488, -0.73246542194752573, + -0.24459108331744986, 1.2887282666627884, 0.17358692122261488, + -0.52337641967749327, -0.10519841513742838, -0.52337641967749327, + -0.035502081047417484, -0.24459108331744986, 0.1038905871326043, + -0.3839837514974716, -0.52337641967749327, -0.17489474922743928, + 2.0553879416529073, -0.31428741740746075, -0.31428741740746075, + 1.2887282666627884, 2.4038696121029615, -0.24459108331744986, + -0.3839837514974716, 0.17358692122261488, -0.035502081047417484, + -0.45368008558748252, 0.31297958940263665, -0.035502081047417484, + 0.1038905871326043, -0.31428741740746075, 3.7280999598131674, + -0.31428741740746075, -0.59307275376750412, -0.31428741740746075, + -0.10519841513742838, 0.24328325531262576, -0.035502081047417484, + -0.17489474922743928, -0.73246542194752573, -0.17489474922743928, + -0.52337641967749327, -0.17489474922743928, -0.24459108331744986, + 0.17358692122261488, -0.24459108331744986, 0.17358692122261488, + -0.34913558445246634, 0.17358692122261488, -0.035502081047417484, + -0.59307275376750412, -0.24459108331744986, -0.3839837514974716, + -0.17489474922743928, -0.24459108331744986, -0.73246542194752573, + 0.034194253042593402, -0.52337641967749327, -0.17489474922743928, + -0.10519841513742838, -0.3839837514974716, 2.0553879416529073, + -0.45368008558748252, -0.52337641967749327, 0.31297958940263665, + 0.17358692122261488, 0.31297958940263665, -0.3839837514974716, + -0.24459108331744986, -0.24459108331744986, 1.2190319325727772, + -0.17489474922743928, -0.24459108331744986, 1.009942930302745, + 2.6129586143729937, -0.17489474922743928, -0.52337641967749327, + 2.3341732780129507, 0.73115759394270163, -0.24459108331744986, + -0.10519841513742838, 0.034194253042593402, -0.31428741740746075, + 2.0553879416529073, 3.3099219552731021, -0.035502081047417484, + -0.52337641967749327, -0.45368008558748252, -0.17489474922743928, + -0.52337641967749327, -0.31428741740746075, 0.17358692122261488, + -0.035502081047417484, 0.034194253042593402, -0.8021617560375367, + 2.8220476166430264, 0.59176492576267992, -0.035502081047417484, + -0.31428741740746075, 0.59176492576267992, -0.45368008558748252, + -0.45368008558748252, -0.45368008558748252, -0.45368008558748252, + -0.87185809012754745, -0.45368008558748252, -0.10519841513742838, + -0.31428741740746075, -0.10519841513742838, 1.4629691018878155, + -0.035502081047417484, -0.31428741740746075, -0.59307275376750412, + -0.45368008558748252, -0.66276908785751487, -0.3839837514974716, + -0.10519841513742838, -0.035502081047417484, -0.52337641967749327, + -0.45368008558748252, -0.24459108331744986, -0.10519841513742838, + -0.17489474922743928, -0.10519841513742838, -0.31428741740746075, + -0.3839837514974716, -0.45368008558748252, -0.31428741740746075, + -0.59307275376750412, -0.45368008558748252, -0.52337641967749327, + 2.125084275742918, 0.034194253042593402, -0.24459108331744986, + -0.31428741740746075, 2.4038696121029615, -0.10519841513742838, + -0.52337641967749327, -0.8021617560375367, -0.59307275376750412, + -0.45368008558748252, 0.20843508826762047, -0.66276908785751487, + -0.66276908785751487, -0.17489474922743928, 2.125084275742918, + -0.73246542194752573, -0.52337641967749327, -0.3839837514974716, + -0.52337641967749327, -0.45368008558748252, 0.31297958940263665, + 0.38267592349264756, -0.10519841513742838, -0.59307275376750412, + 0.034194253042593402, -0.52337641967749327, 4.2159742984432427, + 2.2644769439229395, 1.0796392643927555, -0.3839837514974716, + -0.45368008558748252, -0.52337641967749327, -0.59307275376750412, + -0.3839837514974716, 1.009942930302745, -0.52337641967749327, + -0.24459108331744986, -0.24459108331744986, -0.31428741740746075, + 0.034194253042593402, -0.3839837514974716, -0.3839837514974716, + -0.52337641967749327, -0.31428741740746075, -0.66276908785751487, + 3.8674926279931885, -0.8021617560375367, -0.10519841513742838, + -0.035502081047417484, -0.45368008558748252, -0.24459108331744986, + 0.52206859167266906, 0.034194253042593402, -0.035502081047417484, + -0.3839837514974716, -0.24459108331744986, -0.52337641967749327, + -0.035502081047417484, 0.31297958940263665, -0.24459108331744986, + -0.31428741740746075, -0.10519841513742838, -0.3839837514974716, + -0.24459108331744986, -0.73246542194752573, -0.52337641967749327 + }; + return Eigen::Map(values, + sizeof(values) / sizeof(values[0])); +} diff --git a/test/test_numerical_invariants.cpp b/test/test_numerical_invariants.cpp index a29345f..163a90c 100644 --- a/test/test_numerical_invariants.cpp +++ b/test/test_numerical_invariants.cpp @@ -1,6 +1,7 @@ #define CATCH_CONFIG_MAIN #include "../include/kde1d.hpp" #include "catch.hpp" +#include "data/tied_gap_sample.hpp" namespace { @@ -206,3 +207,77 @@ TEST_CASE("interpolation preserves nonuniform-grid reference values", CHECK(grid.integrate(queries).head(13).isApprox( expected_normalized_integral, 1e-12)); } + +TEST_CASE("unbounded fits are scale equivariant", + "[numerical-invariants][scale]") +{ + Eigen::VectorXd probabilities = + Eigen::VectorXd::LinSpaced(500, 0.5 / 500.0, 1.0 - 0.5 / 500.0); + Eigen::VectorXd observations = kde1d::stats::qnorm(probabilities); + Eigen::VectorXd points = Eigen::VectorXd::LinSpaced(9, -3.0, 3.0); + + for (size_t degree = 0; degree < 3; ++degree) { + kde1d::Kde1d reference(NAN, NAN, "continuous", 1.0, NAN, degree); + reference.fit(observations); + + for (double scale : { 1e-4, 1e4 }) { + INFO("degree=" << degree << ", scale=" << scale); + kde1d::Kde1d scaled(NAN, NAN, "continuous", 1.0, NAN, degree); + scaled.fit(observations * scale); + + CHECK(scaled.get_bandwidth() / scale == + Approx(reference.get_bandwidth()).epsilon(1e-10)); + CHECK((scaled.pdf(points * scale) * scale) + .isApprox(reference.pdf(points), 1e-10)); + } + } +} + +TEST_CASE("the density is positive between the observations", + "[numerical-invariants][gap]") +{ + // Two clumps twenty bandwidths apart: the binned kernel is truncated well + // before the middle of the gap. + Eigen::VectorXd observations(200); + observations.head(100) = Eigen::VectorXd::LinSpaced(100, -0.05, 0.05); + observations.tail(100) = Eigen::VectorXd::LinSpaced(100, 9.95, 10.05); + Eigen::VectorXd points = Eigen::VectorXd::LinSpaced(401, -0.05, 10.05); + + for (size_t degree = 0; degree < 3; ++degree) { + INFO("degree=" << degree); + kde1d::Kde1d fit(NAN, NAN, "continuous", 1.0, 0.5, degree); + fit.fit(observations); + Eigen::VectorXd density = fit.pdf(points); + REQUIRE(density.array().isFinite().all()); + CHECK(density.minCoeff() > 0.0); + CHECK(std::isfinite(fit.get_loglik())); + } +} + +TEST_CASE("tied data with a gap select the plug-in bandwidth", + "[numerical-invariants][bandwidth][gap]") +{ + // The pilot density is zero up to FFT round-off in the empty bins, which + // must not reach the bias functional: the plug-in value is selected, not the + // normal-reference fallback (0.70996). + Eigen::VectorXd observations = tied_gap_sample(); + kde1d::Kde1d fit; + fit.fit(observations); + CHECK(fit.get_bandwidth() == Approx(0.2245107160846).epsilon(1e-8)); + + // Next to the clump of ties at the cap the local quadratic degenerates, and + // a point in the gap is farther than the kernel reaches from any observation. + Eigen::VectorXd points = Eigen::VectorXd::LinSpaced( + 2001, observations.minCoeff(), observations.maxCoeff()); + for (size_t degree = 0; degree < 3; ++degree) { + INFO("degree=" << degree); + kde1d::Kde1d refit(NAN, NAN, "continuous", 1.0, NAN, degree); + refit.fit(observations); + Eigen::VectorXd density = refit.pdf(points); + REQUIRE(density.array().isFinite().all()); + CHECK(density.minCoeff() > 0.0); + CHECK(refit.pdf(Eigen::VectorXd::Constant(1, 5.6796))(0) > 0.0); + CHECK(std::isfinite(refit.get_loglik())); + CHECK(std::isfinite(refit.get_edf())); + } +}