From e9b0fa2a8f3e7083eec217afa7cadb0e3e4426ca Mon Sep 17 00:00:00 2001 From: tnagler Date: Tue, 11 Mar 2025 16:35:16 +0100 Subject: [PATCH 1/2] fix EDF computations (and make them easier to trace) --- include/kde1d/kde1d.hpp | 38 ++++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/include/kde1d/kde1d.hpp b/include/kde1d/kde1d.hpp index 47bd962..e46d970 100644 --- a/include/kde1d/kde1d.hpp +++ b/include/kde1d/kde1d.hpp @@ -127,7 +127,8 @@ class Kde1d const Eigen::VectorXd& weights); double calculate_infl(const size_t& n, const double& f0, - const double& b, + const double& f1, + const double& f2, const double& bandwidth, const double& s, const double& weight); @@ -322,7 +323,7 @@ Kde1d::fit(const Eigen::VectorXd& x, const Eigen::VectorXd& weights) // calculate effective degrees of freedom interp::InterpolationGrid infl_grid( - grid_points, fitted.col(1).cwiseMin(2.0).cwiseMax(0), 0); + 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); @@ -571,6 +572,7 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, fft::KdeFFT kde_fft( x, bandwidth_, grid_points(0), grid_points(m - 1), weights); Eigen::VectorXd f0 = kde_fft.kde_drv(0); + Eigen::VectorXd f1(f0.size()), f2(f0.size()); Eigen::VectorXd wbin = Eigen::VectorXd::Ones(m); if (weights.size()) { @@ -592,11 +594,11 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, return res; // degree > 0 - Eigen::VectorXd f1 = kde_fft.kde_drv(1); + f1 = kde_fft.kde_drv(1); Eigen::VectorXd S = Eigen::VectorXd::Constant(f0.size(), bandwidth_); Eigen::VectorXd b = f1.cwiseQuotient(f0); if (degree_ == 2) { - Eigen::VectorXd f2 = kde_fft.kde_drv(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(); @@ -608,9 +610,8 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, res.col(0) = res.col(0).array() * (-0.5 * b.array().pow(2) * S.array()).exp(); for (size_t k = 0; k < m; k++) { - // TODO: weights res(k, 1) = - calculate_infl(x.size(), f0(k), b(k), bandwidth_, S(k), wbin(k)); + calculate_infl(x.size(), f0(k), f1(k), f2(k), bandwidth_, S(k), wbin(k)); if (std::isnan(res(k, 0))) res.row(k).setZero(); } @@ -623,35 +624,36 @@ Kde1d::fit_lp(const Eigen::VectorXd& x, inline double Kde1d::calculate_infl(const size_t& n, const double& f0, - const double& b, + const double& f1, + const double& f2, const double& bandwidth, const double& s, const double& weight) { double M_inverse00; - double bandwidth2 = std::pow(bandwidth, 2); - double b2 = std::pow(b, 2); + double B = bandwidth * bandwidth; if (degree_ == 0) { M_inverse00 = 1 / f0; } else if (degree_ == 1) { Eigen::Matrix2d M; M(0, 0) = f0; - M(0, 1) = bandwidth2 * b * f0; + M(0, 1) = B * f1; M(1, 0) = M(0, 1); - M(1, 1) = f0 * bandwidth2 + f0 * bandwidth2 * bandwidth2 * b2; + M(1, 1) = f0 * B + B * f1 * f1 * B / f0; M_inverse00 = M.inverse()(0, 0); } else { Eigen::Matrix3d M; M(0, 0) = f0; - M(0, 1) = f0 * b; + M(0, 1) = B * f1; M(1, 0) = M(0, 1); - M(1, 1) = f0 * bandwidth2 + f0 * b2; - M(1, 2) = 0.5 * f0 * (3.0 / s * b + b * b2); + M(1, 1) = B * f2 * B + B * f0; + M(2, 0) = M(1, 1) / 2; + M(0, 2) = M(1, 1) / 2; + double s2 = B * f1 / f0; + M(1, 2) = f0 / 2 * (3 / s * s2 + std::pow(s2, 3)); M(2, 1) = M(1, 2); - M(2, 2) = 0.25 * f0; - M(2, 2) *= 3.0 / std::pow(s, 2) + 6.0 / s * b2 + b2 * b2; - M(0, 2) = M(2, 2); - M(2, 0) = M(2, 2); + M(2, 2) = + f0 / 4 * (3 / (s * s) + 6 / s * std::pow(s2, 2) + std::pow(s2, 4)); M_inverse00 = M.inverse()(0, 0); } From be131251ba8a4614941e6230300d193391066257 Mon Sep 17 00:00:00 2001 From: tnagler Date: Wed, 12 Mar 2025 10:32:24 +0100 Subject: [PATCH 2/2] fix actions script --- .github/workflows/main.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 2a03310..93c2e3c 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -28,7 +28,7 @@ jobs: submodules: recursive - name: Install dependencies id: install-dependencies - uses: vinecopulib/vinecopulib/.github/actions/install-dependencies@better-actions + uses: vinecopulib/vinecopulib/.github/actions/install-dependencies@main with: os: ${{ matrix.cfg.os }} platform: ${{ matrix.cfg.platform }}