Move work on chi^2 distribution. Still isn't quite right.

src/BBS-Numerical-Statistics.adb

@@ -84,9 +84,15 @@ package body BBS.Numerical.Statistics is
  -- would imply zero or fewer data points.
  --
  -- (x^(k/2-1)*exp(-x/2))
- -- The value us --------------------
+ -- The value is --------------------
  -- (2^(k/2)*gamma(k/2))
  --
+ -- For k > about 70, gamma(k/2) will overflow Float. So for more
+ -- degrees of freedom, need to use logarithms to do the calculations.
+ -- Other large numbers from the exponents will also mostly cancel out.
+ --
+ -- The value is exp[(log(x)*(k/2-1) - x/2) - (log(2)*(k/2) + log(gamma(k/2)))]
+ --
  function chi2_pdf(x : f'Base; k : Positive) return f'Base is
  begin
  if x > 0.0 then
@@ -97,6 +103,19 @@ package body BBS.Numerical.Statistics is
  end if;
  end;
  --
+ function chi2_exp(x : f'Base; k : Positive) return f'Base is
+ part1 : f'Base;
+ part2 : f'Base;
+ begin
+ if x > 0.0 then
+ part1 := elem.log(x)*((f'Base(k)/2.0)-1.0) - (x/2.0);
+ part2 := elem.log(2.0)*(f'Base(x)/2.0) + funct.lngamma2n(k);
+ return elem.exp(part1 - part2);
+ else
+ return 0.0;
+ end if;
+ end;
+ --
  function chi2_cdf(a, b : F'Base; k, steps : Positive) return F'Base is
  begin
  deg_freedom := k;
@@ -105,6 +124,6 @@ package body BBS.Numerical.Statistics is
  --
  function partial_chi2(x : f'Base) return f'Base is
  begin
- return chi2_pdf(x, deg_freedom);
+ return chi2_exp(x, deg_freedom);
  end;
 end;

src/BBS-Numerical-Statistics.ads

@@ -41,6 +41,7 @@ package BBS.Numerical.statistics is
  -- would imply zero or fewer data points.
  --
  function chi2_pdf(x : f'Base; k : Positive) return f'Base;
+ function chi2_exp(x : f'Base; k : Positive) return f'Base;
  --
  function chi2_cdf(a, b : F'Base; k, steps : Positive) return F'Base;
  --

test/test_stats.adb

@@ -38,16 +38,18 @@ begin
  Ada.Text_IO.Put_Line(" Normal Chi^2");
  Ada.Text_IO.Put_Line(" X PDF CDF PDF CDF");
  for i in 0 .. 20 loop
- val := real(i)*0.1;
+ val := real(i)*1.0;
  Ada.Text_IO.Put(" ");
- float_io.Put(val, 1, 2, 0);
+ float_io.Put(val, 2, 2, 0);
  Ada.Text_IO.Put(" ");
  float_io.Put(stat.normal_pdf(val), 1, 5, 0);
  Ada.Text_IO.Put(" ");
  float_io.Put(stat.normal_cdf(0.0, val, 20), 1, 5, 0);
  Ada.Text_IO.Put(" ");
  float_io.Put(stat.chi2_pdf(val, 1), 1, 5, 0);
  Ada.Text_IO.Put(" ");
+ float_io.Put(stat.chi2_exp(val, 1), 1, 5, 0);
+ Ada.Text_IO.Put(" ");
  float_io.Put(stat.chi2_cdf(0.0, val, 1, 20), 1, 5, 0);
  Ada.Text_IO.New_Line;
  end loop;