$ontext

  Linear Least Squares Regression

  NIST test data

  Erwin kalvelagen, dec 2004

  Reference:
     Filippelli, A., NIST.


  Model:       Polynomial Class
               11 Parameters (B0,B1,...,B10)

               y = B0 + B1*x + B2*(x**2) + ... + B9*(x**9) + B10*(x**10) + e




               Certified Regression Statistics

                                            Standard Deviation
     Parameter         Estimate                of Estimate

        B0        -1467.48961422980         298.084530995537
        B1        -2772.17959193342         559.779865474950
        B2        -2316.37108160893         466.477572127796
        B3        -1127.97394098372         227.204274477751
        B4        -354.478233703349         71.6478660875927
        B5        -75.1242017393757         15.2897178747400
        B6        -10.8753180355343         2.23691159816033
        B7        -1.06221498588947         0.221624321934227
        B8        -0.670191154593408E-01    0.142363763154724E-01
        B9        -0.246781078275479E-02    0.535617408889821E-03
        B10       -0.402962525080404E-04    0.896632837373868E-05

     Residual
     Standard Deviation   0.334801051324544E-02

     R-Squared            0.996727416185620


               Certified Analysis of Variance Table

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

Regression   10     0.242391619837339     0.242391619837339E-01 2162.43954511489
Residual     71     0.795851382172941E-03 0.112091743968020E-04


$offtext

set i 'cases' /i1*i82/;

table data(i,*)
              y           x
   i1       0.8116   -6.860120914
   i2       0.9072   -4.324130045
   i3       0.9052   -4.358625055
   i4       0.9039   -4.358426747
   i5       0.8053   -6.955852379
   i6       0.8377   -6.661145254
   i7       0.8667   -6.355462942
   i8       0.8809   -6.118102026
   i9       0.7975   -7.115148017
   i10      0.8162   -6.815308569
   i11      0.8515   -6.519993057
   i12      0.8766   -6.204119983
   i13      0.8885   -5.853871964
   i14      0.8859   -6.109523091
   i15      0.8959   -5.79832982
   i16      0.8913   -5.482672118
   i17      0.8959   -5.171791386
   i18      0.8971   -4.851705903
   i19      0.9021   -4.517126416
   i20      0.909    -4.143573228
   i21      0.9139   -3.709075441
   i22      0.9199   -3.499489089
   i23      0.8692   -6.300769497
   i24      0.8872   -5.953504836
   i25      0.89     -5.642065153
   i26      0.891    -5.031376979
   i27      0.8977   -4.680685696
   i28      0.9035   -4.329846955
   i29      0.9078   -3.928486195
   i30      0.7675   -8.56735134
   i31      0.7705   -8.363211311
   i32      0.7713   -8.107682739
   i33      0.7736   -7.823908741
   i34      0.7775   -7.522878745
   i35      0.7841   -7.218819279
   i36      0.7971   -6.920818754
   i37      0.8329   -6.628932138
   i38      0.8641   -6.323946875
   i39      0.8804   -5.991399828
   i40      0.7668   -8.781464495
   i41      0.7633   -8.663140179
   i42      0.7678   -8.473531488
   i43      0.7697   -8.247337057
   i44      0.77     -7.971428747
   i45      0.7749   -7.676129393
   i46      0.7796   -7.352812702
   i47      0.7897   -7.072065318
   i48      0.8131   -6.774174009
   i49      0.8498   -6.478861916
   i50      0.8741   -6.159517513
   i51      0.8061   -6.835647144
   i52      0.846    -6.53165267
   i53      0.8751   -6.224098421
   i54      0.8856   -5.910094889
   i55      0.8919   -5.598599459
   i56      0.8934   -5.290645224
   i57      0.894    -4.974284616
   i58      0.8957   -4.64454848
   i59      0.9047   -4.290560426
   i60      0.9129   -3.885055584
   i61      0.9209   -3.408378962
   i62      0.9219   -3.13200249
   i63      0.7739   -8.726767166
   i64      0.7681   -8.66695597
   i65      0.7665   -8.511026475
   i66      0.7703   -8.165388579
   i67      0.7702   -7.886056648
   i68      0.7761   -7.588043762
   i69      0.7809   -7.283412422
   i70      0.7961   -6.995678626
   i71      0.8253   -6.691862621
   i72      0.8602   -6.392544977
   i73      0.8809   -6.067374056
   i74      0.8301   -6.684029655
   i75      0.8664   -6.378719832
   i76      0.8834   -6.065855188
   i77      0.8898   -5.752272167
   i78      0.8964   -5.132414673
   i79      0.8963   -4.811352704
   i80      0.9074   -4.098269308
   i81      0.9119   -3.66174277
   i82      0.9228   -3.2644011
;

set j /j0*j10/;

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

parameter x(i,j);
x(i,j) = power(data(i,'x'),v(j));
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;

*
* compare against known solution
*

parameter Bcert /
        j0        -1467.48961422980
        j1        -2772.17959193342
        j2        -2316.37108160893
        j3        -1127.97394098372
        j4        -354.478233703349
        j5        -75.1242017393757
        j6        -10.8753180355343
        j7        -1.06221498588947
        j8        -0.670191154593408E-01
        j9        -0.246781078275479E-02
        j10       -0.402962525080404E-04
/;

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";



*
* plot results
*
* first we need to make sure x comes before y.
* we had before declared 'y' before 'x' so we introduce
* 'x0' and 'y0' where we make sure 'x0' comes first.
* if you load plotdata in the gdxviewer you will see that
* indeed 'x0','y0' are ordered correctly. If this step
* was not performed, we would have seen an inverted
* graph.
*

parameter plotdata(i,*);
plotdata(i,'x0') = EPS+data(i,'x');
plotdata(i,'y0') = EPS+data(i,'y');

set k/point1*point200/;
scalar minx, maxx, stepx;
minx = smin(i,data(i,'x'));
maxx = smax(i,data(i,'x'));
stepx = (maxx-minx)/(card(k)-1);
parameter xfit(k), yfit(k);
xfit(k) = minx+stepx*(ord(k)-1);
yfit(k) = sum(j, b.l(j)*power(xfit(k),v(j)));
parameter fitdata(k,*);
fitdata(k,'x0') =  EPS+xfit(k);
fitdata(k,'y0') =  EPS+yfit(k);


execute_unload 'chartdata.gdx',plotdata,fitdata;

$onecho > filip_gch.gch
[CHART]
VERID=GAMSIDE Chart(s) V1
GDXFILE=chartdata.gdx
TITLE=plotdata

[SERIES1]
SYMBOL=plotdata
TYPE=scatter2d

[SERIES2]
SYMBOL=fitdata
TYPE=function
$offecho


execute '=idecmds FileOpen filip_gch.gch'