$ontext

  Linear Least Squares Regression

  NIST test data

  Erwin kalvelagen, dec 2004

  Reference:
               http://www.itl.nist.gov/div898/strd/lls/lls.shtml


               Wampler, R. H. (1970).
               A Report of the Accuracy of Some Widely-Used Least
               Squares Computer Programs.
               Journal of the American Statistical Association, 65, pp. 549-565.


Model:         Polynomial Class
               6 Parameters (B0,B1,...,B5)

               y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)

               Certified Regression Statistics

                                          Standard Deviation
     Parameter          Estimate             of Estimate

        B0        1.00000000000000         21523262.4678170
        B1        1.00000000000000         23635517.3469681
        B2        1.00000000000000         7793435.24331583
        B3        1.00000000000000         1014755.07550350
        B4        1.00000000000000         56456.6512170752
        B5        1.00000000000000         1123.24854679312

     Residual
     Standard Deviation   23601450.2379268

     R-Squared            0.224668921574940E-02


               Certified Analysis of Variance Table

Source of Degrees of    Sums of               Mean
Variation  Freedom      Squares              Squares           F Statistic

Regression    5    18814317208116.7      3762863441623.33 6.7552445824012241E-03
Residual     15    0.835542680000000E+16 557028453333333.

$offtext

set i 'cases' /i1*i21/;

table data(i,*)
                   y     x
   i1        7590001     0
   i2      -20479994     1
   i3       20480063     2
   i4      -20479636     3
   i5       25231365     4
   i6      -20476094     5
   i7       20489331     6
   i8      -20460392     7
   i9       18417449     8
   i10     -20413570     9
   i11      20591111    10
   i12     -20302844    11
   i13      18651453    12
   i14     -20077766    13
   i15      21059195    14
   i16     -19666384    15
   i17      26348481    16
   i18     -18971402    17
   i19      22480719    18
   i20     -17866340    19
   i21      10958421    20
;


set j /j0*j5/;

set j1(j); j1(j)$(ord(j)>1) = yes;
parameter v(j); v(j) = ord(j)-1;

parameter x(i,j);
x(i,'j0') = 1;
x(i,j1) = power(data(i,'x'),v(j1));
display x;

variables
   b(j)  'coefficients to estimate'
   sse   'sum of squared errors'
;

equation
   fit(i)    'equation to fit'
   sumsq
;

sumsq..   sse =n= 0;
fit(i)..  data(i,'y') =e= sum(j, b(j)*x(i,j));

option lp = ls;
model leastsq /fit,sumsq/;
solve leastsq using lp minimizing sse;
option decimals=8;
display b.l;

parameter Bcert(j);
Bcert(j) = 1;

scalar err "Sum of squared errors in estimates";
err = sum(j, sqr(bcert(j)-b.l(j)));
display err;
abort$(err>0.0001) "Solution not accurate";