Updates on Chi squared and Lower incomplete gamma function (#122)

* feat: Calculate CDF for Xi square using a special case when df = 2.

This change introduces the calculation of the CDF for the xi square
distribution by executing the special case when the degrees of
freedom is 2. This is to make it on par with calculations made in R
and some python packages. It derived from work made by @oliver-czulo
in #115.

* fix: Always use at least 10_000 iterations when calculating lower gamma.

This change fixes a small bug where we were calculating the
lower incomplete gamma function via simpson's rule with lower
than 10_000 iterations, which is not ideal. This was found while
doing some debugging on #115.

lib/math.rb

@@ -45,8 +45,11 @@ def self.simpson_rule(a, b, n, &block)
  end
 
  def self.lower_incomplete_gamma_function(s, x)
+ base_iterator = x.round(1)
+ base_iterator += 1 if x < 1.0 && !x.zero?
+
  # The greater the iterations, the better. That's why we are iterating 10_000 * x times
- iterator = (10_000 * x.round(1)).round
+ iterator = (10_000 * base_iterator).round
  iterator = 100_000 if iterator.zero?
 
  self.simpson_rule(0, x.to_r, iterator) do |t|

lib/statistics/distribution/chi_squared.rb

@@ -10,8 +10,13 @@ def initialize(k)
  end
 
  def cumulative_function(value)
- k = degrees_of_freedom/2.0
- Math.lower_incomplete_gamma_function(k, value/2.0)/Math.gamma(k)
+ if degrees_of_freedom == 2
+ # Special case where DF = 2 https://en.wikipedia.org/wiki/Chi-squared_distribution#Cumulative_distribution_function
+ 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)
+ end
  end
 
  def density_function(value)

spec/statistics/distribution/chi_squared_spec.rb

@@ -10,6 +10,20 @@
  expect(chi_sq.cumulative_function(number).round(4)).to eq results[index]
  end
  end
+
+ context 'when the degrees of freedom is 2 for the chi-squared distribution' do
+ it 'performs the probability calculation using the special case instead' do
+ # Results gathered from R 4.0.3
+ # # > pchisq(1:10, df = 2)
+ # # [1] 0.3935 0.6321 0.7769 0.8647 0.9179 0.9502 0.9698 0.9817 0.9889 0.9933
+ results = [0.3935, 0.6321, 0.7769, 0.8647, 0.9179, 0.9502, 0.9698, 0.9817, 0.9889, 0.9933]
+ chi_sq = described_class.new(2)
+
+ (1..10).each_with_index do |number, index|
+ expect(chi_sq.cumulative_function(number).round(4)).to eq results[index]
+ end
+ end
+ end
  end
 
  describe '#density_function' do

spec/statistics/math_spec.rb

@@ -76,6 +76,16 @@
  described_class.lower_incomplete_gamma_function(lower, upper)
  end
 
+ it "uses the simpson's rule with iterations bigger than 10_000 when upper is 0.0 < x < 1.0" do
+ lower = 0
+ upper = rand(0.01..0.99)
+ iteration = 10_000 * (1 + upper.round(1))
+
+ expect(described_class).to receive(:simpson_rule).with(lower, upper.to_r, iteration)
+
+ described_class.lower_incomplete_gamma_function(lower, upper)
+ end
+
  it 'returns the expected calculation' do
  results = [0.6322, 0.594, 1.1536, 3.3992, 13.4283]