Start work on Chi^2 distribution.

src/BBS-Numerical-Statistics.adb

@@ -1,17 +1,20 @@
 with Ada.Text_IO;
 with Ada.Numerics;
 with Ada.Numerics.Generic_Elementary_Functions;
+with BBS.Numerical.functions_real;
 with BBS.Numerical.Integration_real;
 package body BBS.Numerical.Statistics is
  package elem is new Ada.Numerics.Generic_Elementary_Functions(f'Base);
+ use type f'Base;
  package float_io is new Ada.Text_IO.Float_IO(f'Base);
  package integ is new BBS.Numerical.Integration_real(f'Base);
+ package funct is new BBS.Numerical.functions_real(f'Base);
  -- -----------------------------------------------------------------
  --
  -- Compute the mean of an array of data
  --
  function mean(d : data_array) return F'Base is
- sum : F := 0.0;
+ sum : F'Base := 0.0;
  begin
  for i in d'Range loop
  sum := sum + d(i);
@@ -39,8 +42,8 @@ package body BBS.Numerical.Statistics is
  -- instead of mean() if you need both values.
  --
  procedure variance(d : data_array; var : out F'Base; mean : out F'Base) is
- sum2 : F := 0.0;
- sum : F := 0.0;
+ sum2 : F'Base := 0.0;
+ sum : F'Base := 0.0;
  begin
  for i in d'Range loop
  sum := sum + d(i);
@@ -54,15 +57,15 @@ package body BBS.Numerical.Statistics is
  --
  -- Normal distribution
  --
- function normal(p : F'Base) return F'Base is
+ function normal_pdf(p : F'Base) return F'Base is
  scale : constant F'Base := 1.0/elem.Sqrt(2.0*Ada.Numerics.Pi);
  begin
  return scale*elem.Exp(-p*p/2.0);
  end;
  --
  -- Normal distribution with mean and sigma
  --
- function normal(p, mean, sigma : F'Base) return F'Base is
+ function normal_pdf(p, mean, sigma : F'Base) return F'Base is
  scale : constant F'Base := 1.0/(sigma*elem.Sqrt(2.0*Ada.Numerics.Pi));
  begin
  return scale*elem.Exp(-(p-mean)*(p-mean)/(2.0*sigma*sigma));
@@ -71,8 +74,19 @@ package body BBS.Numerical.Statistics is
  -- Compute the area under a Normal curve from a to b using the
  -- specified number of steps of Simpson's integration.
  --
- function normal_area(a, b : F'Base; steps : Positive) return F'Base is
+ function normal_cdf(a, b : F'Base; steps : Positive) return F'Base is
  begin
- return integ.simpson(normal'Access, a, b, steps);
+ return integ.simpson(normal_pdf'Access, a, b, steps);
+ end;
+ --
+ -- Chi squared distribution.
+ -- Note that the degrees of freedom (k) must be positive otherwise that
+ -- would imply zero or fewer data points.
+ --
+ -- The value is (x^(k/2)*exp(-x/2))/(2^(k/2)*gamma(k/2))
+ --
+ function chi2_pdf(x : f'Base; k : Positive) return f'Base is
+ begin
+ return (elem."**"(x, (f'Base(k)/2.0))*elem.Exp(-x/2.0))/(elem."**"(2.0, (f'Base(k)/2.0)*funct.gamma2n(k)));
  end;
 end;

src/BBS-Numerical-Statistics.ads

@@ -22,14 +22,23 @@ package BBS.Numerical.statistics is
  --
  -- Standard Normal distribution
  --
- function normal(p : F'Base) return F'Base;
+ function normal_pdf(p : F'Base) return F'Base;
  --
  -- Normal distribution with mean and sigma
  --
- function normal(p, mean, sigma : F'Base) return F'Base;
+ function normal_pdf(p, mean, sigma : F'Base) return F'Base;
  --
  -- Compute the area under a Normal curve from a to b using the
- -- specified number of steps of Simpson's integration.
+ -- specified number of steps of Simpson's integration. Note that the
+ -- normal PDF is symetrical so the value from 0 to infinity is exactly
+ -- 0.5. Calculation of "a" to infinity can be done by taking 0.5
+ -- minus the value from 0 to "a".
  --
- function normal_area(a, b : F'Base; steps : Positive) return F'Base;
+ function normal_cdf(a, b : F'Base; steps : Positive) return F'Base;
+ --
+ -- Chi square distribution.
+ -- Note that the degrees of freedom (k) must be positive otherwise that
+ -- would imply zero or fewer data points.
+ --
+ function chi2_pdf(x : f'Base; k : Positive) return f'Base;
 end;

src/BBS-Numerical-functions_real.adb

@@ -0,0 +1,37 @@
+with Ada.Text_IO;
+with Ada.Numerics;
+with Ada.Numerics.Generic_Elementary_Functions;
+package body BBS.Numerical.functions_real is
+ package float_io is new Ada.Text_IO.Float_IO(f'Base);
+ package elem is new Ada.Numerics.Generic_Elementary_Functions(f'Base);
+ --
+ -- Compute the gamma function of a positive number divided by two.
+ -- This is initially used by the chi-squared statistics function, but
+ -- may have other applications. The full gamma function will be
+ -- implemented when needed.
+ --
+ -- Note that only gamma of positive and half positive integers are needed,
+ -- i.e. 1/2, 1, 3/2, 2, 5/2, ...
+ -- This will be easier to implement than the full gamma functions since
+ -- gamma(1/2) is sqrt(pi), gamma(1) is 1, and gamma(a+1) is a*gamma(a).
+ --
+ function gamma2n(n : Positive) return f'Base is
+ base : f'Base := elem.sqrt(Ada.Numerics.Pi);
+ begin
+ --
+ -- If n is even, then base is 1, otherwise it is sqrt(pi)
+ --
+ if 2*(n/2) = n then
+ base := 1.0;
+ for i in 2 .. (n/2) loop
+ base := base*f'Base(i-1);
+ end loop;
+ else
+ base := elem.sqrt(Ada.Numerics.Pi);
+ for i in 1 .. (n/2) loop
+ base := base*(f'Base(i)-0.5);
+ end loop;
+ end if;
+ return base;
+ end;
+end BBS.Numerical.functions_real;

src/BBS-Numerical-functions_real.ads

@@ -0,0 +1,20 @@
+--
+-- This is a collection of functions that are used by other numerical
+-- but aren't provided by the Ada elementary functions package. Functions
+-- that are more specific to a sub group would be in that package. For
+-- example, the normal distribution function is in the Statistics package.
+generic
+ type F is digits <>;
+package BBS.Numerical.functions_real is
+ -- Compute the gamma function of a positive number divided by two.
+ -- This is initially used by the chi-squared statistics function, but
+ -- may have other applications. The full gamma function will be
+ -- implemented when needed.
+ --
+ -- Note that only gamma of positive and half positive integers are needed,
+ -- i.e. 1/2, 1, 3/2, 2, 5/2, ...
+ -- This will be easier to implement than the full gamma functions since
+ -- gamma(1/2) is sqrt(pi), gamma(1) is 1, and gamma(a+1) is a*gamma(a).
+ --
+ function gamma2n(n : Positive) return f'Base;
+end BBS.Numerical.functions_real;

test/test_stats.adb

@@ -3,13 +3,15 @@ with Ada.Text_IO;
 with BBS.Numerical;
 with BBS.Numerical.regression;
 with BBS.Numerical.Statistics;
+with BBS.Numerical.functions_real;
 
 procedure test_stats is
  subtype real is Float;
  package linreg is new BBS.Numerical.regression(real);
  package stat is new BBS.Numerical.Statistics(real);
  package float_io is new Ada.Text_IO.Float_IO(real);
  package elem is new Ada.Numerics.Generic_Elementary_Functions(real);
+ package funct is new BBS.Numerical.functions_real(real);
 
  data : linreg.data_array :=
  ((x => 1.0, y => 1.0),
@@ -32,15 +34,19 @@ begin
  Ada.Text_IO.Put(", variance = ");
  float_io.Put(res.variance, 2, 3, 0);
  Ada.Text_IO.New_Line;
- Ada.Text_IO.Put_Line("Normal Curve Areas");
- for i in 0 .. 50 loop
- val := real(i)*0.1;
+ Ada.Text_IO.Put_Line("Probability Distributions");
+ Ada.Text_IO.Put_Line(" Normal Chi^2");
+ Ada.Text_IO.Put_Line(" X PDF CDF PDF");
+ for i in 0 .. 20 loop
+ val := real(i)*0.01;
  Ada.Text_IO.Put(" ");
  float_io.Put(val, 1, 2, 0);
  Ada.Text_IO.Put(" ");
- float_io.Put(stat.normal(val), 1, 5, 0);
+ float_io.Put(stat.normal_pdf(val), 1, 5, 0);
  Ada.Text_IO.Put(" ");
- float_io.Put(stat.normal_area(0.0, val, 200), 1, 5, 0);
+ 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.New_Line;
  end loop;
 end test_stats;