A program to fit data points to polynomials

Ted Stefanik ted at adelie.Adelie.COM
Wed Aug 9 09:00:30 AEST 1989


Here is a program that allows you to fit a set of points to the closest
polynomial.  It also has an imbedded matrix inverter/solver, which I couldn't
find elsewhere on my ULTRIX system.

Note that solving for a polynomial above about the fifth order can easily
cause floating point exceptions.

#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
#	Makefile
#	fit.c
#	fit.h
#	mat.c
# This archive created: Tue Aug  8 18:49:44 1989
export PATH; PATH=/bin:/usr/bin:$PATH
echo shar: "extracting 'Makefile'" '(148 characters)'
if test -f 'Makefile'
then
	echo shar: "will not over-write existing file 'Makefile'"
else
sed 's/^X#//' << \SHAR_EOF > 'Makefile'
X#D      = -g
X#CFLAGS = $(D) $(FLOAT)
X#FLOAT  =
X#MATH   = -lm
X#
X#
X#fit: fit.o mat.o
X#	cc fit.o mat.o $(D) $(FLOAT) -o fit $(MATH)
X#
X#fit.o: fit.h
X#mat.o: fit.h
SHAR_EOF
fi
echo shar: "extracting 'fit.c'" '(3032 characters)'
if test -f 'fit.c'
then
	echo shar: "will not over-write existing file 'fit.c'"
else
sed 's/^X#//' << \SHAR_EOF > 'fit.c'
X##include "fit.h"
X#
X#static void ParseArgs();
X#static void Usage();
X#
X#static void ReadPts();
X#static void PolyFit();
X#static void Results();
X#static void Stats();
X#
X#
X#static points Points;
X#static int    MaxPt;
X#static int    PtCnt;
X#static int    Order;
X#
X#static matrix Mat;
X#static double Det;
X#
X#
X#main(argc, argv)
X#   int  argc;
X#   char **argv;
X#{
X#   ParseArgs(argc, argv);
X#
X#   ReadPts();
X#   if (PtCnt <= 0)
X#      Usage();
X#   PolyFit();
X#   Results();
X#
X#   Stats();
X#}
X#
X#
X#static void ParseArgs(argc, argv)
X#   int argc;
X#   char **argv;
X#{
X#   if (argc != 3)
X#      Usage(argv[0]);
X#
X#   MaxPt = atoi(argv[1]);
X#   Order = atoi(argv[2]);
X#
X#   if (MaxPt <= 0 || Order <= 0)
X#      Usage(argv[0]);
X#}
X#
X#
X#static void Usage(prog)
X#   char *prog;
X#{
X#   fprintf(stderr,
X#           "Usage:  %s max-number-of-data-points order-of-polynomial\n", prog);
X#   fprintf(stderr,
X#           "Note:   points are read from stdin (a coordinate point/line)\n");
X#   exit(1);
X#}
X#
X#
X#static void ReadPts()
X#{
X#   char line[1024];
X#
X#   Points = (points)Calloc(MaxPt, sizeof(point));
X#
X#   for (PtCnt = 0;  PtCnt < MaxPt;  PtCnt++)
X#   {
X#      if (gets(line) == NULL)
X#         return;
X#      sscanf(line, "%f %f", &Points[PtCnt].X, &Points[PtCnt].Y);
X#   }
X#}
X#
X#
X#static void PolyFit()
X#{
X#   register int    i, p, r, c;
X#   register matrow xvec;
X#
X#   Mat = MatInit(Order + 1);
X#   xvec = (matrow)Calloc(Order * 2 + 1, sizeof(matelm));
X#
X#   for (p = 0;  p < PtCnt;  p++)
X#   {
X#      for (xvec[0] = 1.0, i = 1;  i <= Order * 2;  i++)
X#         xvec[i] = xvec[i - 1] * Points[p].X;
X#      for (r = 0;  r < Order + 1;  r++)
X#      {
X#         Mat[r][Order + 1] += Points[p].Y * xvec[r];
X#         for (c = 0;  c < Order + 1;  c++)
X#            Mat[r][c] += xvec[r + c];
X#      }
X#   }
X#
X#   Det = MatInv(Order + 1, Mat);
X#}
X#
X#
X#static void Results()
X#{
X#   register int i;
X#
X#   printf("%5d = Number of data points input\n", PtCnt);
X#   printf("%5d = Order of polynomial fit\n\n",   Order);
X#
X##ifndef NO_DETERM
X#   printf("%e = Determinate of matrix solution\n\n", Det);
X##endif
X#
X#   printf("y = %f", Mat[0][Order + 1]);
X#   for (i = 1;  i < Order + 1;  i++)
X#      printf("  +  %f x^%d", Mat[i][Order + 1], i);
X#}
X#
X#
X#static void Stats()
X#{
X#   register int i, j;
X#   double   x, y;
X#   double   sum, sum2, minerr, maxerr;
X#
X#   sum = sum2 = maxerr = 0;
X#   minerr = HUGE;
X#
X#   for (i = 0;  i < PtCnt;  i++)
X#   {
X#      y = Mat[0][Order + 1];
X#      for (x = 1, j = 1;  j < Order + 1;  j++)
X#      {
X#         x *= Points[i].X;
X#         y += x * Mat[j][Order + 1];
X#      }
X#
X#      y -= Points[i].Y;
X#      if (y < 0)
X#         y = -y;
X#
X#      sum  += y;
X#      sum2 += y * y;
X#
X#      if (y < minerr)
X#         minerr = y;
X#      if (y > maxerr)
X#         maxerr = y;
X#   }
X#
X#   printf("\n\nFit statistics:\n");
X#
X#   printf("%30.15f = Standard Error of Estimate\n", sqrt(sum2 / (PtCnt - 2)));
X#   printf("%30.15f = Average Error\n", sum / PtCnt);
X#   printf("%30.15f = Error's Standard Deviation\n", 
X#             sqrt((PtCnt * sum2 - sum * sum) / PtCnt / (PtCnt - 1)));
X#   printf("%30.15f = Minimum Error\n", minerr);
X#   printf("%30.15f = Maximum Error\n", maxerr);
X#}
SHAR_EOF
fi
echo shar: "extracting 'fit.h'" '(314 characters)'
if test -f 'fit.h'
then
	echo shar: "will not over-write existing file 'fit.h'"
else
sed 's/^X#//' << \SHAR_EOF > 'fit.h'
X##include <stdio.h>
X##include <math.h>
X#
X#struct point
X#{
X#   double X;
X#   double Y;
X#};
X#
X#typedef struct point point;
X#
X#typedef point *points;
X#
X#typedef double matelm;
X#typedef matelm *matrow;
X#typedef matrow *matrix;
X#
X#typedef char   *pointer;
X#
X#
X#extern char *calloc();
X#
X#matrix  MatInit();
X#double  MatInv();
X#pointer Calloc();
SHAR_EOF
fi
echo shar: "extracting 'mat.c'" '(1330 characters)'
if test -f 'mat.c'
then
	echo shar: "will not over-write existing file 'mat.c'"
else
sed 's/^X#//' << \SHAR_EOF > 'mat.c'
X##include "fit.h"
X#
X#matrix MatInit(n)
X#   register int    n;
X#{
X#   register int    i;
X#   register matrix mat;
X#
X#   mat = (matrix)Calloc(n, sizeof(matrow));
X#
X#   for (i = 0;  i < n;  i++)
X#      mat[i] = (matrow)Calloc(n + 1, sizeof(matelm));
X#
X#   return (mat);
X#}
X#
X#
X#double MatInv(n, mat)
X#   register int    n;
X#   register matrix mat;
X#{
X#   int    pivRow;
X#   double piv, det, div, pivTmp;
X#   matrow tmpRow;
X#   register int    p, r, c;
X#
X#   for (det = 1.0, p = 0;  p < n;  p++)
X#   {
X#      for (piv = mat[pivRow = p][p], r = p;  r < n;  r++)
X#         if (((pivTmp = mat[r][p]) < 0 ? -pivTmp : pivTmp) > piv)
X#            piv = mat[pivRow = r][p];
X#
X##ifndef NO_DETERM
X#      if ((det *= piv) == 0.0)
X#         break;
X##endif
X#
X#      if (pivRow != p)
X#      {
X##ifndef NO_DETERM
X#         det = - det;
X##endif
X#         tmpRow = mat[p];
X#         mat[p] = mat[pivRow];
X#         mat[pivRow] = tmpRow;
X#      }
X#
X#      for (c = p;  c <= n;  c++)
X#         mat[p][c] /= piv;
X#
X#      for (r = 0;  r < n;  r++)
X#         if (r != p)
X#            for (c = n;  c >= p;  c--)
X#               mat[r][c] -= mat[r][p] * mat[p][c];
X#   }
X#
X##ifndef NO_DETERM
X#   return (det);
X##endif
X#}
X#
X#
X#pointer Calloc(n, s)
X#   unsigned n, s;
X#{
X#   register pointer ret;
X#
X#   if ((ret = (pointer)calloc(n, s)) == NULL)
X#   {
X#      perror("Allocating matrix");
X#      exit(1);
X#   }
X#
X#   return (ret);
X#}
SHAR_EOF
fi
exit 0
#	End of shell archive
-- 
LIVE:	Ted Stefanik, (617) 965-8480 x14
USPS:	Adelie Corporation, 288 Walnut St., Newtonville, MA  02160
UUCP:	..!{wanginst!infinet | decvax!cca}!emacs!adelie!ted
ARPA:	emacs!adelie!ted-unix.ARPA



More information about the Alt.sources mailing list