From 8125633dbbb59234eed66e2d4a5405c059978906 Mon Sep 17 00:00:00 2001 From: tnagler Date: Tue, 25 Aug 2026 13:46:29 +0200 Subject: [PATCH 1/3] Handle bounded discrete densities --- docs/overview-continuous.dox | 4 +- docs/overview-discrete.dox | 58 +++++++++ include/kde1d/kde1d.hpp | 236 ++++++++++++++++++++++------------- test/test.cpp | 94 ++++++++++++++ 4 files changed, 306 insertions(+), 86 deletions(-) diff --git a/docs/overview-continuous.dox b/docs/overview-continuous.dox index 4d3dc0e..fa181c0 100644 --- a/docs/overview-continuous.dox +++ b/docs/overview-continuous.dox @@ -41,8 +41,8 @@ it does not report an endpoint-expert bandwidth. `boundary_repair=true` makes each known finite endpoint eligible for repair; it does not force an expert to be used. Repair is considered only for -continuous samples with at least 16 cleaned observations and effective sample -size +continuous or discrete samples with at least 16 cleaned observations and +effective sample size @f[ n_e=\frac{(\sum_i w_i)^2}{\sum_i w_i^2} diff --git a/docs/overview-discrete.dox b/docs/overview-discrete.dox index 11850f1..19211c3 100644 --- a/docs/overview-discrete.dox +++ b/docs/overview-discrete.dox @@ -1,4 +1,62 @@ /** @page overview-discrete Discrete Distributions +Discrete fits model nonnegative integer-valued observations. Both the data and +any finite `xmin` or `xmax` must therefore be nonnegative integers. A finite +bound describes the support of the discrete variable itself; probabilities +outside those integer bounds are zero. + +@section discrete-jitter Jittering and support transforms + +Before fitting, tied observations are spread deterministically within their +unit cells by adding conditionally equidistant noise in @f$(-1/2,1/2)@f$. +Consequently, if the discrete support has a finite lower endpoint @f$a@f$ or +upper endpoint @f$b@f$, the continuous jitter support has the corresponding +endpoint + +@f[ +a_J=a-1/2, +\qquad +b_J=b+1/2. +@f] + +The local-likelihood bulk estimator applies the same support transforms +described in @ref continuous-bulk to these jitter-adjusted endpoints. Thus +`xmin=0`, for example, invokes a one-sided transformation at @f$-1/2@f$; +finite `xmin` and `xmax` invoke a two-sided transformation over the full +jitter cells. The public bounds remain unchanged and integer-valued. + +@section discrete-pmf From jitter density to probability mass + +The fitted interpolation grid represents a continuous jitter density, but the +discrete probability mass function is formed only from its values at integer +levels: + +@f[ +\widehat p(k)= +\frac{\widehat f_J(k)}{\sum_{j\in\mathcal S}\widehat f_J(j)}, +\qquad k\in\mathcal S. +@f] + +Here @f$\mathcal S@f$ is the nonnegative integer support, restricted by any +finite public bounds. Noninteger arguments and integer values outside this +support have probability zero. The CDF and quantile function use the same +normalized masses. + +This distinction also explains why an unconstrained jitter density can give a +good discrete estimate despite having a nonzero continuous integral outside +the integer support: that integral is not assigned to any discrete level. +Only the density ordinates at integer centers enter the PMF. Supplying bounds +uses the known support during fitting and prevents the transformed bulk from +extrapolating through the ends of the boundary cells. + +@section discrete-repair Adaptive boundary repair + +When `boundary_repair=true`, each supplied discrete endpoint is eligible for +the adaptive boundary expert described in @ref continuous-repair and +@ref continuous-fusion. Classification and repair operate on the jittered +observations and the half-unit endpoint. The expert is not forced: an endpoint +whose local behavior appears vanishing, exploding, or ambiguous retains the +transformed bulk estimate. This guards finite-looking boundary cells without +replacing a stable bulk tail based only on the existence of a support bound. */ diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index 1311efa..9f51c75 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -18,10 +18,11 @@ enum class VarType //! Local-polynomial density estimation in one dimension. //! -//! Continuous fits use a transformed local-likelihood estimate as their bulk -//! component. When `boundary_repair` is enabled, finite support endpoints may -//! additionally use a nonnegative local-linear boundary kernel. See -//! @ref overview-continuous for the transforms, +//! Continuous and bounded discrete fits use a transformed local-likelihood +//! estimate as their bulk component. When `boundary_repair` is enabled, +//! finite support endpoints may additionally use a nonnegative local-linear +//! boundary kernel. See +//! @ref overview-continuous and @ref overview-discrete for the transforms, //! endpoint classifier, bandwidths, fusion weights, and EDF approximation. class Kde1d { @@ -200,6 +201,11 @@ class Kde1d void check_inputs(const Eigen::VectorXd& x, const Eigen::VectorXd& weights = Eigen::VectorXd()) const; void check_boundaries(const Eigen::VectorXd& x) const; + void check_discrete_data(const Eigen::VectorXd& x) const; + double get_effective_xmin() const; + double get_effective_xmax() const; + double get_discrete_support_min() const; + double get_discrete_support_max() const; Eigen::VectorXd pdf_continuous(const Eigen::VectorXd& x) const; Eigen::VectorXd cdf_continuous(const Eigen::VectorXd& x) const; Eigen::VectorXd quantile_continuous(const Eigen::VectorXd& x) const; @@ -405,6 +411,7 @@ inline void Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) { check_inputs(x, weights); + check_discrete_data(x); check_boundaries(x); // preprocessing for nans and jittering @@ -412,14 +419,6 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) Eigen::VectorXd w = weights; tools::remove_nans(xx, w); - const bool use_boundary_repair = boundary_repair_ && - type_ == VarType::continuous && - xx.size() >= 16 && - (!std::isnan(xmin_) || !std::isnan(xmax_)); - Eigen::VectorXd boundary_observations; - if (use_boundary_repair) - boundary_observations = xx; - if (w.size() > 0) { w /= w.mean(); } @@ -446,15 +445,24 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) xx = stats::equi_jitter(xx); } - if (type_ != VarType::discrete && - (std::isnan(xmin_) != std::isnan(xmax_))) { + const double effective_xmin = get_effective_xmin(); + const double effective_xmax = get_effective_xmax(); + const bool use_boundary_repair = + boundary_repair_ && + (type_ == VarType::continuous || type_ == VarType::discrete) && + xx.size() >= 16 && + (!std::isnan(effective_xmin) || !std::isnan(effective_xmax)); + const Eigen::VectorXd boundary_observations = + use_boundary_repair ? xx : Eigen::VectorXd(); + + if (std::isnan(effective_xmin) != std::isnan(effective_xmax)) { // Scaling the median boundary distance makes the transform equivariant // under changes of measurement units while remaining robust to outliers. Eigen::VectorXd distances; - if (std::isnan(xmin_)) - distances = xmax_ - xx.array(); + if (std::isnan(effective_xmin)) + distances = effective_xmax - xx.array(); else - distances = xx.array() - xmin_; + distances = xx.array() - effective_xmin; if (w.size() == 0) { boundary_scale_ = stats::median(distances); } else { @@ -484,7 +492,7 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) grid_points = finalize_grid(grid_points); Eigen::VectorXd influences = fitted.col(1); - if (std::isnan(xmin_) && !std::isnan(xmax_)) + if (std::isnan(effective_xmin) && !std::isnan(effective_xmax)) influences.reverseInPlace(); // construct interpolation grid @@ -560,8 +568,8 @@ inline Eigen::VectorXd Kde1d::pdf_discrete(const Eigen::VectorXd& x) const { auto fhat = pdf_continuous(x); - auto lb = std::floor(grid_.get_grid_min()); - auto ub = std::ceil(grid_.get_grid_max()); + const double lb = get_discrete_support_min(); + const double ub = get_discrete_support_max(); Eigen::VectorXd lvs = Eigen::VectorXd::LinSpaced(static_cast(ub - lb + 1), lb, ub); @@ -569,8 +577,13 @@ Kde1d::pdf_discrete(const Eigen::VectorXd& x) const (x.array() >= lb) && (x.array() <= ub) && (x.array() == x.array().round()); fhat = fhat.array() * selected.cast().array(); - // normalize - fhat /= grid_.interpolate(lvs).sum(); + // Normalize density ordinates over the integer support. The continuous + // jitter density between integer levels (and outside an unspecified upper + // support) does not represent probability mass of the discrete model. + const double normalizer = pdf_continuous(lvs).sum(); + if (!(normalizer > 0.0) || !std::isfinite(normalizer)) + throw std::runtime_error("failed to normalize the discrete density"); + fhat /= normalizer; return fhat; } @@ -615,8 +628,8 @@ Kde1d::cdf_continuous(const Eigen::VectorXd& x) const inline Eigen::VectorXd Kde1d::cdf_discrete(const Eigen::VectorXd& x) const { - auto lb = std::floor(grid_.get_grid_min()); - auto ub = std::ceil(grid_.get_grid_max()); + const double lb = get_discrete_support_min(); + const double ub = get_discrete_support_max(); Eigen::VectorXd lvs = Eigen::VectorXd::LinSpaced(static_cast(ub - lb + 1), lb, ub); @@ -677,8 +690,8 @@ Kde1d::quantile_continuous(const Eigen::VectorXd& x) const inline Eigen::VectorXd Kde1d::quantile_discrete(const Eigen::VectorXd& x) const { - auto lb = std::floor(grid_.get_grid_min()); - auto ub = std::ceil(grid_.get_grid_max()); + const double lb = get_discrete_support_min(); + const double ub = get_discrete_support_max(); auto nlevels = static_cast(ub - lb + 1); Eigen::VectorXd lvs = Eigen::VectorXd::LinSpaced(nlevels, lb, ub); @@ -864,26 +877,24 @@ Kde1d::calculate_infl(const size_t& n, inline Eigen::VectorXd Kde1d::boundary_transform(const Eigen::VectorXd& x, bool inverse) { - if (type_ == VarType::discrete) { - return x; // no transform for discrete variables - } - + const double xmin = get_effective_xmin(); + const double xmax = get_effective_xmax(); Eigen::VectorXd x_new = x; if (!inverse) { - if (!std::isnan(xmin_) && !std::isnan(xmax_)) { + if (!std::isnan(xmin) && !std::isnan(xmax)) { // two boundaries -> probit transform - auto rng = xmax_ - xmin_; - x_new = (x.array() - xmin_ + 5e-5 * rng) / (1.0001 * rng); + auto rng = xmax - xmin; + x_new = (x.array() - xmin + 5e-5 * rng) / (1.0001 * rng); x_new = stats::qnorm(x_new); - } else if (!std::isnan(xmin_)) { + } else if (!std::isnan(xmin)) { // left boundary -> regularized fourth-root transform - x_new = 4.0 * (((boundary_offset_ + x.array() - xmin_) / + x_new = 4.0 * (((boundary_offset_ + x.array() - xmin) / boundary_scale_) .pow(0.25) - std::pow(boundary_offset_ / boundary_scale_, 0.25)); - } else if (!std::isnan(xmax_)) { + } else if (!std::isnan(xmax)) { // right boundary -> reflected regularized fourth-root transform - x_new = 4.0 * (((boundary_offset_ + xmax_ - x.array()) / + x_new = 4.0 * (((boundary_offset_ + xmax - x.array()) / boundary_scale_) .pow(0.25) - std::pow(boundary_offset_ / boundary_scale_, 0.25)); @@ -891,20 +902,20 @@ Kde1d::boundary_transform(const Eigen::VectorXd& x, bool inverse) // no boundary -> no transform } } else { - if (!std::isnan(xmin_) && !std::isnan(xmax_)) { + if (!std::isnan(xmin) && !std::isnan(xmax)) { // two boundaries -> probit transform - auto rng = xmax_ - xmin_; - x_new = stats::pnorm(x).array() * 1.0001 * rng + xmin_ - 5e-5 * rng; - } else if (!std::isnan(xmin_)) { + auto rng = xmax - xmin; + x_new = stats::pnorm(x).array() * 1.0001 * rng + xmin - 5e-5 * rng; + } else if (!std::isnan(xmin)) { // left boundary -> inverse regularized fourth-root transform x_new = boundary_scale_ * (x.array() / 4.0 + std::pow(boundary_offset_ / boundary_scale_, 0.25)) .pow(4) + - xmin_ - boundary_offset_; - } else if (!std::isnan(xmax_)) { + xmin - boundary_offset_; + } else if (!std::isnan(xmax)) { // right boundary -> inverse reflected regularized fourth-root transform - x_new = xmax_ + boundary_offset_ - + x_new = xmax + boundary_offset_ - boundary_scale_ * (x.array() / 4.0 + std::pow(boundary_offset_ / boundary_scale_, 0.25)) @@ -925,26 +936,24 @@ Kde1d::boundary_transform(const Eigen::VectorXd& x, bool inverse) inline Eigen::VectorXd Kde1d::boundary_correct(const Eigen::VectorXd& x, const Eigen::VectorXd& fhat) { - if (type_ == VarType::discrete) { - return fhat; // no transform for discrete variables - } - + const double xmin = get_effective_xmin(); + const double xmax = get_effective_xmax(); Eigen::VectorXd corr_term(fhat.size()); - if (!std::isnan(xmin_) && !std::isnan(xmax_)) { + if (!std::isnan(xmin) && !std::isnan(xmax)) { // two boundaries -> probit transform - auto rng = xmax_ - xmin_; - corr_term = (x.array() - xmin_ + 5e-5 * rng) / (xmax_ - xmin_ + 1e-4 * rng); + auto rng = xmax - xmin; + corr_term = (x.array() - xmin + 5e-5 * rng) / (xmax - xmin + 1e-4 * rng); corr_term = stats::dnorm(stats::qnorm(corr_term)); - corr_term *= (xmax_ - xmin_ + 1e-4 * rng); + corr_term *= (xmax - xmin + 1e-4 * rng); corr_term = 1.0 / corr_term.array(); - } else if (!std::isnan(xmin_)) { + } else if (!std::isnan(xmin)) { // left boundary -> fourth-root-transform Jacobian - corr_term = ((boundary_offset_ + x.array() - xmin_) / boundary_scale_) + corr_term = ((boundary_offset_ + x.array() - xmin) / boundary_scale_) .pow(-0.75) / boundary_scale_; - } else if (!std::isnan(xmax_)) { + } else if (!std::isnan(xmax)) { // right boundary -> reflected fourth-root-transform Jacobian - corr_term = ((boundary_offset_ + xmax_ - x.array()) / boundary_scale_) + corr_term = ((boundary_offset_ + xmax - x.array()) / boundary_scale_) .pow(-0.75) / boundary_scale_; } else { @@ -953,7 +962,7 @@ Kde1d::boundary_correct(const Eigen::VectorXd& x, const Eigen::VectorXd& fhat) } Eigen::VectorXd f_corr = fhat.cwiseProduct(corr_term); - if (std::isnan(xmin_) && !std::isnan(xmax_)) + if (std::isnan(xmin) && !std::isnan(xmax)) f_corr.reverseInPlace(); return f_corr; @@ -965,32 +974,32 @@ Kde1d::boundary_correct(const Eigen::VectorXd& x, const Eigen::VectorXd& fhat) inline Eigen::VectorXd Kde1d::construct_grid_points(const Eigen::VectorXd& x) { + const double xmin = get_effective_xmin(); + const double xmax = get_effective_xmax(); Eigen::VectorXd rng(2); rng << x.minCoeff(), x.maxCoeff(); - if (type_ == VarType::discrete) { + if (type_ == VarType::discrete && std::isnan(xmin) && std::isnan(xmax)) { // Discrete estimates use jittered observations without transformation. - } else if (std::isnan(xmin_) && std::isnan(xmax_)) { + } else if (std::isnan(xmin) && std::isnan(xmax)) { rng(0) -= 4 * bandwidth_; rng(1) += 4 * bandwidth_; - } else if (!std::isnan(xmin_) && !std::isnan(xmax_)) { + } else if (!std::isnan(xmin) && !std::isnan(xmax)) { Eigen::VectorXd boundaries(2); - boundaries << xmin_, xmax_; + boundaries << xmin, xmax; rng = boundary_transform(boundaries); } else { rng(0) = boundary_transform(Eigen::VectorXd::Constant( - 1, std::isnan(xmin_) ? xmax_ : xmin_))(0); + 1, std::isnan(xmin) ? xmax : xmin))(0); rng(1) += 4 * bandwidth_; } auto zgrid = Eigen::VectorXd::LinSpaced(grid_size_ + 1, rng(0), rng(1)); Eigen::VectorXd grid_points = boundary_transform(zgrid, true); // Avoid round-off at finite support boundaries before evaluating the fit. - if (type_ != VarType::discrete) { - if (!std::isnan(xmin_)) - grid_points(0) = xmin_; - if (!std::isnan(xmax_)) - grid_points(std::isnan(xmin_) ? 0 : grid_size_) = xmax_; - } + if (!std::isnan(xmin)) + grid_points(0) = xmin; + if (!std::isnan(xmax)) + grid_points(std::isnan(xmin) ? 0 : grid_size_) = xmax; return grid_points; } @@ -1000,14 +1009,10 @@ Kde1d::construct_grid_points(const Eigen::VectorXd& x) inline Eigen::VectorXd Kde1d::finalize_grid(Eigen::VectorXd& grid_points) { - if (std::isnan(xmin_) && !std::isnan(xmax_)) + const double xmin = get_effective_xmin(); + const double xmax = get_effective_xmax(); + if (std::isnan(xmin) && !std::isnan(xmax)) grid_points.reverseInPlace(); - if (type_ == VarType::discrete) { - if (!std::isnan(xmin_)) - grid_points(0) = xmin_; - if (!std::isnan(xmax_)) - grid_points(grid_points.size() - 1) = xmax_; - } return grid_points; } @@ -1080,10 +1085,12 @@ Kde1d::prepare_endpoint(const Eigen::VectorXd& x, EndpointSide side) const { EndpointData endpoint; + const double xmin = get_effective_xmin(); + const double xmax = get_effective_xmax(); if (side == EndpointSide::lower) - endpoint.dist = x.array() - xmin_; + endpoint.dist = x.array() - xmin; else - endpoint.dist = xmax_ - x.array(); + endpoint.dist = xmax - x.array(); endpoint.weights.resize(x.size()); const Eigen::VectorXi order = tools::get_order(endpoint.dist); const Eigen::VectorXd unsorted_dist = endpoint.dist; @@ -1342,6 +1349,8 @@ Kde1d::repair_boundaries(const Eigen::VectorXd& x, Eigen::VectorXd& influences, const Eigen::VectorXd& weights) { + const double xmin = get_effective_xmin(); + const double xmax = get_effective_xmax(); const Eigen::VectorXd case_weights = weights.size() > 0 ? weights : Eigen::VectorXd::Ones(x.size()); const double n_eff = @@ -1351,9 +1360,9 @@ Kde1d::repair_boundaries(const Eigen::VectorXd& x, EndpointData lower_endpoint; EndpointData upper_endpoint; - if (!std::isnan(xmin_)) + if (!std::isnan(xmin)) lower_endpoint = prepare_endpoint(x, case_weights, EndpointSide::lower); - if (!std::isnan(xmax_)) + if (!std::isnan(xmax)) upper_endpoint = prepare_endpoint(x, case_weights, EndpointSide::upper); if (!lower_endpoint.finite && !upper_endpoint.finite) return; @@ -1380,7 +1389,7 @@ Kde1d::repair_boundaries(const Eigen::VectorXd& x, // Reflection invariance gives both endpoints the same boundary bandwidth. const double bandwidth_fraction = - (!std::isnan(xmin_) && !std::isnan(xmax_)) ? 1.0 : 0.75; + (!std::isnan(xmin) && !std::isnan(xmax)) ? 1.0 : 0.75; const EndpointData& bandwidth_endpoint = lower_endpoint.finite ? lower_endpoint : upper_endpoint; const double h = @@ -1390,13 +1399,13 @@ Kde1d::repair_boundaries(const Eigen::VectorXd& x, BoundaryComponent upper; if (lower_endpoint.finite && n_lower > 0) { lower = fit_boundary_component(lower_endpoint.dist, - grid.head(n_lower).array() - xmin_, + grid.head(n_lower).array() - xmin, h, lower_endpoint.weights); } if (upper_endpoint.finite && n_upper > 0) { Eigen::VectorXd upper_eval_dist = - (xmax_ - grid.tail(n_upper).array()).reverse(); + (xmax - grid.tail(n_upper).array()).reverse(); upper = fit_boundary_component( upper_endpoint.dist, upper_eval_dist, h, upper_endpoint.weights); } @@ -1426,7 +1435,8 @@ Kde1d::select_bandwidth(const Eigen::VectorXd& x, } bandwidth *= multiplier; - if (type_ == VarType::discrete) { + if (type_ == VarType::discrete && std::isnan(xmin_) && + std::isnan(xmax_)) { bandwidth = std::max(bandwidth, 0.5 / 5); } @@ -1436,8 +1446,66 @@ Kde1d::select_bandwidth(const Eigen::VectorXd& x, inline void Kde1d::check_xmin_xmax(const double& xmin, const double& xmax) const { - if (!std::isnan(xmax) && !std::isnan(xmax) && (xmin > xmax)) + if (!std::isnan(xmin) && !std::isnan(xmax) && (xmin > xmax)) throw std::invalid_argument("xmin must be smaller than xmax"); + + if (type_ == VarType::discrete) { + const auto invalid_bound = [](double bound) { + return !std::isnan(bound) && + (!std::isfinite(bound) || bound < 0.0 || + bound != std::floor(bound)); + }; + if (invalid_bound(xmin) || invalid_bound(xmax)) + throw std::invalid_argument( + "discrete bounds must be nonnegative integers"); + } +} + +inline void +Kde1d::check_discrete_data(const Eigen::VectorXd& x) const +{ + if (type_ != VarType::discrete) + return; + + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (!std::isnan(x(i)) && + (!std::isfinite(x(i)) || x(i) < 0.0 || + x(i) != std::floor(x(i)))) { + throw std::invalid_argument( + "discrete data must be nonnegative integers"); + } + } +} + +inline double +Kde1d::get_effective_xmin() const +{ + if (type_ == VarType::discrete && !std::isnan(xmin_)) + return xmin_ - 0.5; + return xmin_; +} + +inline double +Kde1d::get_effective_xmax() const +{ + if (type_ == VarType::discrete && !std::isnan(xmax_)) + return xmax_ + 0.5; + return xmax_; +} + +inline double +Kde1d::get_discrete_support_min() const +{ + return std::isnan(xmin_) ? 0.0 : xmin_; +} + +inline double +Kde1d::get_discrete_support_max() const +{ + const double lower = get_discrete_support_min(); + return std::isnan(xmax_) + ? std::max(lower, std::ceil(grid_.get_grid_max())) + : xmax_; } inline void diff --git a/test/test.cpp b/test/test.cpp index f1a3ab6..a095585 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -434,6 +434,100 @@ TEST_CASE("discrete CDF stays within the unit interval", "[discrete]") CHECK(discrete.cdf(Eigen::VectorXd::Constant(1, 8.0))(0) == 1.0); } +TEST_CASE("discrete bounds describe integer support", "[discrete][boundary]") +{ + Eigen::VectorXd observations(200); + for (Eigen::Index i = 0; i < observations.size(); ++i) + observations(i) = static_cast(i % 4); + + Kde1d bounded(0.0, 3.0, "discrete", 1.0, NAN, 2, 100, true); + bounded.fit(observations); + const Eigen::VectorXd grid = bounded.get_grid_points(); + CHECK(grid(0) == Approx(-0.5)); + CHECK(grid(grid.size() - 1) == Approx(3.5)); + CHECK((grid.tail(grid.size() - 1) - grid.head(grid.size() - 1)) + .minCoeff() > 0.0); + + const Eigen::VectorXd levels = Eigen::VectorXd::LinSpaced(4, 0.0, 3.0); + CHECK(bounded.pdf(levels).sum() == Approx(1.0)); + CHECK(bounded.cdf(Eigen::VectorXd::Constant(1, -0.1))(0) == 0.0); + CHECK(bounded.cdf(Eigen::VectorXd::Constant(1, 3.0))(0) == 1.0); + Eigen::VectorXd endpoints(2); + endpoints << 0.0, 1.0; + CHECK(bounded.quantile(endpoints)(0) == 0.0); + CHECK(bounded.quantile(endpoints)(1) == 3.0); + + Kde1d left_bounded(0.0, NAN, "discrete", 1.0, NAN, 2, 100); + left_bounded.fit(observations); + const Eigen::VectorXd left_grid = left_bounded.get_grid_points(); + CHECK(left_grid(0) == Approx(-0.5)); + CHECK((left_grid.tail(left_grid.size() - 1) - + left_grid.head(left_grid.size() - 1)) + .minCoeff() > 0.0); + const double left_upper = std::ceil(left_grid(left_grid.size() - 1)); + const Eigen::VectorXd left_levels = Eigen::VectorXd::LinSpaced( + static_cast(left_upper + 1.0), 0.0, left_upper); + CHECK(left_bounded.pdf(left_levels).sum() == Approx(1.0)); + + Kde1d right_bounded(NAN, 3.0, "discrete", 1.0, NAN, 2, 100); + right_bounded.fit(observations); + const Eigen::VectorXd right_grid = right_bounded.get_grid_points(); + CHECK(right_grid(right_grid.size() - 1) == Approx(3.5)); + CHECK((right_grid.tail(right_grid.size() - 1) - + right_grid.head(right_grid.size() - 1)) + .minCoeff() > 0.0); + CHECK(right_bounded.pdf(levels).sum() == Approx(1.0)); + CHECK(right_bounded.pdf(Eigen::VectorXd::Constant(1, -1.0))(0) == 0.0); + + Kde1d unbounded(NAN, NAN, "discrete"); + unbounded.fit(observations); + CHECK(unbounded.pdf(Eigen::VectorXd::Constant(1, -1.0))(0) == 0.0); + CHECK(unbounded.cdf(Eigen::VectorXd::Constant(1, -0.1))(0) == 0.0); + CHECK(unbounded.quantile(Eigen::VectorXd::Constant(1, 0.0))(0) == 0.0); + + Kde1d singleton(0.0, 0.0, "discrete"); + singleton.fit(Eigen::VectorXd::Zero(40)); + CHECK(singleton.pdf(Eigen::VectorXd::Constant(1, 0.0))(0) == 1.0); + CHECK(singleton.get_grid_points()(0) == Approx(-0.5)); + CHECK(singleton.get_grid_points()(singleton.get_actual_grid_size() - 1) == + Approx(0.5)); +} + +TEST_CASE("bounded discrete fits use adaptive boundary experts", + "[discrete][boundary-expert]") +{ + Eigen::VectorXd observations(200); + for (Eigen::Index i = 0; i < observations.size(); ++i) + observations(i) = static_cast(i % 4); + + Kde1d repaired(0.0, 3.0, "discrete", 1.0, NAN, 2, 100, true); + repaired.fit(observations); + Kde1d bulk(0.0, 3.0, "discrete", 1.0, NAN, 2, 100, false); + bulk.fit(observations); + + CHECK_FALSE(repaired.get_values().isApprox(bulk.get_values(), 1e-8)); + CHECK(repaired.get_values().array().isFinite().all()); + CHECK(repaired.get_values().minCoeff() >= 0.0); +} + +TEST_CASE("discrete inputs are nonnegative integers", "[discrete][input-checks]") +{ + CHECK_THROWS(Kde1d(-1.0, NAN, "discrete")); + CHECK_THROWS(Kde1d(0.5, NAN, "discrete")); + CHECK_THROWS(Kde1d(NAN, 3.5, "discrete")); + + Kde1d fit(NAN, NAN, "discrete"); + Eigen::VectorXd negative(2); + negative << 0.0, -1.0; + CHECK_THROWS(fit.fit(negative)); + Eigen::VectorXd fractional(2); + fractional << 0.0, 1.5; + CHECK_THROWS(fit.fit(fractional)); + Eigen::VectorXd infinite(2); + infinite << 0.0, std::numeric_limits::infinity(); + CHECK_THROWS(fit.fit(infinite)); +} + TEST_CASE("likelihood summaries match fitted densities", "[likelihood]") { SECTION("weighted likelihood") From 57fb8ab78ad772651be01966ccefbfdab94efed8 Mon Sep 17 00:00:00 2001 From: tnagler Date: Tue, 25 Aug 2026 13:59:44 +0200 Subject: [PATCH 2/3] Allow negative discrete support --- docs/overview-discrete.dox | 16 +++++++-------- include/kde1d/kde1d.hpp | 14 +++++-------- test/test.cpp | 40 ++++++++++++++++++++++++++++---------- 3 files changed, 43 insertions(+), 27 deletions(-) diff --git a/docs/overview-discrete.dox b/docs/overview-discrete.dox index 19211c3..dfe13a7 100644 --- a/docs/overview-discrete.dox +++ b/docs/overview-discrete.dox @@ -1,9 +1,9 @@ /** @page overview-discrete Discrete Distributions -Discrete fits model nonnegative integer-valued observations. Both the data and -any finite `xmin` or `xmax` must therefore be nonnegative integers. A finite -bound describes the support of the discrete variable itself; probabilities -outside those integer bounds are zero. +Discrete fits model integer-valued observations, including negative integers. +Any finite `xmin` or `xmax` must also be an integer. A finite bound describes +the support of the discrete variable itself; probabilities outside those +integer bounds are zero. @section discrete-jitter Jittering and support transforms @@ -37,10 +37,10 @@ levels: \qquad k\in\mathcal S. @f] -Here @f$\mathcal S@f$ is the nonnegative integer support, restricted by any -finite public bounds. Noninteger arguments and integer values outside this -support have probability zero. The CDF and quantile function use the same -normalized masses. +Here @f$\mathcal S@f$ is the integer support inferred from the interpolation +grid and restricted by any finite public bounds. Noninteger arguments and +integer values outside this support have probability zero. The CDF and +quantile function use the same normalized masses. This distinction also explains why an unconstrained jitter density can give a good discrete estimate despite having a nonzero continuous integral outside diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index 9f51c75..d23295d 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -1452,12 +1452,10 @@ Kde1d::check_xmin_xmax(const double& xmin, const double& xmax) const if (type_ == VarType::discrete) { const auto invalid_bound = [](double bound) { return !std::isnan(bound) && - (!std::isfinite(bound) || bound < 0.0 || - bound != std::floor(bound)); + (!std::isfinite(bound) || bound != std::floor(bound)); }; if (invalid_bound(xmin) || invalid_bound(xmax)) - throw std::invalid_argument( - "discrete bounds must be nonnegative integers"); + throw std::invalid_argument("discrete bounds must be integers"); } } @@ -1469,10 +1467,8 @@ Kde1d::check_discrete_data(const Eigen::VectorXd& x) const for (Eigen::Index i = 0; i < x.size(); ++i) { if (!std::isnan(x(i)) && - (!std::isfinite(x(i)) || x(i) < 0.0 || - x(i) != std::floor(x(i)))) { - throw std::invalid_argument( - "discrete data must be nonnegative integers"); + (!std::isfinite(x(i)) || x(i) != std::floor(x(i)))) { + throw std::invalid_argument("discrete data must be integers"); } } } @@ -1496,7 +1492,7 @@ Kde1d::get_effective_xmax() const inline double Kde1d::get_discrete_support_min() const { - return std::isnan(xmin_) ? 0.0 : xmin_; + return std::isnan(xmin_) ? std::floor(grid_.get_grid_min()) : xmin_; } inline double diff --git a/test/test.cpp b/test/test.cpp index a095585..d7fb5fe 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -476,14 +476,37 @@ TEST_CASE("discrete bounds describe integer support", "[discrete][boundary]") CHECK((right_grid.tail(right_grid.size() - 1) - right_grid.head(right_grid.size() - 1)) .minCoeff() > 0.0); - CHECK(right_bounded.pdf(levels).sum() == Approx(1.0)); - CHECK(right_bounded.pdf(Eigen::VectorXd::Constant(1, -1.0))(0) == 0.0); + const double right_lower = std::floor(right_grid(0)); + const Eigen::VectorXd right_levels = Eigen::VectorXd::LinSpaced( + static_cast(3.0 - right_lower + 1.0), right_lower, 3.0); + CHECK(right_bounded.pdf(right_levels).sum() == Approx(1.0)); + CHECK(right_bounded.pdf(Eigen::VectorXd::Constant(1, 4.0))(0) == 0.0); Kde1d unbounded(NAN, NAN, "discrete"); unbounded.fit(observations); - CHECK(unbounded.pdf(Eigen::VectorXd::Constant(1, -1.0))(0) == 0.0); - CHECK(unbounded.cdf(Eigen::VectorXd::Constant(1, -0.1))(0) == 0.0); - CHECK(unbounded.quantile(Eigen::VectorXd::Constant(1, 0.0))(0) == 0.0); + const Eigen::VectorXd unbounded_grid = unbounded.get_grid_points(); + const double unbounded_lower = std::floor(unbounded_grid(0)); + const double unbounded_upper = + std::ceil(unbounded_grid(unbounded_grid.size() - 1)); + const Eigen::VectorXd unbounded_levels = Eigen::VectorXd::LinSpaced( + static_cast(unbounded_upper - unbounded_lower + 1.0), + unbounded_lower, + unbounded_upper); + CHECK(unbounded.pdf(unbounded_levels).sum() == Approx(1.0)); + CHECK(unbounded.quantile(Eigen::VectorXd::Constant(1, 0.0))(0) == + unbounded_lower); + + const Eigen::VectorXd signed_observations = observations.array() - 2.0; + Kde1d signed_bounded(-2.0, 1.0, "discrete", 1.0, NAN, 2, 100); + signed_bounded.fit(signed_observations); + const Eigen::VectorXd signed_grid = signed_bounded.get_grid_points(); + CHECK(signed_grid(0) == Approx(-2.5)); + CHECK(signed_grid(signed_grid.size() - 1) == Approx(1.5)); + const Eigen::VectorXd signed_levels = + Eigen::VectorXd::LinSpaced(4, -2.0, 1.0); + CHECK(signed_bounded.pdf(signed_levels).sum() == Approx(1.0)); + CHECK(signed_bounded.pdf(Eigen::VectorXd::Constant(1, -3.0))(0) == 0.0); + CHECK(signed_bounded.pdf(Eigen::VectorXd::Constant(1, 2.0))(0) == 0.0); Kde1d singleton(0.0, 0.0, "discrete"); singleton.fit(Eigen::VectorXd::Zero(40)); @@ -510,16 +533,13 @@ TEST_CASE("bounded discrete fits use adaptive boundary experts", CHECK(repaired.get_values().minCoeff() >= 0.0); } -TEST_CASE("discrete inputs are nonnegative integers", "[discrete][input-checks]") +TEST_CASE("discrete inputs are integers", "[discrete][input-checks]") { - CHECK_THROWS(Kde1d(-1.0, NAN, "discrete")); + CHECK_NOTHROW(Kde1d(-3.0, -1.0, "discrete")); CHECK_THROWS(Kde1d(0.5, NAN, "discrete")); CHECK_THROWS(Kde1d(NAN, 3.5, "discrete")); Kde1d fit(NAN, NAN, "discrete"); - Eigen::VectorXd negative(2); - negative << 0.0, -1.0; - CHECK_THROWS(fit.fit(negative)); Eigen::VectorXd fractional(2); fractional << 0.0, 1.5; CHECK_THROWS(fit.fit(fractional)); From d07006c30fe0b5dae52ba19d60b532c2b9bafe4b Mon Sep 17 00:00:00 2001 From: tnagler Date: Tue, 25 Aug 2026 14:08:48 +0200 Subject: [PATCH 3/3] Apply bounds to zero-inflated continuous part --- docs/overview-continuous.dox | 17 ++++++++++++-- include/kde1d/kde1d.hpp | 27 +++++++++++++--------- test/test.cpp | 43 ++++++++++++++++++++++++++++++++++++ 3 files changed, 75 insertions(+), 12 deletions(-) diff --git a/docs/overview-continuous.dox b/docs/overview-continuous.dox index fa181c0..b5fdbdb 100644 --- a/docs/overview-continuous.dox +++ b/docs/overview-continuous.dox @@ -37,12 +37,25 @@ An explicitly supplied `bandwidth` is interpreted on the transformed scale. `get_bandwidth()` reports that bulk bandwidth before `multiplier` is applied; it does not report an endpoint-expert bandwidth. +@section continuous-zero-inflated Zero-inflated fits + +For a zero-inflated fit, exact zeros form a separate point mass and the +remaining observations form a continuous hurdle component. Finite `xmin` and +`xmax` constrain only this continuous component: zeros need not lie in the +interval, while every nonzero observation must. After zeros are removed, the +continuous observations use the same support transform, normalization, and +adaptive boundary repair as an ordinary continuous fit. The bounds are used +as supplied, without the half-unit adjustment used for discrete jitter cells. +If all observations are zero, the point mass is one and no continuous +component is fitted. + @section continuous-repair Boundary repair `boundary_repair=true` makes each known finite endpoint eligible for repair; it does not force an expert to be used. Repair is considered only for -continuous or discrete samples with at least 16 cleaned observations and -effective sample size +continuous samples, bounded discrete jitter samples, or nonzero components of +zero-inflated samples with at least 16 cleaned observations and effective +sample size @f[ n_e=\frac{(\sum_i w_i)^2}{\sum_i w_i^2} diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index d23295d..5ab4d85 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -18,11 +18,11 @@ enum class VarType //! Local-polynomial density estimation in one dimension. //! -//! Continuous and bounded discrete fits use a transformed local-likelihood -//! estimate as their bulk component. When `boundary_repair` is enabled, -//! finite support endpoints may additionally use a nonnegative local-linear -//! boundary kernel. See -//! @ref overview-continuous and @ref overview-discrete for the transforms, +//! Continuous fits, bounded discrete fits, and the continuous component of +//! zero-inflated fits use a transformed local-likelihood bulk estimator. When +//! `boundary_repair` is enabled, finite support endpoints may additionally use +//! a nonnegative local-linear boundary kernel. +//! See @ref overview-continuous and @ref overview-discrete for the transforms, //! endpoint classifier, bandwidths, fusion weights, and EDF approximation. class Kde1d { @@ -448,9 +448,7 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) const double effective_xmin = get_effective_xmin(); const double effective_xmax = get_effective_xmax(); const bool use_boundary_repair = - boundary_repair_ && - (type_ == VarType::continuous || type_ == VarType::discrete) && - xx.size() >= 16 && + boundary_repair_ && xx.size() >= 16 && (!std::isnan(effective_xmin) || !std::isnan(effective_xmax)); const Eigen::VectorXd boundary_observations = use_boundary_repair ? xx : Eigen::VectorXd(); @@ -1535,8 +1533,17 @@ Kde1d::check_inputs(const Eigen::VectorXd& x, inline void Kde1d::check_boundaries(const Eigen::VectorXd& x) const { - if ((x.array() < xmin_).any() || (x.array() > xmax_).any()) { - throw std::invalid_argument("x must be contained in [xmin, xmax]."); + for (Eigen::Index i = 0; i < x.size(); ++i) { + if (std::isnan(x(i)) || + (type_ == VarType::zero_inflated && x(i) == 0.0)) + continue; + if ((!std::isnan(xmin_) && x(i) < xmin_) || + (!std::isnan(xmax_) && x(i) > xmax_)) { + const std::string subject = + type_ == VarType::zero_inflated ? "nonzero x" : "x"; + throw std::invalid_argument(subject + + " must be contained in [xmin, xmax]."); + } } } diff --git a/test/test.cpp b/test/test.cpp index d7fb5fe..a04c819 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -1304,3 +1304,46 @@ TEST_CASE("zero-inflated data", "[zero-inflated]") fit.quantile(stats::simulate_uniform(100, { 5 })).cwiseEqual(0).all()); } } + +TEST_CASE("zero-inflated bounds apply to the continuous component", + "[zero-inflated][boundary]") +{ + Eigen::VectorXd positive = Eigen::VectorXd::Zero(200); + positive.tail(160) = Eigen::VectorXd::LinSpaced(160, 1.01, 1.99); + + Kde1d repaired(1.0, 2.0, "zero-inflated", 1.0, NAN, 2, 100, true); + CHECK_NOTHROW(repaired.fit(positive)); + CHECK(repaired.get_prob0() == Approx(0.2)); + CHECK(repaired.get_grid_points()(0) == Approx(1.0)); + CHECK(repaired.get_grid_points()(repaired.get_actual_grid_size() - 1) == + Approx(2.0)); + CHECK(repaired.pdf(Eigen::VectorXd::Constant(1, 0.0))(0) == Approx(0.2)); + CHECK(repaired.pdf(Eigen::VectorXd::Constant(1, 0.5))(0) == 0.0); + CHECK(repaired.cdf(Eigen::VectorXd::Constant(1, -0.1))(0) == 0.0); + CHECK(repaired.cdf(Eigen::VectorXd::Constant(1, 0.0))(0) == Approx(0.2)); + CHECK(repaired.cdf(Eigen::VectorXd::Constant(1, 0.5))(0) == Approx(0.2)); + + Kde1d bulk(1.0, 2.0, "zero-inflated", 1.0, NAN, 2, 100, false); + bulk.fit(positive); + CHECK_FALSE(repaired.get_values().isApprox(bulk.get_values(), 1e-8)); + + Eigen::VectorXd negative = Eigen::VectorXd::Zero(200); + negative.tail(160) = Eigen::VectorXd::LinSpaced(160, -2.99, -2.01); + Kde1d negative_fit(-3.0, -2.0, "zero-inflated", 1.0, NAN, 2, 100); + CHECK_NOTHROW(negative_fit.fit(negative)); + CHECK(negative_fit.get_grid_points()(0) == Approx(-3.0)); + CHECK(negative_fit.get_grid_points()(negative_fit.get_actual_grid_size() - + 1) == Approx(-2.0)); + CHECK(negative_fit.pdf(Eigen::VectorXd::Constant(1, 0.0))(0) == Approx(0.2)); + CHECK(negative_fit.pdf(Eigen::VectorXd::Constant(1, -1.5))(0) == 0.0); + + Eigen::VectorXd invalid(3); + invalid << 0.0, 0.5, 1.5; + Kde1d invalid_fit(1.0, 2.0, "zero-inflated"); + CHECK_THROWS(invalid_fit.fit(invalid)); + + Kde1d point_mass(1.0, 2.0, "zero-inflated"); + CHECK_NOTHROW(point_mass.fit(Eigen::VectorXd::Zero(40))); + CHECK(point_mass.get_prob0() == 1.0); + CHECK(point_mass.pdf(Eigen::VectorXd::Constant(1, 0.0))(0) == 1.0); +}