Use adaptive simpson's integration for chi^2 CDF. This produces bette…

…r results. Probably a midpoint method would work well.
BrentSeidel · Jun 5, 2024 · 1ee0a64 · 1ee0a64
1 parent 90f8e99
commit 1ee0a64
Show file tree

Hide file tree

Showing 2 changed files with 5 additions and 4 deletions.
diff --git a/src/BBS-Numerical-Statistics.adb b/src/BBS-Numerical-Statistics.adb
@@ -117,9 +117,10 @@ package body BBS.Numerical.Statistics is
  end;
  --
  function chi2_cdf(a, b : F'Base; k, steps : Positive) return F'Base is
+ tol : f'Base := 0.0001;
  begin
  deg_freedom := k;
- return integ.simpson(partial_chi2'Access, a, b, steps);
+ return integ.adapt_simpson(partial_chi2'Access, a, b, tol, steps);
  end;
  --
  function partial_chi2(x : f'Base) return f'Base is

diff --git a/test/test_stats.adb b/test/test_stats.adb
@@ -46,11 +46,11 @@ begin
  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);
+ float_io.Put(stat.chi2_pdf(val, 20), 1, 5, 0);
  Ada.Text_IO.Put(" ");
- float_io.Put(stat.chi2_exp(val, 1), 1, 5, 0);
+ float_io.Put(stat.chi2_exp(val, 20), 1, 5, 0);
  Ada.Text_IO.Put(" ");
- float_io.Put(stat.chi2_cdf(0.0, val, 1, 20), 1, 5, 0);
+ float_io.Put(stat.chi2_cdf(0.0, val, 20, 20), 1, 5, 0);
  Ada.Text_IO.New_Line;
  end loop;
 end test_stats;