ref: 7ae7351ace73a70cbbaf5a2da016005485d0d43b
dir: /appl/math/fit.b/
# fit a polynomial to a set of points
# fit -dn [-v]
# where n is the degree of the polynomial
implement Fit;
include "sys.m";
sys: Sys;
include "draw.m";
include "math.m";
maths: Math;
include "bufio.m";
bufio: Bufio;
include "arg.m";
Fit: module
{
init: fn(nil: ref Draw->Context, argv: list of string);
};
MAXPTS: con 512;
MAXDEG: con 16;
EPS: con 0.0000005;
init(nil: ref Draw->Context, argv: list of string)
{
sys = load Sys Sys->PATH;
maths = load Math Math->PATH;
if(maths == nil)
fatal(sys->sprint("cannot load maths library"));
bufio = load Bufio Bufio->PATH;
if(bufio == nil)
fatal(sys->sprint("cannot load bufio"));
main(argv);
}
isn(r: real, n: int): int
{
s := r - real n;
if(s < 0.0)
s = -s;
return s < EPS;
}
fact(n: int): real
{
f := 1.0;
for(i := 1; i <= n; i++)
f *= real i;
return f;
}
comb(n: int, r: int): real
{
f := 1.0;
for(i := 0; i < r; i++)
f *= real (n-i);
return f/fact(r);
}
matalloc(n: int): array of array of real
{
mat := array[n] of array of real;
for(i := 0; i < n; i++)
mat[i] = array[n] of real;
return mat;
}
matsalloc(n: int): array of array of array of real
{
mats := array[n+1] of array of array of real;
for(i := 0; i <= n; i++)
mats[i] = matalloc(i);
return mats;
}
det(mat: array of array of real, n: int, mats: array of array of array of real): real
{
# easy cases first
if(n == 0)
return 1.0;
if(n == 1)
return mat[0][0];
if(n == 2)
return mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
d := 0.0;
s := 1;
m := mats[n-1];
for(k := 0; k < n; k++){
for(i := 0; i < n-1; i++){
for(j := 0; j < n-1; j++){
if(j < k)
m[i][j] = mat[i+1][j];
else
m[i][j] = mat[i+1][j+1];
}
}
d += (real s)*mat[0][k]*det(m, n-1, mats);
s = -s;
}
return d;
}
main(argv: list of string)
{
i, j: int;
x, y, z: real;
fb: ref Bufio->Iobuf;
n := 0;
p := 1;
arg := load Arg Arg->PATH;
if(arg == nil)
fatal(sys->sprint("cannot load %s: %r", Arg->PATH));
arg->init(argv);
verbose := 0;
while((o := arg->opt()) != 0)
case o{
'd' =>
p = int arg->arg();
'v' =>
verbose = 1;
* =>
fatal(sys->sprint("bad option %c", o));
}
args := arg->argv();
arg = nil;
if(args != nil){
s := hd args;
fb = bufio->open(s, bufio->OREAD);
if(fb == nil)
fatal(sys->sprint("cannot open %s", s));
}
else{
fb = bufio->open("/dev/cons", bufio->OREAD);
if(fb == nil)
fatal(sys->sprint("missing data file name"));
}
a := array[p+1] of real;
b := array[p+1] of real;
sx := array[2*p+1] of real;
sxy := array[p+1] of real;
xd := array[MAXPTS] of real;
yd := array[MAXPTS] of real;
while(1){
xs := ss(bufio->fb.gett(" \t\r\n"));
if(xs == nil)
break;
ys := ss(bufio->fb.gett(" \t\r\n"));
if(ys == nil)
fatal(sys->sprint("missing value"));
if(n >= MAXPTS)
fatal(sys->sprint("too many points"));
xd[n] = real xs;
yd[n] = real ys;
n++;
}
if(p < 0)
fatal(sys->sprint("negative power"));
if(p > MAXDEG)
fatal(sys->sprint("power too large"));
if(n < p+1)
fatal(sys->sprint("not enough points"));
# use x-xbar, y-ybar to avoid overflow
for(i = 0; i <= p; i++)
sxy[i] = 0.0;
for(i = 0; i <= 2*p; i++)
sx[i] = 0.0;
xbar := ybar := 0.0;
for(i = 0; i < n; i++){
xbar += xd[i];
ybar += yd[i];
}
xbar = xbar/(real n);
ybar = ybar/(real n);
for(i = 0; i < n; i++){
x = xd[i]-xbar;
y = yd[i]-ybar;
for(j = 0; j <= p; j++)
sxy[j] += y*x**j;
for(j = 0; j <= 2*p; j++)
sx[j] += x**j;
}
mats := matsalloc(p+1);
mat := mats[p+1];
for(i = 0; i <= p; i++)
for(j = 0; j <= p; j++)
mat[i][j] = sx[i+j];
d := det(mat, p+1, mats);
if(isn(d, 0))
fatal(sys->sprint("points not independent"));
for(j = 0; j <= p; j++){
for(i = 0; i <= p; i++)
mat[i][j] = sxy[i];
a[j] = det(mat, p+1, mats)/d;
for(i = 0; i <= p; i++)
mat[i][j] = sx[i+j];
}
if(verbose)
sys->print("\npt actual x actual y predicted y\n");
e := 0.0;
for(i = 0; i < n; i++){
x = xd[i]-xbar;
y = yd[i]-ybar;
z = 0.0;
for(j = 0; j <= p; j++)
z += a[j]*x**j;
z += ybar;
e += (z-yd[i])*(z-yd[i]);
if(verbose)
sys->print("%d. %f %f %f\n", i+1, xd[i], yd[i], z);
}
if(verbose)
sys->print("root mean squared error = %f\n", maths->sqrt(e/(real n)));
for(i = 0; i <= p; i++)
b[i] = 0.0;
b[0] += ybar;
for(i = 0; i <= p; i++)
for(j = 0; j <= i; j++)
b[j] += a[i]*comb(i, j)*(-xbar)**(i-j);
pr := 0;
sys->print("y = ");
for(i = p; i >= 0; i--){
if(!isn(b[i], 0) || (i == 0 && pr == 0)){
if(b[i] < 0.0){
sys->print("-");
b[i] = -b[i];
}
else if(pr)
sys->print("+");
pr = 1;
if(i == 0)
sys->print("%f", b[i]);
else{
if(!isn(b[i], 1))
sys->print("%f*", b[i]);
sys->print("x");
if(i > 1)
sys->print("^%d", i);
}
}
}
sys->print("\n");
}
ss(s: string): string
{
l := len s;
while(l > 0 && (s[0] == ' ' || s[0] == '\t' || s[0] == '\r' || s[0] == '\n')){
s = s[1: ];
l--;
}
while(l > 0 && (s[l-1] == ' ' || s[l-1] == '\t' || s[l-1] == '\r' || s[l-1] == '\n')){
s = s[0: l-1];
l--;
}
return s;
}
fatal(s: string)
{
sys->fprint(sys->fildes(2), "fit: %s\n", s);
exit;
}