alpar@1: /*Extended Yet Another Curve Fitting Solution (The poor man's RMA)
alpar@1: 
alpar@1:   An extension of yacfs.mod adding a Weight parameter:
alpar@1:     When set to 1 the model produces best fit by least squares with all error in y and none in x (YonX);
alpar@1:     When set to zero the model produces best fit by least squares with all error in x and none in y (XonY);
alpar@1:     When set to 0.5 the model assumes equal error in x and y producing results similar to fitting by Reduced Major Axis Analysis.
alpar@1: 
alpar@1:   Nigel_Galloway@operamail.com
alpar@1:   November 5th., 2009
alpar@1: */
alpar@1: set Sample;
alpar@1: param Sx {z in Sample};
alpar@1: param Sy {z in Sample};
alpar@1: param Weight := 0.5;
alpar@1: 
alpar@1: var a;
alpar@1: var b;
alpar@1: var p;
alpar@1: var q;
alpar@1: 
alpar@1: 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];
alpar@1: XonY2 :sum{z in Sample} q*Sy[z] + sum{z in Sample} p = sum{z in Sample} Sx[z];
alpar@1: 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];
alpar@1: YonX2 :sum{z in Sample} a*Sx[z] + sum{z in Sample} b = sum{z in Sample} Sy[z];
alpar@1: 
alpar@1: solve;
alpar@1: 
alpar@1: param W := Weight*a + (1-Weight)*(1/q);
alpar@1: 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);
alpar@1: 
alpar@1: data;
alpar@1: 
alpar@1: param:
alpar@1: Sample:   Sx    Sy :=
alpar@1:   1         0    1
alpar@1:   2       0.5  0.9
alpar@1:   3         1  0.7
alpar@1:   4       1.5  1.5
alpar@1:   5       1.9    2
alpar@1:   6       2.5  2.4
alpar@1:   7         3  3.2
alpar@1:   8       3.5    2
alpar@1:   9         4  2.7
alpar@1:  10       4.5  3.5
alpar@1:  11         5    1
alpar@1:  12       5.5    4
alpar@1:  13         6  3.6
alpar@1:  14       6.6  2.7
alpar@1:  15         7  5.7
alpar@1:  16       7.6  4.6
alpar@1:  17       8.5    6
alpar@1:  18         9  6.8
alpar@1:  19        10  7.3
alpar@1: ;
alpar@1: 
alpar@1: end;