Issue-113: Lower incomplete gamma function cannot be calculated properly for certain chi squared values (#159)

* gemspec: Bump rspec requirements.

* feat: Implement a normalised lower incomplete gamma function.

This change implements the normalised lower incomplete gamma function to
be used for the Chi Squared distribution, improving the speed and
accurary of the results for the Chi Squared tests too.

* fix: Use new implementation of the normalised lower incomplete gamma.

This change fixes issue #113 where it couldn't be calculated due to
overflow operations when using the lower incomplete gamma function that
relies on the composite simpson's rule computation.

* spec: Update Chi Squared Test specs to reflect new changes.

* spec: Add specs for `BigDecimal` calculations.

* spec: Update number expectations to use delta matchers.

This change allow us to remove the use of the `round` function when
comparing float/decimal numbers, giving us a more robust comparison
framework.

* feat: add chi square distribution table.

* refactor: Make minor changes to methods to make it more ruby-friendly.

---------

Co-authored-by: shamangeorge <shamangeorge@fruitopology.net>

Author: estebanz01

lib/math.rb

@@ -57,6 +57,53 @@ def self.lower_incomplete_gamma_function(s, x)
     end
   end
 
+  # Algorithm implementation translated from the ASA147 C++ version https://people.sc.fsu.edu/~jburkardt/cpp_src/asa147/asa147.html
+  # translated from FORTRAN by John Burkardt. Original algorithm written by Chi Leung Lau.
+  # It contains a modification on the error and underflow parameters to use maximum available float number
+  # and it performs the series using `Rational` objects to avoid memory exhaustion and reducing precision errors.
+  #
+  # This algorithm is licensed with MIT license.
+  #
+  # Reference:
+  #
+  #    Chi Leung Lau,
+  #    Algorithm AS 147:
+  #    A Simple Series for the Incomplete Gamma Integral,
+  #    Applied Statistics,
+  #    Volume 29, Number 1, 1980, pages 113-114.
+  def self.normalised_lower_incomplete_gamma_function(s, x)
+    return 0.0 if s.negative? || x.zero? || x.negative?
+
+    # e = 1.0e-308
+    # uflo = 1.0e-47
+    e = Float::MIN
+    uflo = Float::MIN
+
+    lgamma, sign = Math.lgamma(s + 1.0)
+    arg = s * Math.log(x) - (sign * lgamma) - x
+
+    return 0.0 if arg < Math.log(uflo)
+
+    f = Math.exp(arg).to_r
+
+    return 0.0 if f.zero?
+
+    c = 1r
+    value = 1r
+    a = s.to_r
+    rational_x = x.to_r
+
+    loop do
+      a += 1r
+      c = c * (rational_x / a)
+      value += c
+
+      break if c <= (e * value)
+    end
+
+    (value * f).to_f
+  end
+
   def self.beta_function(x, y)
     return 1 if x == 1 && y == 1
 

lib/ruby-statistics/distribution/chi_squared.rb

@@ -15,7 +15,8 @@ def cumulative_function(value)
           1.0 - Math.exp((-1.0 * value / 2.0))
         else
           k = degrees_of_freedom/2.0
-          Math.lower_incomplete_gamma_function(k, value/2.0)/Math.gamma(k)
+          # Math.lower_incomplete_gamma_function(k, value/2.0)/Math.gamma(k)
+          Math.normalised_lower_incomplete_gamma_function(k, value / 2.0)
         end
       end
 
@@ -27,7 +28,7 @@ def density_function(value)
         left_down = (2 ** common) * Math.gamma(common)
         right = (value ** (common - 1)) * Math.exp(-(value/2.0))
 
-        (1.0/left_down) * right
+        right / left_down
       end
 
       def mode

lib/ruby-statistics/distribution/tables/chi_squared.rb

@@ -0,0 +1,88 @@
+module RubyStatistics
+  module Distribution
+    module Tables
+      class ChiSquared
+        # Values retrieved from the following table provided by the University of Arizona.
+        # https://math.arizona.edu/~jwatkins/chi-square-table.pdf
+        TABLE = [
+          [0.000, 0.000, 0.001, 0.004, 0.016, 2.706, 3.841, 5.024, 6.635, 7.879],
+          [0.010, 0.020, 0.051, 0.103, 0.211, 4.605, 5.991, 7.378, 9.210, 10.597],
+          [0.072, 0.115, 0.216, 0.352, 0.584, 6.251, 7.815, 9.348, 11.345, 12.838],
+          [0.207, 0.297, 0.484, 0.711, 1.064, 7.779, 9.488, 11.143, 13.277, 14.860],
+          [0.412, 0.554, 0.831, 1.145, 1.610, 9.236, 11.070, 12.833, 15.086, 16.750],
+          [0.676, 0.872, 1.237, 1.635, 2.204, 10.645, 12.592, 14.449, 16.812, 18.548],
+          [0.989, 1.239, 1.690, 2.167, 2.833, 12.017, 14.067, 16.013, 18.475, 20.278],
+          [1.344, 1.646, 2.180, 2.733, 3.490, 13.362, 15.507, 17.535, 20.090, 21.955],
+          [1.735, 2.088, 2.700, 3.325, 4.168, 14.684, 16.919, 19.023, 21.666, 23.589],
+          [2.156, 2.558, 3.247, 3.940, 4.865, 15.987, 18.307, 20.483, 23.209, 25.188],
+          [2.603, 3.053, 3.816, 4.575, 5.578, 17.275, 19.675, 21.920, 24.725, 26.757],
+          [3.074, 3.571, 4.404, 5.226, 6.304, 18.549, 21.026, 23.337, 26.217, 28.300],
+          [3.565, 4.107, 5.009, 5.892, 7.042, 19.812, 22.362, 24.736, 27.688, 29.819],
+          [4.075, 4.660, 5.629, 6.571, 7.790, 21.064, 23.685, 26.119, 29.141, 31.319],
+          [4.601, 5.229, 6.262, 7.261, 8.547, 22.307, 24.996, 27.488, 30.578, 32.801],
+          [5.142, 5.812, 6.908, 7.962, 9.312, 23.542, 26.296, 28.845, 32.000, 34.267],
+          [5.697, 6.408, 7.564, 8.672, 10.085, 24.769, 27.587, 30.191, 33.409, 35.718],
+          [6.265, 7.015, 8.231, 9.390, 10.865, 25.989, 28.869, 31.526, 34.805, 37.156],
+          [6.844, 7.633, 8.907, 10.117, 11.651, 27.204, 30.144, 32.852, 36.191, 38.582],
+          [7.434, 8.260, 9.591, 10.851, 12.443, 28.412, 31.410, 34.170, 37.566, 39.997],
+          [8.034, 8.897, 10.283, 11.591, 13.240, 29.615, 32.671, 35.479, 38.932, 41.401],
+          [8.643, 9.542, 10.982, 12.338, 14.041, 30.813, 33.924, 36.781, 40.289, 42.796],
+          [9.260, 10.196, 11.689, 13.091, 14.848, 32.007, 35.172, 38.076, 41.638, 44.181],
+          [9.886, 10.856, 12.401, 13.848, 15.659, 33.196, 36.415, 39.364, 42.980, 45.559],
+          [10.520, 11.524, 13.120, 14.611, 16.473, 34.382, 37.652, 40.646, 44.314, 46.928],
+          [11.160, 12.198, 13.844, 15.379, 17.292, 35.563, 38.885, 41.923, 45.642, 48.290],
+          [11.808, 12.879, 14.573, 16.151, 18.114, 36.741, 40.113, 43.195, 46.963, 49.645],
+          [12.461, 13.565, 15.308, 16.928, 18.939, 37.916, 41.337, 44.461, 48.278, 50.993],
+          [13.121, 14.256, 16.047, 17.708, 19.768, 39.087, 42.557, 45.722, 49.588, 52.336],
+          [13.787, 14.953, 16.791, 18.493, 20.599, 40.256, 43.773, 46.979, 50.892, 53.672],
+          [20.707, 22.164, 24.433, 26.509, 29.051, 51.805, 55.758, 59.342, 63.691, 66.766],
+          [27.991, 29.707, 32.357, 34.764, 37.689, 63.167, 67.505, 71.420, 76.154, 79.490],
+          [35.534, 37.485, 40.482, 43.188, 46.459, 74.397, 79.082, 83.298, 88.379, 91.952],
+          [43.275, 45.442, 48.758, 51.739, 55.329, 85.527, 90.531, 95.023, 100.425, 104.215],
+          [51.172, 53.540, 57.153, 60.391, 64.278, 96.578, 101.879, 106.629, 112.329, 116.321],
+          [59.196, 61.754, 65.647, 69.126, 73.291, 107.565, 113.145, 118.136, 124.116, 128.299],
+          [67.328, 70.065, 74.222, 77.929, 82.358, 118.498, 124.342, 129.561, 135.807, 140.169]
+        ].freeze
+
+        ALPHA_HEADER =
+          {
+            0.995 => 0,
+            0.990 => 1,
+            0.975 => 2,
+            0.95 => 3,
+            0.9 => 4,
+            0.1 => 5,
+            0.05 => 6,
+            0.025 => 7,
+            0.01 => 8,
+            0.005 => 9
+          }.freeze
+
+        DEGREES_OF_FREEDOM = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 40, 50, 60, 70, 80, 90, 100].freeze
+
+        # return an array of the alpha values in correct order
+        # i.e. 0.995 -> 0, 0.990 -> 1 etc
+        def self.alpha_values
+          ALPHA_HEADER.keys
+        end
+
+        # Checks if a valid alpha value is passed to look up values
+        def self.valid_alpha?(alpha)
+          self.alpha_values.include?(alpha)
+        end
+
+        # Return a whole column of the distribution table for a certain alpha value
+        def self.alpha_column(alpha)
+          raise "Undefined alpha value." unless self.valid_alpha?(alpha)
+
+          # return an array of hashes for an alpha value
+          # including the degree of freedom for each critical value
+          # i.e. [{df: 1, critival_value: x},...]
+          TABLE.map.with_index do |row, index|
+            { df: DEGREES_OF_FREEDOM[index], critical_value: row[ALPHA_HEADER[alpha]] }
+          end
+        end
+      end
+    end
+  end
+end

ruby-statistics.gemspec

@@ -31,7 +31,7 @@ Gem::Specification.new do |spec|
   spec.required_ruby_version = '>= 3.0'
 
   spec.add_development_dependency "rake", '~> 13.0', '>= 12.0.0'
-  spec.add_development_dependency "rspec", '~> 3.6', '>= 3.6.0'
+  spec.add_development_dependency "rspec", '~> 3.13', '>= 3.10.0'
   spec.add_development_dependency 'byebug', '~> 11.1', '>= 11.1.0'
   spec.add_development_dependency 'pry', '~> 0.14', '>= 0.14.0'
   spec.add_development_dependency 'bigdecimal', '~> 3.1', '>= 3.1.9'