ls02.gms : Test LS (Least Squares) utility - higher difficulty

Description

Run the LS solver on the Filip test problem (Level of Difficulty:
Higher) from NIST and verify that the results are those certified to
be correct.  The problem and results were both taken from the NIST Web site.

Contributor: Erwin Kalvelagen and Steve Dirkse, July 2008.


Small Model of Type : GAMS


Category : GAMS Test library


Main file : ls02.gms

$TITLE Test LS (Least Squares) utility - higher difficulty (LS02,SEQ=395)

$ontext

Run the LS solver on the Filip test problem (Level of Difficulty:
Higher) from NIST and verify that the results are those certified to
be correct.  The problem and results were both taken from the NIST Web site.

Contributor: Erwin Kalvelagen and Steve Dirkse, July 2008.

$offtext


$ontext
Model: y = sum{p=0..10, beta(p) * power(x,p)} + error

Certified Regression Statistics
                                          Standard Deviation
 Parameter        Estimate                   of Estimate
   beta(0)      -1467.48961422980         298.084530995537
   beta(1)      -2772.17959193342         559.779865474950
   beta(2)      -2316.37108160893         466.477572127796
   beta(3)      -1127.97394098372         227.204274477751
   beta(4)      -354.478233703349         71.6478660875927
   beta(5)      -75.1242017393757         15.2897178747400
   beta(6)      -10.8753180355343         2.23691159816033
   beta(7)      -1.06221498588947         0.221624321934227
   beta(8)      -0.670191154593408E-01    0.142363763154724E-01
   beta(9)      -0.246781078275479E-02    0.535617408889821E-03
   beta(10)     -0.402962525080404E-04    0.896632837373868E-05


Residual
Standard Deviation      0.334801051324544E-02
R-Squared               0.996727416185620

Certified Analysis of Variance Table

Source of
Variation       Degrees     Sums of Squares       Mean Squares          F Statistic
                of Freedom 
Regression      10          0.242391619837339     0.242391619837339E-01 2162.43954511489
Residual        71          0.795851382172941E-03 0.112091743968020E-04

$offtext


set i 'cases' /i01*i82/;

table data(i,*)
            y          x
i01       0.8116   -6.860120914
i02       0.9072   -4.324130045
i03       0.9052   -4.358625055
i04       0.9039   -4.358426747
i05       0.8053   -6.955852379
i06       0.8377   -6.661145254
i07       0.8667   -6.355462942
i08       0.8809   -6.118102026
i09       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 p  'coefficients' /
  b0   'beta(0)',
  b1   'beta(1)'
  b2   'beta(2)'
  b3   'beta(3)'
  b4   'beta(4)'
  b5   'beta(5)'
  b6   'beta(6)'
  b7   'beta(7)'
  b8   'beta(8)'
  b9   'beta(9)'
  b10  'beta(10)'
/;

variables
  x(p)
  sse 'sum of squared errors'
  ;

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

sumsq..   sse =n= 0;
fit(i)..  data(i,'y') =e= sum {p, x(p) * power(data(i,'x'),ord(p)-1)};

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

* The values with _ appended are the certified values:
* we test that we reproduce these

Scalars
  d
  tol             / 2e-4 /
  sigma_  'Standard error'     / 0.334801051324544E-02 /
  sigma   'Standard error'
  r2_ 'R-Squared' / 0.996727416185620 /
  r2  'R-Squared'
  df_  'Degrees of freedom' / 71 /
  df   'Degrees of freedom'
  rss_ 'Residual sum of squares' / 0.795851382172941E-03 /
  rss  'Residual sum of squares'
  resvar_  'Residual variance' / 0.112091743968020E-04 /
  resvar   'Residual variance'
  ;

parameters
  estimate_(p) 'Estimated coefficients' /
   'b0'      -1467.48961422980
   'b1'      -2772.17959193342
   'b2'      -2316.37108160893
   'b3'      -1127.97394098372
   'b4'      -354.478233703349
   'b5'      -75.1242017393757
   'b6'      -10.8753180355343
   'b7'      -1.06221498588947
   'b8'      -0.670191154593408E-01
   'b9'      -0.246781078275479E-02
   'b10'     -0.402962525080404E-04
  /
  estimate(p) 'Estimated coefficients'
  se_(p)  'Standard errors' /
   'b0'     298.084530995537
   'b1'     559.779865474950
   'b2'     466.477572127796
   'b3'     227.204274477751
   'b4'     71.6478660875927
   'b5'     15.2897178747400
   'b6'     2.23691159816033
   'b7'     0.221624321934227
   'b8'     0.142363763154724E-01
   'b9'     0.535617408889821E-03
   'b10'    0.896632837373868E-05
  /
  se(p)  'Standard errors'

  tval_(p)  't values' /
   'b0' -4.92306535592814 
   'b1' -4.95226746866985 
   'b2' -4.9656644637948
   'b3' -4.96458068038224 
   'b4' -4.94750585675733 
   'b5' -4.91338055354818 
   'b6' -4.86175589332647 
   'b7' -4.79286292997965 
   'b8' -4.70759656125059 
   'b9' -4.60741336969886 
   'b10' -4.49417546489266
  /
  tval(p)  't values'
  pval_(p)  'p values' /
   'b0' 5.34685204911511E-6 
   'b1' 4.78348738397472E-6 
   'b2' 4.54487872048048E-6 
   'b3' 4.56374709756346E-6 
   'b4' 4.871224064118E-6 
   'b5' 5.54761663806858E-6 
   'b6' 6.74864941174746E-6 
   'b7' 8.75364821695257E-6 
   'b8' 1.20509912449052E-5 
   'b9' 1.74863253050717E-5 
   'b10' 2.65146118181292E-5
  /
  pval(p)  'p values'
  ;


d = smax{p, abs(x.l(p)-estimate_(p))};
abort$[d > tol] 'bad solution x.l', x.l, estimate_, d;

execute_load 'ls', estimate, se, sigma, r2, df, rss, resvar, tval, pval;
abort$[execerror > 0] 'Could not load statistics from GDX';

d = abs(sigma_ - sigma);
abort$[d > tol] 'bad standard error', sigma_, sigma, d;

d = abs(r2_ - r2);
abort$[d > tol] 'bad R-Squared', r2_, r2, d;

d = abs(df_ - df);
abort$[d > tol] 'bad degrees of freedom', df_, df, d;

d = abs(rss_ - rss);
abort$[d > tol] 'bad residual sum of squares', rss_, rss, d;

d = abs(resvar_ - resvar);
abort$[d > tol] 'bad residual variance', resvar_, resvar, d;

d = smax{p, abs(estimate_(p) - estimate(p))};
abort$[d > tol] 'bad estimate', estimate_, estimate, d;

d = smax{p, abs(se_(p) - se(p))};
abort$[d > tol] 'bad standard error', se_, se, d;

d = smax{p, abs(tval_(p) - tval(p))};
abort$[d > tol] 'bad t values', tval_, tval, d;

d = smax{p, abs(pval_(p) - pval(p))};
abort$[d > tol] 'bad p values', pval_, pval, d;