subpack-glpk
Rev | Line | |
---|
[9] | 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.