home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / qtawk / curvefit.exp < prev    next >
Encoding:
Text File  |  1990-04-30  |  11.5 KB  |  425 lines
  1. # Utility to find least squares best fit
  2. # Copyright (C) 1990 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.  
  55. # array to hold strings of equations solved, y indepedent variable
  56.     curve_strx[1] = "x = (y - a)/b";
  57.     curve_strx[2] = "x = b/(y - a)";
  58.     curve_strx[3] = "x = ((1/y) - a)/b";
  59.     curve_strx[4] = "x = a/((1/y) - b)";
  60.     curve_strx[5] = "x = ln(y/a)/b";
  61.     curve_strx[6] = "x = b/ln(y/a)";
  62.     curve_strx[7] = "x = e^((y - a)/b)";
  63.     curve_strx[8] = "x = e^(((1/y) - a)/b)";
  64.     curve_strx[9] = "x = e^(ln(y/a)/b)";
  65.  
  66. # array to hold strings of equations solved, x indepedent variable
  67.     curve_stry[1] = "y = a + bx";
  68.     curve_stry[2] = "y = a + (b/x)";
  69.     curve_stry[3] = "y = 1/(a + bx)";
  70.     curve_stry[4] = "y = x/(a + bx)";
  71.     curve_stry[5] = "y = ae^(bx)";
  72.     curve_stry[6] = "y = ae^(b/x)";
  73.     curve_stry[7] = "y = a + b*ln(x)";
  74.     curve_stry[8] = "y = 1/(a + b*ln(x))";
  75.     curve_stry[9] = "y = ax^b";
  76.  
  77. # array to hold strings of equations to be solved for x with "execute" function
  78.     curve_x_exe[1] = "(y - ca)/cb;";
  79.     curve_x_exe[2] = "cb / (y - ca);";
  80.     curve_x_exe[3] = "((1/y) - ca)/cb;";
  81.     curve_x_exe[4] = "ca/((1/y) - cb);";
  82.     curve_x_exe[5] = "log(y/ca)/cb;";
  83.     curve_x_exe[6] = "cb/log(y/ca);";
  84.     curve_x_exe[7] = "exp((y - ca)/cb);";
  85.     curve_x_exe[8] = "exp(((1/y) - ca)/cb);";
  86.     curve_x_exe[9] = "exp(log(y/ca)/cb);";
  87.  
  88. # array to hold strings of equations to be solved for y with "execute" function
  89.     curve_y_exe[1] = "ca + cb * z;";
  90.     curve_y_exe[2] = "ca + (cb / z);";
  91.     curve_y_exe[3] = "1/(ca + cb * z);";
  92.     curve_y_exe[4] = "x/(ca + cb * z);";
  93.     curve_y_exe[5] = "ca * exp(cb * z);";
  94.     curve_y_exe[6] = "ca * exp(cb / z);";
  95.     curve_y_exe[7] = "ca + cb * log(z);";
  96.     curve_y_exe[8] = "1/(ca + cb * log(z));";
  97.     curve_y_exe[9] = "ca * z ^ cb;";
  98.  
  99. # array to hold functions to be used to solve for 'a' coeffiecent when
  100. # forcing through a point. use "execute" function to solve
  101.     fcoefa[1] = "force_y - coef_b * force_x;";
  102.     fcoefa[2] = "force_y - coef_b * ifx;";
  103.     fcoefa[3] = "ify - coef_b * force_x;";
  104.     fcoefa[4] = "ify - coef_b * ifx;";
  105.     fcoefa[5] = "force_y / exp(coef_b * force_x);";
  106.     fcoefa[6] = "force_y / exp(coef_b * ifx);";
  107.     fcoefa[7] = "force_y - coef_b * lnfx;";
  108.     fcoefa[8] = "ify - coef_b * lnfx;";
  109.     fcoefa[9] = "force_y / force_x ^ coef_b;";
  110.  
  111.     stderr = "stderr";
  112.     copyright;
  113. }
  114.  
  115. INITIAL {
  116.     init_sums;
  117. }
  118.  
  119. # solve current data to specified or best fit
  120. Solve {
  121.     local cftn = TRUE;
  122.  
  123.     if ( NF > 1 && $2 ~~ /{_i}/ ) {
  124.     curve_no = $2;
  125.     cftn = specified_curve(curve_no,TRUE);
  126.     } else curve_no = best_fit;
  127.     if ( cftn ) print_curve(curve_no,a[curve_no],b[curve_no],r[curve_no]);
  128.     next;
  129. }
  130.  
  131. # new curve fit problem, re-initialize
  132. New {
  133.     init_sums;
  134.     next;
  135. }
  136.  
  137. # force fit through a given point, must re-initialize
  138. Force {
  139.     init_sums;
  140.     if ( NF > 1 && $2 ~~ input_data ) force_x = $2 + 0.0;
  141.     if ( NF > 1 && $3 ~~ input_data ) force_y = $3 + 0.0;
  142.     if ( force_x <= 0 ) {
  143.     ifx = lnfx = 0.0;
  144.     if ( force_x == 0 ) crv[2] = crv[4] = crv[6] = FALSE;
  145.       else {
  146.         crv[7] = crv[8] = crv[9] = FALSE;
  147.         ifx = 1/force_x;
  148.     }
  149.       } else {
  150.     lnfx = log(force_x);
  151.     ifx = 1/force_x;
  152.     }
  153.     if ( force_y <= 0 ) {
  154.     ify = lnfy = 0.0;
  155.     if ( force_y == 0 ) crv[3] = crv[4] = crv[8] = FALSE;
  156.       else {
  157.         crv[5] = crv[6] = crv[9] = FALSE;
  158.         ify = 1/force_y;
  159.     }
  160.       } else {
  161.     lnfy = log(force_y);
  162.     ify = 1/force_y;
  163.     }
  164.     force_pt = TRUE;
  165.     next;
  166. }
  167.  
  168. # predict x from y
  169. Pred_x {
  170.     local y = $2 + 0.0;
  171.     local fmts = "y = " ∩ $2 ∩ ", X = %.14g\n";
  172.  
  173.     if ( curve_no ) {
  174.     print "Predict: " ∩ curve_strx[curve_no];
  175.     printf(fmts,pred_x(curve_no,y));
  176.     }
  177.     next;
  178. }
  179.  
  180. # predict y from x
  181. Pred_y {
  182.     local x = $2 + 0.0;
  183.     local fmts = "x = " ∩ $2 ∩ ", Y = %.14g\n";
  184.  
  185.     if ( curve_no ) {
  186.     print "Predict: " ∩ curve_stry[curve_no];
  187.     printf(fmts,pred_y(curve_no,x));
  188.     }
  189.     next;
  190. }
  191.  
  192. # input data
  193. data_line {
  194.     accumulate_point($1 + 0.0,$2 + 0.0);
  195. # to print input data, un-comment following line
  196. #    print "x = " ∩ $1 ∩ "y = " ∩ $2;
  197.     ocurve = FALSE;
  198. }
  199.  
  200. FINAL {
  201.     if ( !curve_no ) curve_no = best_fit;
  202.     if ( !ocurve ) print_curve(curve_no,a[curve_no],b[curve_no],r[curve_no]);
  203. }
  204.  
  205. function print_curve(cn,a,b,r2) {
  206.     local sofmt = OFMT;
  207.  
  208.     OFMT = "%.14g";
  209.     print "Fit Curve Number: " ∩ cn;
  210.     print curve_stry[cn];
  211.     print "a = " ∩ a;
  212.     print "b = " ∩ b;
  213.     print "Goodness: " ∩ r2;
  214.     ocurve = TRUE;
  215.     OFMT = sofmt;
  216. }
  217.  
  218. function init_sums() {
  219.     sum_x = 0.0;
  220.     sum_y = 0.0;
  221.     sum_x2 = 0.0;
  222.     sum_y2 = 0.0;
  223.     sum_xy = 0.0;
  224.     sum_ix = 0.0;
  225.     sum_iy = 0.0;
  226.     sum_ix2 = 0.0;
  227.     sum_iy2 = 0.0;
  228.     sum_ixy = 0.0;
  229.     sum_lnx = 0.0;
  230.     sum_lny = 0.0;
  231.     sum_lnx2 = 0.0;
  232.     sum_lny2 = 0.0;
  233.     sum_lnxy = 0.0;
  234.     sum_xlny = 0.0;
  235.     sum_ixlny = 0.0;
  236.     sum_ylnx = 0.0;
  237.     sum_iylnx = 0.0;
  238.     sum_yix = 0.0;
  239.     sum_xiy = 0.0;
  240.     num_points = curve_no = 0;
  241.     ocurve = FALSE;
  242.     for ( i = 1 ; i < 10 ; i++ ) crv[i] = TRUE;
  243.  
  244.     force_pt = FALSE;
  245.     force_x = 0.0;
  246.     force_y = 0.0;
  247.     lnfx = 0.0;
  248.     lnfy = 0.0;
  249.     ifx = 0.0;
  250.     ify = 0.0;
  251. }
  252.  
  253. function accumulate_point(x,y) {
  254.     local invx, lnx;
  255.     local invy, lny;
  256.  
  257.     if ( x <= 0 ) {
  258.     invx = -ifx;
  259.     lnx = 0.0;
  260.     if ( x == 0 ) crv[2] = crv[4] = crv[6] = FALSE;
  261.       else {
  262.         crv[7] = crv[8] = crv[9] = FALSE;
  263.         invx = (1/x) - ifx;
  264.     }
  265.       } else {
  266.     invx = (1/x) - ifx;
  267.     lnx = log(x) - lnfx;
  268.     }
  269.     if ( y <= 0 ) {
  270.     invy = -ify;
  271.     lny = 0.0;
  272.     if ( y == 0 ) crv[3] = crv[4] = crv[8] = FALSE;
  273.       else {
  274.         crv[5] = crv[6] = crv[9] = FALSE;
  275.         invy = (1/y) - ify;
  276.     }
  277.       } else {
  278.     invy = (1/y) - ify;
  279.     lny = log(y) - lnfy;
  280.     }
  281.     x -= force_x;
  282.     y -= force_y;
  283.     accum_pt(x,invx,lnx,y,invy,lny);
  284.     if ( force_pt ) accum_pt(-x,-invx,-lnx,-y,-invy,-lny);
  285. }
  286.  
  287. function accum_pt(x,ix,lnx,y,iy,lny) {
  288.     local x2 = x * x, y2 = y * y;
  289.     local ix2 = ix * ix, iy2 = iy * iy;
  290.  
  291.     sum_x += x;
  292.     sum_y += y;
  293.     sum_x2 += x2;
  294.     sum_y2 += y2;
  295.     sum_xy += x * y;
  296.     sum_ix += ix;
  297.     sum_ix2 += ix2;
  298.     sum_ixlny += lny * ix;
  299.     sum_yix += y * ix;
  300.     sum_iy += iy;
  301.     sum_iy2 += iy2;
  302.     sum_iylnx += lnx * iy;
  303.     sum_xiy += x * iy;
  304.     sum_ixy += ix * iy;
  305.     sum_lnx += lnx;
  306.     sum_lny += lny;
  307.     sum_lnx2 += lnx * lnx;
  308.     sum_lny2 += lny * lny;
  309.     sum_lnxy += lnx * lny;
  310.     sum_xlny += x * lny;
  311.     sum_ylnx += y * lnx;
  312.     num_points++;
  313. }
  314.  
  315. function best_fit() {
  316.     local curve_no = FALSE, i, lr = 0;
  317.  
  318.     for ( i = 1 ; i < 10 ; i++ ) {
  319.     if ( specified_curve(i,FALSE) ) {
  320.         if ( r[i] > lr ) {
  321.         curve_no = i;
  322.         lr = r[i];
  323.         }
  324.     }
  325.     }
  326.     return curve_no;
  327. }
  328.  
  329. function specified_curve(curve_no,msg) {
  330.     local fit = FALSE, swap_ab;
  331.  
  332.     if ( crv[curve_no] ) {
  333.     switch ( curve_no ) {
  334.         case 1:
  335.         calc_a_b(sum_x,sum_y,sum_xy,sum_x2,sum_y2,num_points,1);
  336.         break;
  337.         case 2:
  338.         calc_a_b(sum_ix,sum_y,sum_yix,sum_ix2,sum_y2,num_points,2);
  339.         break;
  340.         case 3:
  341.         calc_a_b(sum_x,sum_iy,sum_xiy,sum_x2,sum_iy2,num_points,3);
  342.         break;
  343.         case 4:
  344.         calc_a_b(sum_ix,sum_iy,sum_ixy,sum_ix2,sum_iy2,num_points,4);
  345.         swap_ab = coef_a;
  346.         coef_a = coef_b;
  347.         coef_b = swap_ab;
  348.         break;
  349.         case 5:
  350.         calc_a_b(sum_x,sum_lny,sum_xlny,sum_x2,sum_lny2,num_points,5);
  351.         if ( !force_pt ) coef_a = exp(coef_a);
  352.         break;
  353.         case 6:
  354.         calc_a_b(sum_ix,sum_lny,sum_ixlny,sum_ix2,sum_lny2,num_points,6);
  355.         if ( !force_pt ) coef_a = exp(coef_a);
  356.         break;
  357.         case 7:
  358.         calc_a_b(sum_lnx,sum_y,sum_ylnx,sum_lnx2,sum_y2,num_points,7);
  359.         break;
  360.         case 8:
  361.         calc_a_b(sum_lnx,sum_iy,sum_iylnx,sum_lnx2,sum_iy2,num_points,8);
  362.         break;
  363.         case 9:
  364.         calc_a_b(sum_lnx,sum_lny,sum_lnxy,sum_lnx2,sum_lny2,num_points,9);
  365.         if ( !force_pt ) coef_a = exp(coef_a);
  366.         break;
  367.         default:
  368.         print "Improper Curve Specified";
  369.         print "File: "FILENAME;
  370.         print "Line #: "FNR;
  371.         exit 20;
  372.         break;
  373.     }
  374.     a[curve_no]  = coef_a;
  375.     b[curve_no]  = coef_b;
  376.     r[curve_no] = rr;
  377.     fit = TRUE;
  378.     } else if ( msg ) print "Data Cannot Fit Curve Selected !";
  379.     return fit;
  380. }
  381.  
  382. function calc_a_b(sum_x,sum_y,sum_xy,sum_x2,sum_y2,n,eqn) {
  383.     local sxy = sum_xy - sum_x * sum_y / n;
  384.     local ssx2 = sum_x2 - sum_x * sum_x / n;
  385.     local ssy2 = sum_y2 - sum_y * sum_y / n;
  386.  
  387.     coef_b = sxy / ssx2;
  388.     if ( force_pt ) coef_a = execute(fcoefa[eqn],TRUE);
  389.       else coef_a = (sum_y - coef_b * sum_x) / n;
  390.     rr = sqrt((sxy * sxy)/(ssx2 * ssy2));
  391. }
  392.  
  393. # function to compute x from y and curve number
  394. function pred_x(cn,y) {
  395.     local sofmt = OFMT;
  396.     local curve = curve_x_exe[cn];
  397.     local ca = a[cn], cb = b[cn];
  398.     local x;
  399.  
  400.     OFMT = "%.14g";
  401.     gsub(/y/,y,curve);
  402.     x = execute(curve,TRUE);
  403.     OFMT = sofmt;
  404.     return x;
  405. }
  406.  
  407. # function to compute y from x and curve number
  408. function pred_y(cn,x) {
  409.     local sofmt = OFMT;
  410.     local curve = curve_y_exe[cn];
  411.     local ca = a[cn], cb = b[cn];
  412.     local y;
  413.  
  414.     OFMT = "%.14g";
  415.     gsub(/z/,x,curve);
  416.     y = execute(curve,TRUE);
  417.     OFMT = sofmt;
  418.     return y;
  419. }
  420.  
  421. # function display copyright
  422. function copyright() {
  423.     fprintf(stderr,"Least Squares Best Fit.\nCopyright (C) 1990 Terry D. Boldt, All Rights Reserved.\n");
  424. }
  425.