From ef80ad2abbe47c5d73b22669a0d5f994fc3f362b Mon Sep 17 00:00:00 2001 From: tnagler Date: Fri, 11 Jul 2025 15:41:54 +0200 Subject: [PATCH 1/5] prevent /0 nans in edf computation --- include/kde1d/kde1d.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index e46d970..82a1c0d 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -583,7 +583,10 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, grid_points(m - 1), m - 1, Eigen::VectorXd::Ones(x.size())); - wbin = wcount.cwiseQuotient(count); + wbin = (count.array() == 0).select( + Eigen::VectorXd::Zero(m), + wcount.cwiseQuotient(count) + ); } Eigen::MatrixXd res(f0.size(), 2); From a8d0c9ed3393a9218995d15bc195eb89752ea3e9 Mon Sep 17 00:00:00 2001 From: tnagler Date: Fri, 11 Jul 2025 15:43:52 +0200 Subject: [PATCH 2/5] fix loglik for zero-inflated data (obs=0 was omitted) --- include/kde1d/kde1d.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index 82a1c0d..7bbb8ca 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -319,13 +319,16 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) if (type_ == VarType::discrete) { xx = xx.array().round(); } - loglik_ = (this->pdf(xx, false).array().log()).sum(); + loglik_ = (this->pdf(xx, false).array().log().array() * w.array()).sum(); + if (prob0_ > 0) { + loglik_ += x.size() * prob0_ * std::log(prob0_); + } // calculate effective degrees of freedom interp::InterpolationGrid infl_grid( - grid_points, fitted.col(1).cwiseMin(3.0).cwiseMax(0), 0); - Eigen::VectorXd influences = infl_grid.interpolate(xx).array() * (1 - prob0_); - edf_ = influences.sum() + (prob0_ > 0); + grid_points, fitted.col(1).cwiseMin(3.0).cwiseMax(0), 0); + Eigen::VectorXd influences = infl_grid.interpolate(xx).array(); + edf_ = influences.sum() + static_cast(prob0_ > 0); // store bandwidth in standardized format bandwidth_ = bandwidth_ / multiplier_; From 0eda0f0cb83741513dbd340b3c33b4de79845715 Mon Sep 17 00:00:00 2001 From: tnagler Date: Fri, 11 Jul 2025 15:49:26 +0200 Subject: [PATCH 3/5] fix implicit cast --- include/kde1d/kde1d.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index 7bbb8ca..bef06d5 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -321,7 +321,7 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) } loglik_ = (this->pdf(xx, false).array().log().array() * w.array()).sum(); if (prob0_ > 0) { - loglik_ += x.size() * prob0_ * std::log(prob0_); + loglik_ += static_cast(x.size()) * prob0_ * std::log(prob0_); } // calculate effective degrees of freedom From 81860f244452fce4a84a3ae06b8c8fc11c9a54d4 Mon Sep 17 00:00:00 2001 From: tnagler Date: Fri, 11 Jul 2025 16:12:11 +0200 Subject: [PATCH 4/5] safeguard weighted loglik computation --- include/kde1d/kde1d.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index bef06d5..adae0a8 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -270,8 +270,9 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) Eigen::VectorXd w = weights; tools::remove_nans(xx, w); - if (w.size() > 0) + if (w.size() > 0) { w /= w.mean(); + } if (type_ == VarType::zero_inflated) { if (w.size() == 0) @@ -319,6 +320,10 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) if (type_ == VarType::discrete) { xx = xx.array().round(); } + + if (w.size() == 0) { + w = Eigen::VectorXd::Ones(xx.size()); + } loglik_ = (this->pdf(xx, false).array().log().array() * w.array()).sum(); if (prob0_ > 0) { loglik_ += static_cast(x.size()) * prob0_ * std::log(prob0_); @@ -587,7 +592,7 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, m - 1, Eigen::VectorXd::Ones(x.size())); wbin = (count.array() == 0).select( - Eigen::VectorXd::Zero(m), + Eigen::VectorXd::Zero(count.size()), wcount.cwiseQuotient(count) ); } From 2dc6f31ba38c0aad37287e4920a51792344349be Mon Sep 17 00:00:00 2001 From: tnagler Date: Mon, 21 Jul 2025 16:08:34 +0200 Subject: [PATCH 5/5] add comment for loglik --- include/kde1d/kde1d.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index adae0a8..663f792 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -324,8 +324,12 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) if (w.size() == 0) { w = Eigen::VectorXd::Ones(xx.size()); } + loglik_ = (this->pdf(xx, false).array().log().array() * w.array()).sum(); if (prob0_ > 0) { + // For zero inflated data, all observations with value 0 have been removed, + // so their likelihood contribution is missing. There were n * prob0_ such + // observations, each with log-likelihood contribution log(prob0_). loglik_ += static_cast(x.size()) * prob0_ * std::log(prob0_); }