ENH: add geninvgauss and genhyperbolic pdf kernels

Author: fbourgey

include/xsf/cpu/stats.h

@@ -159,5 +159,81 @@ namespace cpu {
         return std::min(std::max(cdf, 0.0), 1.0);
     }
 
+    // Logarithm of the generalized inverse Gaussian probability density function.
+    inline double geninvgauss_logpdf(double x, double p, double b) {
+        if (x <= 0) {
+            return -std::numeric_limits<double>::infinity();
+        }
+
+        double z = cyl_bessel_ke(p, b);
+        if (std::isinf(z)) {
+            return std::numeric_limits<double>::quiet_NaN();
+        }
+
+        return -std::log(2.0) - std::log(z) + b + (p - 1.0) * std::log(x) - b * (x + 1.0 / x) / 2.0;
+    }
+
+    inline float geninvgauss_logpdf(float x, float p, float b) {
+        return static_cast<float>(
+            geninvgauss_logpdf(static_cast<double>(x), static_cast<double>(p), static_cast<double>(b))
+        );
+    }
+
+    // Generalized inverse Gaussian probability density function.
+    inline double geninvgauss_pdf(double x, double p, double b) { return std::exp(geninvgauss_logpdf(x, p, b)); }
+
+    inline float geninvgauss_pdf(float x, float p, float b) {
+        return static_cast<float>(
+            geninvgauss_pdf(static_cast<double>(x), static_cast<double>(p), static_cast<double>(b))
+        );
+    }
+
+    namespace detail {
+
+        inline double genhyperbolic_log_norming_constant(double p, double a, double b) {
+            double t1, t2, t3, t4, t5, t6;
+
+            t1 = (a + b) * (a - b);
+            t2 = p * 0.5 * std::log(t1);
+            t3 = 0.5 * std::log(2 * M_PI);
+            t4 = (p - 0.5) * std::log(a);
+            t5 = std::sqrt(t1);
+            t6 = std::log(cyl_bessel_ke(p, t5)) - t5;
+
+            return t2 - t3 - t4 - t6;
+        }
+
+    } // namespace detail
+
+    // Logarithm of the generalized hyperbolic probability density function.
+    inline double genhyperbolic_logpdf(double x, double p, double a, double b) {
+        double t1, t2, t3, t4, t5;
+
+        t1 = detail::genhyperbolic_log_norming_constant(p, a, b);
+        t2 = std::sqrt(1.0 + x * x);
+        t3 = (p - 0.5) * std::log(t2);
+        t4 = std::log(cyl_bessel_ke(p - 0.5, a * t2)) - a * t2;
+        t5 = b * x;
+
+        return t1 + t3 + t4 + t5;
+    }
+
+    inline float genhyperbolic_logpdf(float x, float p, float a, float b) {
+        return static_cast<float>(genhyperbolic_logpdf(
+            static_cast<double>(x), static_cast<double>(p), static_cast<double>(a), static_cast<double>(b)
+        ));
+    }
+
+    // Generalized hyperbolic probability density function.
+    inline double genhyperbolic_pdf(double x, double p, double a, double b) {
+        return std::exp(genhyperbolic_logpdf(x, p, a, b));
+    }
+
+    inline float genhyperbolic_pdf(float x, float p, float a, float b) {
+        return static_cast<float>(genhyperbolic_pdf(
+            static_cast<double>(x), static_cast<double>(p), static_cast<double>(a), static_cast<double>(b)
+        ));
+    }
+
 } // namespace cpu
 } // namespace xsf

tests/xsf_tests/test_genhyperbolic.cpp

@@ -0,0 +1,87 @@
+#include "../testing_utils.h"
+
+#include <array>
+#include <cmath>
+#include <limits>
+
+#include <xsf/cpu/stats.h>
+
+TEST_CASE("generalized hyperbolic density", "[genhyperbolic][xsf_tests]") {
+    // Mirrors https://github.com/scipy/scipy/blob/v1.18.1/scipy/stats/tests/test_distributions.py#L962-L982
+    constexpr std::array<double, 10> expected_pdf{
+        2.94895678275316e-13, 1.75746848647696e-10, 9.48149804073045e-08, 4.17862521692026e-05, 0.0103947630463822,
+        0.240864958986839,    0.162833527161649,    0.0374609592899472,   0.00634894847327781,  0.000941920705790324,
+    };
+    constexpr double p = 2.0;
+    constexpr double a = 3.0;
+    constexpr double b = 1.5;
+    constexpr double loc = 0.5;
+    constexpr double scale = 1.5;
+
+    const std::vector<double> x_values = linspace(-10.0, 10.0, expected_pdf.size());
+    for (int i = 0; i < expected_pdf.size(); ++i) {
+        double x = x_values[i];
+        double standardized_x = (x - loc) / scale;
+        double pdf = xsf::cpu::genhyperbolic_pdf(standardized_x, p, a, b) / scale;
+        CAPTURE(i, x, standardized_x, pdf, expected_pdf[i]);
+        REQUIRE(xsf::extended_relative_error(pdf, expected_pdf[i]) <= 1e-13);
+    }
+}
+
+TEST_CASE("generalized hyperbolic Student's t limit", "[genhyperbolic][xsf_tests]") {
+    // Mirrors https://github.com/scipy/scipy/blob/v1.18.1/scipy/stats/tests/test_distributions.py#L1081-L1097
+    constexpr std::array<double, 10> lower{
+        -31.820519750798752, -3.640296435003648, -2.9487520061496664, -2.7322159231076126, -2.6271579957950233,
+        -2.565225109345393,  -2.524412431436311, -2.495498782411984,  -2.4739457843399073, -2.4572615423796655,
+    };
+    constexpr std::array<double, 10> upper{
+        31.820140377530542, 3.640296435001152,  2.9487520061496824, 2.732215923107273, 2.627157995794383,
+        2.5652251093455574, 2.5244124314366156, 2.4954987824114236, 2.473945784340125, 2.4572615423800115,
+    };
+    constexpr double alpha_epsilon = std::numeric_limits<float>::epsilon();
+
+    const std::vector<double> degrees_of_freedom = linspace(1.0, 30.0, lower.size());
+    for (int j = 0; j < degrees_of_freedom.size(); ++j) {
+        double df = degrees_of_freedom[j];
+        double p = -df / 2.0;
+        double a = df * df * alpha_epsilon;
+        double scale = std::sqrt(df);
+
+        for (const double x : linspace(lower[j], upper[j], 50)) {
+            double pdf = xsf::cpu::genhyperbolic_pdf(x / scale, p, a, 0.0) / scale;
+            double expected_pdf = std::tgamma((df + 1.0) / 2.0) / (std::sqrt(df * M_PI) * std::tgamma(df / 2.0)) *
+                                  std::pow(1.0 + x * x / df, -(df + 1.0) / 2.0);
+            CAPTURE(j, df, x, pdf, expected_pdf);
+            REQUIRE(xsf::extended_relative_error(pdf, expected_pdf) <= 1e-6);
+        }
+    }
+}
+
+TEST_CASE("generalized hyperbolic Cauchy limit", "[genhyperbolic][xsf_tests]") {
+    // Mirrors https://github.com/scipy/scipy/blob/v1.18.1/scipy/stats/tests/test_distributions.py#L1099-L1114
+    constexpr double p = -0.5;
+    constexpr double a = std::numeric_limits<float>::epsilon();
+    constexpr double lower = -31.820519750798752;
+    constexpr double upper = 31.820140377530542;
+
+    for (const double x : linspace(lower, upper, 50)) {
+        double pdf = xsf::cpu::genhyperbolic_pdf(x, p, a, 0.0);
+        double expected_pdf = 1.0 / (M_PI * (1.0 + x * x));
+        CAPTURE(x, pdf, expected_pdf);
+        REQUIRE(xsf::extended_relative_error(pdf, expected_pdf) <= 1e-6);
+    }
+}
+
+TEST_CASE("generalized hyperbolic Laplace limit", "[genhyperbolic][xsf_tests]") {
+    // Mirrors https://github.com/scipy/scipy/blob/v1.18.1/scipy/stats/tests/test_distributions.py#L1116-L1135
+    constexpr double scale = std::numeric_limits<float>::epsilon();
+
+    for (const double loc : linspace(-10.0, 10.0, 10)) {
+        for (const double x : linspace(-20.0, 20.0, 50)) {
+            double pdf = xsf::cpu::genhyperbolic_pdf((x - loc) / scale, 1.0, scale, 0.0) / scale;
+            double expected_pdf = 0.5 * std::exp(-std::abs(x - loc));
+            CAPTURE(loc, x, pdf, expected_pdf);
+            REQUIRE(xsf::extended_relative_error(pdf, expected_pdf) <= 1e-11);
+        }
+    }
+}

tests/xsf_tests/test_geninvgauss.cpp

@@ -0,0 +1,26 @@
+#include "../testing_utils.h"
+
+#include <array>
+#include <cmath>
+
+#include <xsf/cpu/stats.h>
+
+TEST_CASE("generalized inverse Gaussian density", "[geninvgauss][xsf_tests]") {
+    // Mirrors https://github.com/scipy/scipy/blob/v1.18.1/scipy/stats/tests/test_distributions.py#L941-L950
+    constexpr std::array<double, 10> expected_pdf{
+        2.081176820e-21, 4.488660034e-01, 3.747774338e-01, 2.693297528e-01, 1.905637275e-01,
+        1.351476913e-01, 9.636538981e-02, 6.909040154e-02, 4.978006801e-02, 3.602084467e-02,
+    };
+    const std::vector<double> x = linspace(0.01, 5.0, 10);
+    for (int i = 0; i < x.size(); ++i) {
+        double pdf = xsf::cpu::geninvgauss_pdf(x[i], 0.5, 1.0);
+        CAPTURE(i, x[i], pdf, expected_pdf[i]);
+        REQUIRE(xsf::extended_relative_error(pdf, expected_pdf[i]) <= 1e-8);
+    }
+}
+
+TEST_CASE("generalized inverse Gaussian PDF support", "[geninvgauss][xsf_tests]") {
+    // Mirrors https://github.com/scipy/scipy/blob/v1.18.1/scipy/stats/tests/test_distributions.py#L952-L957
+    REQUIRE(xsf::cpu::geninvgauss_pdf(0.0, 0.5, 0.5) == 0.0);
+    REQUIRE(xsf::cpu::geninvgauss_pdf(2e6, 50.0, 2.0) == 0.0);
+}