home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 11 Util / 11-Util.zip / qtawkos2.zip / QTAUTL.ZIP / curvefit.exp < prev    next >
Text File  |  1994-11-21  |  13KB  |  432 lines
  1. # Utility to find least squares best fit
  2. # Copyright (C) 1990 - 1994 Terry D. Boldt, All Rights Reserved
  3. #
  4. # Least Square curve fit to best curve. May also specify point through which
  5. # desired curve must pass.
  6. # In general linearize and fit:              Y = A + BX
  7. # If forcing through a specified point, Fit: Y = BX
  8. # Fit the following:
  9. # Linear:
  10. #   1) y = a + bx       Y = y,   A = a, B = b, X = x
  11. #   2) y = a + (b/x)    Y = y,   A = a, B = b, X = 1/x                 x ╪ 0
  12. #   3) y = 1/(a + bx)   Y = 1/y, A = a, B = b, X = x            y ╪ 0
  13. #   4) y = x/(a + bx)   Y = 1/y, A = b, B = a, X = 1/x          y ╪ 0, x ╪ 0
  14. # Exponential:
  15. #   5) y = ae^(bx)      Y = ln(y), A = ln(a), B = b, X = x      y > 0
  16. #   6) y = ae^(b/x)     Y = ln(y), A = ln(a), B = b, X = 1/x    y > 0, x ╪ 0
  17. # Logarithmic:
  18. #   7) y = a + b*ln(x)      Y = y,   A = a, B = b, X = ln(x)           x > 0
  19. #   8) y = 1/(a + b*ln(x))  Y = 1/y, A = a, B = b, X = ln(x)    y ╪ 0, x > 0
  20. # Power:
  21. #   9) y = ax^b         Y = ln(y), A = ln(a), B = b, X = ln(x)  y > 0, x > 0
  22. #
  23. # Input File Format: (Data and Commands)
  24. # x_data y_data # comment
  25. # solve         # will find best fit to data input
  26. # solve n       # 1 <= n <= 9, will fit data to curve n
  27. # new           # will re-initialize
  28. # force         # force curve to fit through origin, (0,0)
  29. # force x y     # force curve to fit through point (x,y)
  30. # pred_x y      # predict x given current curve and y
  31. # pred_y x      # predict y given current curve and x
  32. #
  33. # Reference on forcing to a specified point:
  34. #  Edward B. Hands, "Forcing Regression Through a Given Point Using Any
  35. #  Familiar Computational Routine", U.S. Army, Corps of Engineers, Coastal
  36. #  Engineering Research Center, Technical Paper No. 83-1, March 1983
  37. #  Several errors in the equations are present and corrected in this
  38. #  implementation. The technique translates the origin of the co-ordinate
  39. #  system to the point to be fit, and fits the equation: y' = b*x'
  40. #  then translates back to the original system, recovering 'a' in the
  41. #  second translation
  42. #
  43. BEGIN {
  44. # command and data regular expressions
  45.     input_data = /({_i}|{_f}|{_r})/;
  46.     data_comment = /({_w}+#.*)?$/;
  47.     data = /{input_data}{_w}+{input_data}{data_comment}/;
  48.     data_line = /^{_w}*{data}/;
  49.     Solve = /^{_w}*[Ss]olve({_w}|$)/;
  50.     New = /^{_w}*[Nn]ew({_w}|$)/;
  51.     Force = /^{_w}*[Ff]orce({_w}+{data})?/;
  52.     Pred_x = /^{_w}*[Pp]red_x{_w}+{input_data}{data_comment}/;
  53.     Pred_y = /^{_w}*[Pp]red_y{_w}+{input_data}{data_comment}/;
  54.     Comment = /^{_w}*#/;
  55.  
  56. # array to hold strings of equations solved, y independent variable
  57.     curve_strx[1] = "x = (y - a)/b";
  58.     curve_strx[2] = "x = b/(y - a)";
  59.     curve_strx[3] = "x = ((1/y) - a)/b";
  60.     curve_strx[4] = "x = a/((1/y) - b)";
  61.     curve_strx[5] = "x = ln(y/a)/b";
  62.     curve_strx[6] = "x = b/ln(y/a)";
  63.     curve_strx[7] = "x = e^((y - a)/b)";
  64.     curve_strx[8] = "x = e^(((1/y) - a)/b)";
  65.     curve_strx[9] = "x = e^(ln(y/a)/b)";
  66.  
  67. # array to hold strings of equations solved, x independent variable
  68.     curve_stry[1] = "y = a + bx";
  69.     curve_stry[2] = "y = a + (b/x)";
  70.     curve_stry[3] = "y = 1/(a + bx)";
  71.     curve_stry[4] = "y = x/(a + bx)";
  72.     curve_stry[5] = "y = ae^(bx)";
  73.     curve_stry[6] = "y = ae^(b/x)";
  74.     curve_stry[7] = "y = a + b*ln(x)";
  75.     curve_stry[8] = "y = 1/(a + b*ln(x))";
  76.     curve_stry[9] = "y = ax^b";
  77.  
  78. # array to hold strings of equations to be solved for x with "execute" function
  79.     curve_x_exe[1] = "(y - ca)/cb;";
  80.     curve_x_exe[2] = "cb / (y - ca);";
  81.     curve_x_exe[3] = "((1/y) - ca)/cb;";
  82.     curve_x_exe[4] = "ca/((1/y) - cb);";
  83.     curve_x_exe[5] = "log(y/ca)/cb;";
  84.     curve_x_exe[6] = "cb/log(y/ca);";
  85.     curve_x_exe[7] = "exp((y - ca)/cb);";
  86.     curve_x_exe[8] = "exp(((1/y) - ca)/cb);";
  87.     curve_x_exe[9] = "exp(log(y/ca)/cb);";
  88.  
  89. # array to hold strings of equations to be solved for y with "execute" function
  90.     curve_y_exe[1] = "ca + cb * z;";
  91.     curve_y_exe[2] = "ca + (cb / z);";
  92.     curve_y_exe[3] = "1/(ca + cb * z);";
  93.     curve_y_exe[4] = "x/(ca + cb * z);";
  94.     curve_y_exe[5] = "ca * exp(cb * z);";
  95.     curve_y_exe[6] = "ca * exp(cb / z);";
  96.     curve_y_exe[7] = "ca + cb * log(z);";
  97.     curve_y_exe[8] = "1/(ca + cb * log(z));";
  98.     curve_y_exe[9] = "ca * z ^ cb;";
  99.  
  100. # array to hold functions to be used to solve for 'a' coeffiecent when
  101. # forcing through a point. use "execute" function to solve
  102.     fcoefa[1] = "force_y - coef_b * force_x;";
  103.     fcoefa[2] = "force_y - coef_b * ifx;";
  104.     fcoefa[3] = "ify - coef_b * force_x;";
  105.     fcoefa[4] = "ify - coef_b * ifx;";
  106.     fcoefa[5] = "force_y / exp(coef_b * force_x);";
  107.     fcoefa[6] = "force_y / exp(coef_b * ifx);";
  108.     fcoefa[7] = "force_y - coef_b * lnfx;";
  109.     fcoefa[8] = "ify - coef_b * lnfx;";
  110.     fcoefa[9] = "force_y / force_x ^ coef_b;";
  111.  
  112.     stderr = "stderr";
  113.     copyright;
  114. }
  115.  
  116. INITIAL {
  117.     init_sums;
  118. }
  119.  
  120. # comment line - pass through
  121. GROUP Comment {
  122.     print;
  123. #    next;
  124. }
  125.  
  126. # solve current data to specified or best fit
  127. GROUP Solve {
  128.     local cftn = TRUE;
  129.  
  130.     if ( NF > 1 && $2 ~~ /{_i}/ ) {
  131.         curve_no = $2;
  132.         cftn = specified_curve(curve_no,TRUE);
  133.     } else curve_no = best_fit;
  134.     if ( cftn ) print_curve(curve_no,a[curve_no],b[curve_no],r[curve_no]);
  135. #    next;
  136. }
  137.  
  138. # new curve fit problem, re-initialize
  139. GROUP New {
  140.     init_sums;
  141. #    next;
  142. }
  143.  
  144. # force fit through a given point, must re-initialize
  145. GROUP Force {
  146.     init_sums;
  147.     if ( NF > 1 && $2 ~~ input_data ) force_x = $2 + 0.0;
  148.     if ( NF > 1 && $3 ~~ input_data ) force_y = $3 + 0.0;
  149.     if ( force_x <= 0 ) {
  150.         ifx = lnfx = 0.0;
  151.         if ( force_x == 0 ) crv[2] = crv[4] = crv[6] = FALSE;
  152.           else {
  153.             crv[7] = crv[8] = crv[9] = FALSE;
  154.             ifx = 1/force_x;
  155.         }
  156.       } else {
  157.         lnfx = log(force_x);
  158.         ifx = 1/force_x;
  159.     }
  160.     if ( force_y <= 0 ) {
  161.         ify = lnfy = 0.0;
  162.         if ( force_y == 0 ) crv[3] = crv[4] = crv[8] = FALSE;
  163.           else {
  164.             crv[5] = crv[6] = crv[9] = FALSE;
  165.             ify = 1/force_y;
  166.         }
  167.       } else {
  168.         lnfy = log(force_y);
  169.         ify = 1/force_y;
  170.     }
  171.     force_pt = TRUE;
  172. #    next;
  173. }
  174.  
  175. # predict x from y
  176. GROUP Pred_x {
  177.     local y = $2 + 0.0;
  178.     local fmts = "y = " ∩ $2 ∩ ", X = %.14g\n";
  179.  
  180.     if ( curve_no ) {
  181.         print "Predict: " ∩ curve_strx[curve_no];
  182.         printf(fmts,pred_x(curve_no,y));
  183.     }
  184. #    next;
  185. }
  186.  
  187. # predict y from x
  188. GROUP Pred_y {
  189.     local x = $2 + 0.0;
  190.     local fmts = "x = " ∩ $2 ∩ ", Y = %.14g\n";
  191.  
  192.     if ( curve_no ) {
  193.         print "Predict: " ∩ curve_stry[curve_no];
  194.         printf(fmts,pred_y(curve_no,x));
  195.     }
  196. #    next;
  197. }
  198.  
  199. # input data
  200. GROUP data_line {
  201.     accumulate_point($1 + 0.0,$2 + 0.0);
  202. # to print input data, un-comment following line
  203. #    print "x = " ∩ $1 ∩ "y = " ∩ $2;
  204.     ocurve = FALSE;
  205. }
  206.  
  207. FINAL {
  208.     if ( !curve_no ) curve_no = best_fit;
  209.     if ( !ocurve ) print_curve(curve_no,a[curve_no],b[curve_no],r[curve_no]);
  210. }
  211.  
  212. function print_curve(cn,a,b,r2) {
  213.     local sofmt = OFMT;
  214.  
  215.     OFMT = "%.14g";
  216.     print "Fit Curve Number: " ∩ cn;
  217.     print curve_stry[cn];
  218.     print "a = " ∩ a;
  219.     print "b = " ∩ b;
  220.     print "Goodness: " ∩ r2;
  221.     ocurve = TRUE;
  222.     OFMT = sofmt;
  223. }
  224.  
  225. function init_sums() {
  226.     sum_x = 0.0;
  227.     sum_y = 0.0;
  228.     sum_x2 = 0.0;
  229.     sum_y2 = 0.0;
  230.     sum_xy = 0.0;
  231.     sum_ix = 0.0;
  232.     sum_iy = 0.0;
  233.     sum_ix2 = 0.0;
  234.     sum_iy2 = 0.0;
  235.     sum_ixy = 0.0;
  236.     sum_lnx = 0.0;
  237.     sum_lny = 0.0;
  238.     sum_lnx2 = 0.0;
  239.     sum_lny2 = 0.0;
  240.     sum_lnxy = 0.0;
  241.     sum_xlny = 0.0;
  242.     sum_ixlny = 0.0;
  243.     sum_ylnx = 0.0;
  244.     sum_iylnx = 0.0;
  245.     sum_yix = 0.0;
  246.     sum_xiy = 0.0;
  247.     num_points = curve_no = 0;
  248.     ocurve = FALSE;
  249.     for ( i = 1 ; i < 10 ; i++ ) crv[i] = TRUE;
  250.  
  251.     force_pt = FALSE;
  252.     force_x = 0.0;
  253.     force_y = 0.0;
  254.     lnfx = 0.0;
  255.     lnfy = 0.0;
  256.     ifx = 0.0;
  257.     ify = 0.0;
  258. }
  259.  
  260. function accumulate_point(x,y) {
  261.     local invx, lnx;
  262.     local invy, lny;
  263.  
  264.     if ( x <= 0 ) {
  265.         invx = -ifx;
  266.         lnx = 0.0;
  267.         if ( x == 0 ) crv[2] = crv[4] = crv[6] = FALSE;
  268.           else {
  269.             crv[7] = crv[8] = crv[9] = FALSE;
  270.             invx = (1/x) - ifx;
  271.         }
  272.       } else {
  273.         invx = (1/x) - ifx;
  274.         lnx = log(x) - lnfx;
  275.     }
  276.     if ( y <= 0 ) {
  277.         invy = -ify;
  278.         lny = 0.0;
  279.         if ( y == 0 ) crv[3] = crv[4] = crv[8] = FALSE;
  280.           else {
  281.             crv[5] = crv[6] = crv[9] = FALSE;
  282.             invy = (1/y) - ify;
  283.         }
  284.       } else {
  285.         invy = (1/y) - ify;
  286.         lny = log(y) - lnfy;
  287.     }
  288.     x -= force_x;
  289.     y -= force_y;
  290.     accum_pt(x,invx,lnx,y,invy,lny);
  291.     if ( force_pt ) accum_pt(-x,-invx,-lnx,-y,-invy,-lny);
  292. }
  293.  
  294. function accum_pt(x,ix,lnx,y,iy,lny) {
  295.     local x2 = x * x, y2 = y * y;
  296.     local ix2 = ix * ix, iy2 = iy * iy;
  297.  
  298.     sum_x += x;
  299.     sum_y += y;
  300.     sum_x2 += x2;
  301.     sum_y2 += y2;
  302.     sum_xy += x * y;
  303.     sum_ix += ix;
  304.     sum_ix2 += ix2;
  305.     sum_ixlny += lny * ix;
  306.     sum_yix += y * ix;
  307.     sum_iy += iy;
  308.     sum_iy2 += iy2;
  309.     sum_iylnx += lnx * iy;
  310.     sum_xiy += x * iy;
  311.     sum_ixy += ix * iy;
  312.     sum_lnx += lnx;
  313.     sum_lny += lny;
  314.     sum_lnx2 += lnx * lnx;
  315.     sum_lny2 += lny * lny;
  316.     sum_lnxy += lnx * lny;
  317.     sum_xlny += x * lny;
  318.     sum_ylnx += y * lnx;
  319.     num_points++;
  320. }
  321.  
  322. function best_fit() {
  323.     local curve_no = FALSE, i, lr = 0;
  324.  
  325.     for ( i = 1 ; i < 10 ; i++ ) {
  326.         if ( specified_curve(i,FALSE) ) {
  327.             if ( r[i] > lr ) {
  328.                 curve_no = i;
  329.                 lr = r[i];
  330.             }
  331.         }
  332.     }
  333.     return curve_no;
  334. }
  335.  
  336. function specified_curve(curve_no,msg) {
  337.     local fit = FALSE, swap_ab;
  338.  
  339.     if ( crv[curve_no] ) {
  340.         switch ( curve_no ) {
  341.             case 1:
  342.                 calc_a_b(sum_x,sum_y,sum_xy,sum_x2,sum_y2,num_points,1);
  343.                 break;
  344.             case 2:
  345.                 calc_a_b(sum_ix,sum_y,sum_yix,sum_ix2,sum_y2,num_points,2);
  346.                 break;
  347.             case 3:
  348.                 calc_a_b(sum_x,sum_iy,sum_xiy,sum_x2,sum_iy2,num_points,3);
  349.                 break;
  350.             case 4:
  351.                 calc_a_b(sum_ix,sum_iy,sum_ixy,sum_ix2,sum_iy2,num_points,4);
  352.                 swap_ab = coef_a;
  353.                 coef_a = coef_b;
  354.                 coef_b = swap_ab;
  355.                 break;
  356.             case 5:
  357.                 calc_a_b(sum_x,sum_lny,sum_xlny,sum_x2,sum_lny2,num_points,5);
  358.                 if ( !force_pt ) coef_a = exp(coef_a);
  359.                 break;
  360.             case 6:
  361.                 calc_a_b(sum_ix,sum_lny,sum_ixlny,sum_ix2,sum_lny2,num_points,6);
  362.                 if ( !force_pt ) coef_a = exp(coef_a);
  363.                 break;
  364.             case 7:
  365.                 calc_a_b(sum_lnx,sum_y,sum_ylnx,sum_lnx2,sum_y2,num_points,7);
  366.                 break;
  367.             case 8:
  368.                 calc_a_b(sum_lnx,sum_iy,sum_iylnx,sum_lnx2,sum_iy2,num_points,8);
  369.                 break;
  370.             case 9:
  371.                 calc_a_b(sum_lnx,sum_lny,sum_lnxy,sum_lnx2,sum_lny2,num_points,9);
  372.                 if ( !force_pt ) coef_a = exp(coef_a);
  373.                 break;
  374.             default:
  375.                 print "Improper Curve Specified";
  376.                 print "File: "FILENAME;
  377.                 print "Line #: "FNR;
  378.                 exit 20;
  379.                 break;
  380.         }
  381.         a[curve_no]  = coef_a;
  382.         b[curve_no]  = coef_b;
  383.         r[curve_no] = rr;
  384.         fit = TRUE;
  385.     } else if ( msg ) print "Data Cannot Fit Curve Selected !";
  386.     return fit;
  387. }
  388.  
  389. function calc_a_b(sum_x,sum_y,sum_xy,sum_x2,sum_y2,n,eqn) {
  390.     local sxy = sum_xy - sum_x * sum_y / n;
  391.     local ssx2 = sum_x2 - sum_x * sum_x / n;
  392.     local ssy2 = sum_y2 - sum_y * sum_y / n;
  393.  
  394.     coef_b = sxy / ssx2;
  395.     if ( force_pt ) coef_a = execute(fcoefa[eqn],TRUE);
  396.       else coef_a = (sum_y - coef_b * sum_x) / n;
  397.     rr = sqrt((sxy * sxy)/(ssx2 * ssy2));
  398. }
  399.  
  400. # function to compute x from y and curve number
  401. function pred_x(cn,y) {
  402.     local sofmt = OFMT;
  403.     local curve = curve_x_exe[cn];
  404.     local ca = a[cn], cb = b[cn];
  405.     local x;
  406.  
  407.     OFMT = "%.14g";
  408.     gsub(/y/,y,curve);
  409.     x = execute(curve,TRUE);
  410.     OFMT = sofmt;
  411.     return x;
  412. }
  413.  
  414. # function to compute y from x and curve number
  415. function pred_y(cn,x) {
  416.     local sofmt = OFMT;
  417.     local curve = curve_y_exe[cn];
  418.     local ca = a[cn], cb = b[cn];
  419.     local y;
  420.  
  421.     OFMT = "%.14g";
  422.     gsub(/z/,x,curve);
  423.     y = execute(curve,TRUE);
  424.     OFMT = sofmt;
  425.     return y;
  426. }
  427.  
  428. # function display copyright
  429. function copyright() {
  430.     fprintf(stderr,"Least Squares Best Fit.\nCopyright (C) 1990, 1991 Terry D. Boldt, All Rights Reserved.\n");
  431. }
  432.