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())); + } +}