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