$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        0.000000000000000
        B1        1.00000000000000        0.000000000000000
        B2        1.00000000000000        0.000000000000000
        B3        1.00000000000000        0.000000000000000
        B4        1.00000000000000        0.000000000000000
        B5        1.00000000000000        0.000000000000000

     Residual
     Standard Deviation   0.000000000000000

     R-Squared            1.00000000000000


               Certified Analysis of Variance Table

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

Regression    5      18814317208116.7     3762863441623.33       Infinity
Residual     15      0.000000000000000    0.000000000000000



$offtext

set i 'cases' /i1*i21/;

table data(i,*)
                 y     x
   i1            1     0
   i2            6     1
   i3           63     2
   i4          364     3
   i5         1365     4
   i6         3906     5
   i7         9331     6
   i8        19608     7
   i9        37449     8
   i10       66430     9
   i11      111111    10
   i12      177156    11
   i13      271453    12
   i14      402234    13
   i15      579195    14
   i16      813616    15
   i17     1118481    16
   i18     1508598    17
   i19     2000719    18
   i20     2613660    19
   i21     3368421    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";