Add simple linear regression and some testing.

numerical.gpr → BBS_Numerical.gpr

@@ -1,6 +1,6 @@
 with "../BBS-Ada/bbs.gpr";
 
-library project Numerical is
+library project BBS_Numerical is
 
  for Languages use ("Ada");
  for Library_Name use "Numerical";
@@ -9,5 +9,5 @@ library project Numerical is
  for Library_Dir use "lib";
  for Library_Kind use "Static";
 
-end Numerical;
+end BBS_Numerical;
 

README.md

@@ -30,3 +30,7 @@ The implemented algorithms are:
 * Mueller
 
 ## Quaternions
+Some basic operations are implemented.
+
+## Regression
+Simple linear regression is implemented.

src/BBS-Numerical-regression.adb

@@ -0,0 +1,36 @@
+with Ada.Text_IO;
+package body BBS.Numerical.regression is
+ package float_io is new Ada.Text_IO.Float_IO(f);
+ --
+ -- Compute the simple linear regression to an array of points.
+ --
+ function simple_linear(d : data_array) return simple_linreg_result is
+ Sx : f'Base := 0.0; -- Sum of x
+ Sy : f'Base := 0.0; -- Sum of y
+ Sxy : f'Base := 0.0; -- Sum of x*y
+ Sxx : f'Base := 0.0; -- Sum of x*x
+ SSxy : f'Base := 0.0; -- (Sum of x*y) squared
+ SSxx : f'Base := 0.0; -- (Sum of x*x) squared
+ err : f'Base; -- Errors
+ SSe : f'Base := 0.0; -- Sum of err*err
+ a : f'Base; -- Y intercept
+ b : f'Base; -- Slope
+ n : constant f'Base := f'Base(d'Length);
+ begin
+ for i in d'Range loop
+ Sx := Sx + d(i).x;
+ Sy := Sy + d(i).y;
+ Sxy := Sxy + d(i).x*d(i).y;
+ Sxx := Sxx + d(i).x*d(i).x;
+ end loop;
+ SSxy := Sxy - (Sx*Sy)/n;
+ SSxx := Sxx - (Sx*Sx)/n;
+ a := Sy/n - (Sx/n)*(SSxy/SSxx);
+ b := SSxy/SSxx;
+ for i in d'Range loop
+ err := d(i).y - (a + b*d(i).x);
+ SSe := SSe + err*err;
+ end loop;
+ return (a => a, b => b, SSe => SSe, variance => SSe/(n-2.0));
+ end;
+end BBS.Numerical.regression;

src/BBS-Numerical-regression.ads

@@ -0,0 +1,23 @@
+generic
+ type F is digits <>;
+package BBS.Numerical.regression is
+
+ type point is record
+ x : f'Base;
+ y : f'Base;
+ end record;
+ type data_array is array (Integer range <>) of point;
+
+ --
+ -- Holds the return data from the simple linear regression function.
+ -- It contains coefficients for the line y = a + bX.
+ type simple_linreg_result is record
+ a : f'Base; -- Y intercept
+ b : f'Base; -- Slope
+ SSe : f'Base; -- Sum of square errors
+ variance : f'Base; -- Variance
+ end record;
+
+ function simple_linear(d : data_array) return simple_linreg_result;
+
+end BBS.Numerical.regression;

test/test.adb

@@ -0,0 +1,33 @@
+with Ada.Text_IO;
+with BBS.Numerical;
+with BBS.Numerical.complex_real;
+with BBS.Numerical.integration_real;
+with BBS.Numerical.regression;
+with BBS.Numerical.roots_real;
+with BBS.Numerical.Statistics;
+
+procedure test is
+ package linreg is new BBS.Numerical.regression(Float);
+ package float_io is new Ada.Text_IO.Float_IO(Float);
+
+ data : linreg.data_array :=
+ ((x => 1.0, y => 1.0),
+ (x => 2.0, y => 1.0),
+ (x => 3.0, y => 2.0),
+ (x => 4.0, y => 2.0),
+ (x => 5.0, y => 4.0));
+ res : linreg.simple_linreg_result;
+begin
+ Ada.Text_IO.Put_Line("Testing some of the numerical routines.");
+ res := linreg.simple_linear(data);
+ Ada.Text_IO.Put_Line("Linear regression result:");
+ Ada.Text_IO.Put("a = ");
+ float_io.Put(res.a, 2, 3, 0);
+ Ada.Text_IO.Put(", b = ");
+ float_io.Put(res.b, 2, 3, 0);
+ Ada.Text_IO.Put(", SSe = ");
+ float_io.Put(res.SSe, 2, 3, 0);
+ Ada.Text_IO.Put(", variance = ");
+ float_io.Put(res.variance, 2, 3, 0);
+ Ada.Text_IO.New_Line;
+end test;

test/test.gpr

@@ -0,0 +1,12 @@
+with "../BBS_Numerical.gpr";
+
+project test is
+ for languages use ("Ada");
+ for Source_Dirs use (".");
+ for Object_Dir use "../obj";
+ for Main use ("test.adb");
+ for Exec_Dir use ".";
+ package compiler is
+ for switches ("Ada") use ("-g", "-gnateE");
+ end compiler;
+end test;