From 0bae39ef9b399ad602760198356bf86d708664b3 Mon Sep 17 00:00:00 2001
From: Vincent Gao <gaobing1230@gmail.com>
Date: Fri, 31 Jul 2026 21:39:02 +0200
Subject: [PATCH] Fix gamma CDF normalization and incomplete gamma fraction

GammaDistribution#cdf was scaled by scale/Gamma(shape), so it returned values outside [0,1] (e.g. 3.0 for shape=2, scale=3) and saturated below 1.0 whenever the two did not cancel. incomplete_gamma is already regularized.

gamma_fraction returned after its first loop iteration and never applied the Lentz update, leaving a 0.5-1.5% error on the x >= a+1 branch. The update now runs every iteration with the XMININ guards, mirroring beta_fraction, and the result is returned after the loop.

Tests cover the [0,1] range for scale != Gamma(shape) and mpmath reference values on both incomplete-gamma branches.
---
 lib/rubystats/gamma_distribution.rb |  2 +-
 lib/rubystats/modules.rb            | 14 +++++-----
 test/tc_gamma.rb                    | 40 ++++++++++++++++++++++++++---
 3 files changed, 45 insertions(+), 11 deletions(-)

diff --git a/lib/rubystats/gamma_distribution.rb b/lib/rubystats/gamma_distribution.rb
index 22e6d37..faff930 100644
--- a/lib/rubystats/gamma_distribution.rb
+++ b/lib/rubystats/gamma_distribution.rb
@@ -36,7 +36,7 @@ def get_pdf(x)
     # return the probability that a stochastic variable x is less then X, i.e. P(x<X).
     def get_cdf(x)
       check_range(x,0.0,MAX_VALUE)
-      @scale * incomplete_gamma(@shape, x/@scale) / Math.gamma(@shape)
+      incomplete_gamma(@shape, x/@scale)
     end
 
     # Private method to obtain single inverse CDF value.
diff --git a/lib/rubystats/modules.rb b/lib/rubystats/modules.rb
index d0ee34c..956bfde 100644
--- a/lib/rubystats/modules.rb
+++ b/lib/rubystats/modules.rb
@@ -414,19 +414,19 @@ def gamma_fraction(a, x)
           an = -i * (i - a)
           b += 2.0
           d  = an * d + b
+          if d.abs < XMININ
+            d = XMININ
+          end
           c  = b + an / c
           if c.abs < XMININ
             c = XMININ
-            if d.abs < XMININ
-              c = XMININ
-              d   = 1.0 / d
-              del = d * c
-              h  *= del
-            end
           end
+          d   = 1.0 / d
+          del = d * c
+          h  *= del
         end
-        return Math.exp(-x + a * Math.log(x) - log_gamma(a)) * h
       end
+      Math.exp(-x + a * Math.log(x) - log_gamma(a)) * h
     end
 
     # Beta function.
diff --git a/test/tc_gamma.rb b/test/tc_gamma.rb
index 9140964..6b48abd 100644
--- a/test/tc_gamma.rb
+++ b/test/tc_gamma.rb
@@ -16,8 +16,8 @@ def test_simple
     assert_in_delta(0.0, gamma.pdf(0.0), 0.000001)
     assert_in_delta(0.1353353, gamma.pdf(4.0), 0.000001)
     assert_in_delta(0.001134998, gamma.pdf(20.0), 0.000001)
-    assert_in_delta(0.5768099, gamma.cdf(6.0), 0.01)
-    assert_in_delta(0.986246, gamma.cdf(16.0), 0.01)
+    assert_in_delta(0.5768099, gamma.cdf(6.0), 0.0000001)
+    assert_in_delta(0.986246, gamma.cdf(16.0), 0.0000001)
     
    # make sure the mean of RNG values is close to the expected mean
     rngs = []
@@ -36,4 +36,38 @@ def test_integer_input
     pdff = Rubystats::GammaDistribution.new(shapef,scalef).pdf(xf)
     assert_in_delta pdfi, pdff, 0.00000001   
   end
-end
\ No newline at end of file
+
+  def test_cdf_stays_in_unit_interval
+    # cdf is P(shape, x/scale); for scale != Gamma(shape) the old
+    # normalization returned values up to 5.0 and saturated below 1.0
+    [[2.0,3.0],[0.5,4.0],[5.0,10.0]].each do |shape, scale|
+      gamma = Rubystats::GammaDistribution.new(shape, scale)
+      prev = 0.0
+      [0.0,0.5,1.0,2.0,5.0,10.0,50.0,500.0,1.0e6].each do |x|
+        cdf = gamma.cdf(x)
+        assert(cdf >= 0.0 && cdf <= 1.0, "cdf(#{x}) = #{cdf} not in [0,1] for shape=#{shape} scale=#{scale}")
+        assert(cdf >= prev, "cdf(#{x}) = #{cdf} < cdf(previous) = #{prev} for shape=#{shape} scale=#{scale}")
+        prev = cdf
+      end
+      assert_in_delta(1.0, gamma.cdf(1.0e6), 0.0000000001)
+    end
+  end
+
+  def test_cdf_matches_reference_values
+    # P(shape, x/scale) from mpmath at 50 dps; covers both the series and
+    # continued-fraction branches of the incomplete gamma
+    [[2.0,3.0,1.0,0.0446249192349477],
+     [2.0,3.0,5.0,0.496331725766502],
+     [2.0,3.0,10.0,0.84541269549524],
+     [3.0,2.0,6.0,0.576809918873156],
+     [3.0,2.0,16.0,0.986246032255997],
+     [5.0,1.0,8.0,0.900367599512954],
+     [10.0,1.0,15.0,0.93014633930059],
+     [20.0,1.0,30.0,0.978126531558609],
+     [0.5,1.0,1.5,0.91673548333645],
+     [1.5,1.0,2.6,0.842275549603337],
+     [7.5,2.0,11.0,0.247405629288778]].each do |shape, scale, x, expected|
+      assert_in_delta(expected, Rubystats::GammaDistribution.new(shape, scale).cdf(x), 0.000000001)
+    end
+  end
+end