Add normal distribution and area under curve to statistics package.

src/BBS-Numerical-Statistics.adb

@@ -1,11 +1,16 @@
---with Ada.Text_IO;
---with Ada.Float_Text_IO;
+with Ada.Text_IO;
+with Ada.Numerics;
 with Ada.Numerics.Generic_Elementary_Functions;
+with BBS.Numerical.Integration_real;
 package body BBS.Numerical.Statistics is
+ package elem is new Ada.Numerics.Generic_Elementary_Functions(f'Base);
+ package float_io is new Ada.Text_IO.Float_IO(f'Base);
+ package integ is new BBS.Numerical.Integration_real(f'Base);
+ -- -----------------------------------------------------------------
  --
  -- Compute the mean of an array of data
  --
- function mean(d : data_array) return F is
+ function mean(d : data_array) return F'Base is
  sum : F := 0.0;
  begin
  for i in d'Range loop
@@ -16,7 +21,7 @@ package body BBS.Numerical.Statistics is
  --
  -- Compute the limits of an array of data
  --
- procedure limits(d : data_array; min : out F; max : out F) is
+ procedure limits(d : data_array; min : out F'Base; max : out F'Base) is
  begin
  min := d(d'First);
  max := d(d'First);
@@ -33,7 +38,7 @@ package body BBS.Numerical.Statistics is
  -- Compute the variance (and mean) of an array of data. Use this,
  -- instead of mean() if you need both values.
  --
- procedure variance(d : data_array; var : out F; mean : out F) is
+ procedure variance(d : data_array; var : out F'Base; mean : out F'Base) is
  sum2 : F := 0.0;
  sum : F := 0.0;
  begin
@@ -42,6 +47,32 @@ package body BBS.Numerical.Statistics is
  sum2 := sum2 + (d(i)*d(i));
  end loop;
  mean := sum/F(d'Length);
- var := (sum2 - sum*sum/F(d'Length))/F(d'Length - 1);
+ var := (sum2 - sum*sum/F'Base(d'Length))/F'Base(d'Length - 1);
+ end;
+ -- -----------------------------------------------------------------
+ -- Distributions
+ --
+ -- Normal distribution
+ --
+ function normal(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
+ 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));
+ end;
+ --
+ -- 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
+ begin
+ return integ.simpson(normal'Access, a, b, steps);
  end;
 end;

src/BBS-Numerical-Statistics.ads

@@ -1,18 +1,35 @@
 generic
  type F is digits <>;
 package BBS.Numerical.statistics is
- type data_array is array (Integer range <>) of F;
+ type data_array is array (Integer range <>) of F'Base;
+ -- -----------------------------------------------------------------
+ -- Statistics of a sample of data
  --
  -- Compute the mean of an array of data
  --
- function mean(d : data_array) return F;
+ function mean(d : data_array) return F'Base;
  --
  -- Compute the median of an array of data
  --
- procedure limits(d : data_array; min : out F; max : out F);
+ procedure limits(d : data_array; min : out F'Base; max : out F'Base);
  --
  -- Compute the variance (and mean) of an array of data. Use this,
  -- instead of mean() if you need both values.
  --
- procedure variance(d : data_array; var : out F; mean : out F);
+ procedure variance(d : data_array; var : out F'Base; mean : out F'Base);
+ -- -----------------------------------------------------------------
+ -- Distributions
+ --
+ -- Standard Normal distribution
+ --
+ function normal(p : F'Base) return F'Base;
+ --
+ -- Normal distribution with mean and sigma
+ --
+ function normal(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.
+ --
+ function normal_area(a, b : F'Base; steps : Positive) return F'Base;
 end;

src/BBS-Numerical-integration_real.adb

@@ -27,17 +27,17 @@ package body BBS.Numerical.integration_real is
  -- segment, the effective number of steps is doubled.
  --
  function simpson(test : test_func; lower, upper : f'Base; steps : Positive) return f'Base is
- step_size : constant f'Base := (upper - lower)/f(2*steps);
+ step_size : constant f'Base := (upper - lower)/f'Base(2*steps);
  ends : constant f'Base := test(lower) + test(upper);
  sum1 : f'Base := 0.0;
  sum2 : f'Base := 0.0;
  base : f'Base := lower;
  begin
  for i in 1 .. steps - 1 loop
  base := base + step_size;
- sum1 := sum1 + test(base);
+ sum1 := sum1 + test(base); -- Odd steps
  base := base + step_size;
- sum2 := sum2+ test(base);
+ sum2 := sum2 + test(base); -- Even steps
  end loop;
  base := base + step_size;
  sum1 := sum1 + test(base);

test/test_stats.adb

@@ -7,6 +7,7 @@ with BBS.Numerical.Statistics;
 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);
 
@@ -17,6 +18,7 @@ procedure test_stats is
  (x => 4.0, y => 2.0),
  (x => 5.0, y => 4.0));
  res : linreg.simple_linreg_result;
+ val : real;
 begin
  Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
  res := linreg.simple_linear(data);
@@ -30,4 +32,15 @@ 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(" ");
+ float_io.Put(val, 1, 2, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(stat.normal(val), 1, 5, 0);
+ Ada.Text_IO.Put(" ");
+ float_io.Put(stat.normal_area(0.0, val, 200), 1, 5, 0);
+ Ada.Text_IO.New_Line;
+ end loop;
 end test_stats;