calculate inverse cdf of the Gamma distribution by Halley's method wi…

…th improved initial values New method is significantly faster due to the better initial values and more economic calculation.
rabauke · Feb 22, 2024 · d858f08 · d858f08
1 parent 1af9373
commit d858f08
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 23 deletions.
diff --git a/trng/gamma_dist.hpp b/trng/gamma_dist.hpp
@@ -118,23 +118,7 @@ namespace trng {
     // inverse cumulative density function
     TRNG_CUDA_ENABLE
     result_type icdf_(result_type x) const {
-      if (x <= math::numeric_limits<result_type>::epsilon())
-        return 0;
-      if (P.kappa() == 1)  // special case of exponential distribution
-        return -math::ln(1 - x) * P.theta();
-      const result_type ln_Gamma_kappa{math::ln_Gamma(P.kappa())};
-      result_type y{P.kappa()}, y_old;
-      int num_iterations{0};
-      do {
-        ++num_iterations;
-        y_old = y;
-        const result_type f0{math::GammaP(P.kappa(), y) - x};
-        const result_type f1{math::exp((P.kappa() - 1) * math::ln(y) - y - ln_Gamma_kappa)};
-        const result_type f2{f1 * (P.kappa() - 1 - y) / y};
-        y -= f0 / f1 * (1 + f0 * f2 / (2 * f1 * f1));
-      } while (num_iterations < 16 &&
-               math::abs((y - y_old) / y) > 16 * math::numeric_limits<result_type>::epsilon());
-      return y * P.theta();
+      return math::inv_GammaP(P.kappa(), x) * P.theta();
     }
 
   public:
@@ -202,7 +186,7 @@ namespace trng {
         return 0;
       if (x == 1)
         return math::numeric_limits<result_type>::infinity();
-      return icdf_(x);
+      return math::inv_GammaP(P.kappa(), x) * P.theta();
     }
   };
 

diff --git a/trng/special_functions.hpp b/trng/special_functions.hpp
@@ -586,10 +586,10 @@ namespace trng {
       template<typename T>
       TRNG_CUDA_ENABLE T inv_GammaP(T a, T p) {
         const T eps{sqrt(numeric_limits<T>::epsilon())};
-        T a1{a - 1};
-        T glna{ln_Gamma(a)};
-        T lna1{ln(a1)};
-        T afac{exp(a1 * (lna1 - 1) - glna)};
+        const T a1{a - 1};
+        const T glna{ln_Gamma(a)};
+        const T lna1{ln(a1)};
+        const T afac{exp(a1 * (lna1 - 1) - glna)};
         T x;
         // initial guess
         if (a > T{1}) {
@@ -603,7 +603,7 @@ namespace trng {
           x = p < t ? pow(p / t, 1 / a) : 1 - ln1p(-(p - t) / (1 - t));
         }
         // refinement by Halley's method
-        for (int i{0}; i < 32; ++i) {
+        for (int i{0}; i < numeric_limits<T>::digits; ++i) {
           if (x <= 0) {
             x = 0;
             break;