$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         215232.624678170
        B1        1.00000000000000         236355.173469681
        B2        1.00000000000000         77934.3524331583
        B3        1.00000000000000         10147.5507550350
        B4        1.00000000000000         564.566512170752
        B5        1.00000000000000         11.2324854679312

     Residual
     Standard Deviation   236014.502379268

     R-Squared            0.957478440825662


               Certified Analysis of Variance Table

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

Regression    5      18814317208116.7      3762863441623.33   67.5524458240122
Residual     15      835542680000.000      55702845333.3333

$offtext

set i 'cases' /i1*i21/;

table data(i,*)
                 y     x
   i1         75901    0
   i2       -204794    1
   i3        204863    2
   i4       -204436    3
   i5        253665    4
   i6       -200894    5
   i7        214131    6
   i8       -185192    7
   i9        221249    8
   i10      -138370    9
   i11       315911   10
   i12       -27644   11
   i13       455253   12
   i14       197434   13
   i15       783995   14
   i16       608816   15
   i17      1370781   16
   i18      1303798   17
   i19      2205519   18
   i20      2408860   19
   i21      3444321   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";