You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

56 lines
1.5 KiB

  1. /*Extended Yet Another Curve Fitting Solution (The poor man's RMA)
  2. An extension of yacfs.mod adding a Weight parameter:
  3. When set to 1 the model produces best fit by least squares with all error in y and none in x (YonX);
  4. When set to zero the model produces best fit by least squares with all error in x and none in y (XonY);
  5. When set to 0.5 the model assumes equal error in x and y producing results similar to fitting by Reduced Major Axis Analysis.
  6. Nigel_Galloway@operamail.com
  7. November 5th., 2009
  8. */
  9. set Sample;
  10. param Sx {z in Sample};
  11. param Sy {z in Sample};
  12. param Weight := 0.5;
  13. var a;
  14. var b;
  15. var p;
  16. var q;
  17. XonY1 :sum{z in Sample} q*Sy[z]*Sy[z] + sum{z in Sample} p*Sy[z] = sum{z in Sample} Sy[z]*Sx[z];
  18. XonY2 :sum{z in Sample} q*Sy[z] + sum{z in Sample} p = sum{z in Sample} Sx[z];
  19. YonX1 :sum{z in Sample} a*Sx[z]*Sx[z] + sum{z in Sample} b*Sx[z] = sum{z in Sample} Sy[z]*Sx[z];
  20. YonX2 :sum{z in Sample} a*Sx[z] + sum{z in Sample} b = sum{z in Sample} Sy[z];
  21. solve;
  22. param W := Weight*a + (1-Weight)*(1/q);
  23. printf "\nbest linear fit is:\n\ty = %f %s %fx\n\n", b*Weight - (1-Weight)*(p/q), if W < 0 then "-" else "+", abs(W);
  24. data;
  25. param:
  26. Sample: Sx Sy :=
  27. 1 0 1
  28. 2 0.5 0.9
  29. 3 1 0.7
  30. 4 1.5 1.5
  31. 5 1.9 2
  32. 6 2.5 2.4
  33. 7 3 3.2
  34. 8 3.5 2
  35. 9 4 2.7
  36. 10 4.5 3.5
  37. 11 5 1
  38. 12 5.5 4
  39. 13 6 3.6
  40. 14 6.6 2.7
  41. 15 7 5.7
  42. 16 7.6 4.6
  43. 17 8.5 6
  44. 18 9 6.8
  45. 19 10 7.3
  46. ;
  47. end;