home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Misc / The_SOURCE#7.LHA / Source / polyfit.lha / polyfit / fit.c next >
Encoding:
C/C++ Source or Header  |  1990-05-04  |  3.0 KB  |  159 lines
  1. #include "fit.h"
  2.  
  3. static void ParseArgs();
  4. static void Usage();
  5.  
  6. static void ReadPts();
  7. static void PolyFit();
  8. static void Results();
  9. static void Stats();
  10.  
  11.  
  12. static points Points;
  13. static int    MaxPt;
  14. static int    PtCnt;
  15. static int    Order;
  16.  
  17. static matrix Mat;
  18. static double Det;
  19.  
  20.  
  21. main(argc, argv)
  22.    int  argc;
  23.    char **argv;
  24. {
  25.    ParseArgs(argc, argv);
  26.  
  27.    ReadPts();
  28.    if (PtCnt <= 0)
  29.       Usage();
  30.    PolyFit();
  31.    Results();
  32.  
  33.    Stats();
  34. }
  35.  
  36.  
  37. static void ParseArgs(argc, argv)
  38.    int argc;
  39.    char **argv;
  40. {
  41.    if (argc != 3)
  42.       Usage(argv[0]);
  43.  
  44.    MaxPt = atoi(argv[1]);
  45.    Order = atoi(argv[2]);
  46.  
  47.    if (MaxPt <= 0 || Order <= 0)
  48.       Usage(argv[0]);
  49. }
  50.  
  51.  
  52. static void Usage(prog)
  53.    char *prog;
  54. {
  55.    fprintf(stderr,
  56.            "Usage:  %s max-number-of-data-points order-of-polynomial\n", prog);
  57.    fprintf(stderr,
  58.            "Note:   points are read from stdin (a coordinate point/line)\n");
  59.    exit(1);
  60. }
  61.  
  62.  
  63. static void ReadPts()
  64. {
  65.    char line[1024];
  66.  
  67.    Points = (points)Calloc(MaxPt, sizeof(point));
  68.  
  69.    for (PtCnt = 0;  PtCnt < MaxPt;  PtCnt++)
  70.    {
  71.       if (gets(line) == NULL)
  72.          return;
  73.       sscanf(line, "%f %f", &Points[PtCnt].X, &Points[PtCnt].Y);
  74.    }
  75. }
  76.  
  77.  
  78. static void PolyFit()
  79. {
  80.    register int    i, p, r, c;
  81.    register matrow xvec;
  82.  
  83.    Mat = MatInit(Order + 1);
  84.    xvec = (matrow)Calloc(Order * 2 + 1, sizeof(matelm));
  85.  
  86.    for (p = 0;  p < PtCnt;  p++)
  87.    {
  88.       for (xvec[0] = 1.0, i = 1;  i <= Order * 2;  i++)
  89.          xvec[i] = xvec[i - 1] * Points[p].X;
  90.       for (r = 0;  r < Order + 1;  r++)
  91.       {
  92.          Mat[r][Order + 1] += Points[p].Y * xvec[r];
  93.          for (c = 0;  c < Order + 1;  c++)
  94.             Mat[r][c] += xvec[r + c];
  95.       }
  96.    }
  97.  
  98.    Det = MatInv(Order + 1, Mat);
  99. }
  100.  
  101.  
  102. static void Results()
  103. {
  104.    register int i;
  105.  
  106.    printf("%5d = Number of data points input\n", PtCnt);
  107.    printf("%5d = Order of polynomial fit\n\n",   Order);
  108.  
  109. #ifndef NO_DETERM
  110.    printf("%e = Determinate of matrix solution\n\n", Det);
  111. #endif
  112.  
  113.    printf("y = %f", Mat[0][Order + 1]);
  114.    for (i = 1;  i < Order + 1;  i++)
  115.       printf("  +  %f x^%d", Mat[i][Order + 1], i);
  116. }
  117.  
  118.  
  119. static void Stats()
  120. {
  121.    register int i, j;
  122.    double   x, y;
  123.    double   sum, sum2, minerr, maxerr;
  124.  
  125.    sum = sum2 = maxerr = 0;
  126.    minerr = HUGE;
  127.  
  128.    for (i = 0;  i < PtCnt;  i++)
  129.    {
  130.       y = Mat[0][Order + 1];
  131.       for (x = 1, j = 1;  j < Order + 1;  j++)
  132.       {
  133.          x *= Points[i].X;
  134.          y += x * Mat[j][Order + 1];
  135.       }
  136.  
  137.       y -= Points[i].Y;
  138.       if (y < 0)
  139.          y = -y;
  140.  
  141.       sum  += y;
  142.       sum2 += y * y;
  143.  
  144.       if (y < minerr)
  145.          minerr = y;
  146.       if (y > maxerr)
  147.          maxerr = y;
  148.    }
  149.  
  150.    printf("\n\nFit statistics:\n");
  151.  
  152.    printf("%30.15f = Standard Error of Estimate\n", sqrt(sum2 / (PtCnt - 2)));
  153.    printf("%30.15f = Average Error\n", sum / PtCnt);
  154.    printf("%30.15f = Error's Standard Deviation\n", 
  155.              sqrt((PtCnt * sum2 - sum * sum) / PtCnt / (PtCnt - 1)));
  156.    printf("%30.15f = Minimum Error\n", minerr);
  157.    printf("%30.15f = Maximum Error\n", maxerr);
  158. }
  159.