subpack-glpk
Line | |
---|
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; |
---|
Note: See
TracBrowser
for help on using the repository browser.