Fix gamma CDF normalization and incomplete gamma fraction

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.

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.

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
+
+  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