octave_packages/optim-1.2.0/test_wpolyfit.m

   1 ## Author: Paul Kienzle <pkienzle@gmail.com>
   2 ## This program is granted to the public domain.
   3
   4 ## Tests for wpolyfit.
   5 ##
   6 ## Test cases are taken from the NIST Statistical Reference Datasets
   7 ##    http://www.itl.nist.gov/div898/strd/
   8
   9 1;
  10
  11 function do_test(n,x,y,p,dp,varargin)
  12   [myp,s] = wpolyfit(x,y,n,varargin{:});
  13   %if length(varargin)==0, [myp,s] = polyfit(x,y,n); else return; end
  14   mydp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
  15   if length(varargin)>0, mydp = [mydp;0]; end %origin
  16   %[svdp,j,svddp] = svdfit(x,y,n);
  17   disp('parameter  certified value  rel. error');
  18   [myp(:), p, abs((myp(:)-p)./p)] %, svdp, abs((svdp-p)./p) ]
  19   disp('p-error    certified value  rel. error');
  20   [mydp(:), dp, abs((mydp(:) - dp)./dp)] %, svdp, abs((svddp - dp)./dp)]
  21   input('Press <Enter> to proceed to the next test');
  22 endfunction
  23
  24 ##          x         y          dy
  25 data = [0.0013852  0.2144023  0.0020470
  26         0.0018469  0.2516856  0.0022868
  27         0.0023087  0.3070443  0.0026362
  28         0.0027704  0.3603186  0.0029670
  29         0.0032322  0.4260864  0.0033705
  30         0.0036939  0.4799956  0.0036983 ];
  31 x=data(:,1); y=data(:,2); dy=data(:,3);
  32 wpolyfit(x,y,dy,1);
  33 disp('computing parameter uncertainty from monte carlo simulation...');
  34 fflush(stdout);
  35 n=100; p=zeros(2,n);
  36 for i=1:n, p(:,i)=(polyfit(x,y+randn(size(y)).*dy,1)).'; end
  37 printf("%15s %15s\n", "Coefficient", "Error");
  38 printf("%15g %15g\n", [mean(p'); std(p')]);
  39 input('Press <Enter> to see some sample regression lines: ');
  40 t = [x(1), x(length(x))];
  41 [p,s] = wpolyfit(x,y,dy,1); dp=sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
  42 hold off;
  43 for i=1:15, plot(t,polyval(p(:)+randn(size(dp)).*dp,t),'-g;;'); hold on; end
  44 errorbar(x,y,dy,"~b;;");
  45 [yf,dyf]=polyconf(p,x,s,0.05,'ci');
  46 plot(x,yf-dyf,"-r;;",x,yf+dyf,'-r;95% confidence interval;')
  47 hold off;
  48 input('Press <Enter> to continue with the tests: ');
  49
  50
  51 ##Procedure:     Linear Least Squares Regression
  52 ##Reference:     Filippelli, A., NIST.
  53 ##Model:         Polynomial Class
  54 ##               11 Parameters (B0,B1,...,B10)
  55 ##
  56 ##               y = B0 + B1*x + B2*(x**2) + ... + B9*(x**9) + B10*(x**10) + e
  57
  58 ##Data:
  59 ##            y          x
  60 data = [    0.8116   -6.860120914
  61             0.9072   -4.324130045
  62             0.9052   -4.358625055
  63             0.9039   -4.358426747
  64             0.8053   -6.955852379
  65             0.8377   -6.661145254
  66             0.8667   -6.355462942
  67             0.8809   -6.118102026
  68             0.7975   -7.115148017
  69             0.8162   -6.815308569
  70             0.8515   -6.519993057
  71             0.8766   -6.204119983
  72             0.8885   -5.853871964
  73             0.8859   -6.109523091
  74             0.8959   -5.79832982
  75             0.8913   -5.482672118
  76             0.8959   -5.171791386
  77             0.8971   -4.851705903
  78             0.9021   -4.517126416
  79             0.909    -4.143573228
  80             0.9139   -3.709075441
  81             0.9199   -3.499489089
  82             0.8692   -6.300769497
  83             0.8872   -5.953504836
  84             0.89     -5.642065153
  85             0.891    -5.031376979
  86             0.8977   -4.680685696
  87             0.9035   -4.329846955
  88             0.9078   -3.928486195
  89             0.7675   -8.56735134
  90             0.7705   -8.363211311
  91             0.7713   -8.107682739
  92             0.7736   -7.823908741
  93             0.7775   -7.522878745
  94             0.7841   -7.218819279
  95             0.7971   -6.920818754
  96             0.8329   -6.628932138
  97             0.8641   -6.323946875
  98             0.8804   -5.991399828
  99             0.7668   -8.781464495
 100             0.7633   -8.663140179
 101             0.7678   -8.473531488
 102             0.7697   -8.247337057
 103             0.77     -7.971428747
 104             0.7749   -7.676129393
 105             0.7796   -7.352812702
 106             0.7897   -7.072065318
 107             0.8131   -6.774174009
 108             0.8498   -6.478861916
 109             0.8741   -6.159517513
 110             0.8061   -6.835647144
 111             0.846    -6.53165267
 112             0.8751   -6.224098421
 113             0.8856   -5.910094889
 114             0.8919   -5.598599459
 115             0.8934   -5.290645224
 116             0.894    -4.974284616
 117             0.8957   -4.64454848
 118             0.9047   -4.290560426
 119             0.9129   -3.885055584
 120             0.9209   -3.408378962
 121             0.9219   -3.13200249
 122             0.7739   -8.726767166
 123             0.7681   -8.66695597
 124             0.7665   -8.511026475
 125             0.7703   -8.165388579
 126             0.7702   -7.886056648
 127             0.7761   -7.588043762
 128             0.7809   -7.283412422
 129             0.7961   -6.995678626
 130             0.8253   -6.691862621
 131             0.8602   -6.392544977
 132             0.8809   -6.067374056
 133             0.8301   -6.684029655
 134             0.8664   -6.378719832
 135             0.8834   -6.065855188
 136             0.8898   -5.752272167
 137             0.8964   -5.132414673
 138             0.8963   -4.811352704
 139             0.9074   -4.098269308
 140             0.9119   -3.66174277
 141             0.9228   -3.2644011];
 142
 143 ##Certified values:
 144 ##                      p                       dP
 145 target = [      -1467.48961422980         298.084530995537
 146                 -2772.17959193342         559.779865474950
 147                 -2316.37108160893         466.477572127796
 148                 -1127.97394098372         227.204274477751
 149                 -354.478233703349         71.6478660875927
 150                 -75.1242017393757         15.2897178747400
 151                 -10.8753180355343         2.23691159816033
 152                 -1.06221498588947         0.221624321934227
 153                 -0.670191154593408E-01    0.142363763154724E-01
 154                 -0.246781078275479E-02    0.535617408889821E-03
 155                 -0.402962525080404E-04    0.896632837373868E-05];
 156 if 1
 157   disp("Filippelli, A.,  NIST.");
 158   do_test(10, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 159 endif
 160
 161 ##Procedure:     Linear Least Squares Regression
 162 ##
 163 ##Reference:     Pontius, P., NIST.
 164 ##               Load Cell Calibration.
 165 ##
 166 ##Model:         Quadratic Class
 167 ##               3 Parameters (B0,B1,B2)
 168 ##               y = B0 + B1*x + B2*(x**2)
 169
 170
 171 ##Data:       y             x
 172 data = [ \
 173          .11019        150000
 174          .21956        300000
 175          .32949        450000
 176          .43899        600000
 177          .54803        750000
 178          .65694        900000
 179          .76562       1050000
 180          .87487       1200000
 181          .98292       1350000
 182         1.09146       1500000
 183         1.20001       1650000
 184         1.30822       1800000
 185         1.41599       1950000
 186         1.52399       2100000
 187         1.63194       2250000
 188         1.73947       2400000
 189         1.84646       2550000
 190         1.95392       2700000
 191         2.06128       2850000
 192         2.16844       3000000
 193          .11052        150000
 194          .22018        300000
 195          .32939        450000
 196          .43886        600000
 197          .54798        750000
 198          .65739        900000
 199          .76596       1050000
 200          .87474       1200000
 201          .98300       1350000
 202         1.09150       1500000
 203         1.20004       1650000
 204         1.30818       1800000
 205         1.41613       1950000
 206         1.52408       2100000
 207         1.63159       2250000
 208         1.73965       2400000
 209         1.84696       2550000
 210         1.95445       2700000
 211         2.06177       2850000
 212         2.16829       3000000 ];
 213
 214 ##               Certified Regression Statistics
 215 ##
 216 ##                                          Standard Deviation
 217 ##                     Estimate             of Estimate
 218 target = [ \
 219               0.673565789473684E-03    0.107938612033077E-03
 220               0.732059160401003E-06    0.157817399981659E-09
 221              -0.316081871345029E-14    0.486652849992036E-16 ];
 222
 223 if 1
 224   disp("Pontius, P., NIST");
 225   do_test(2, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 226 endif
 227
 228
 229
 230 #Procedure:     Linear Least Squares Regression
 231 #Reference:     Eberhardt, K., NIST.
 232 #Model:         Linear Class
 233 #               1 Parameter (B1)
 234 #
 235 #               y = B1*x + e
 236
 237 #Data:     y     x
 238 data =[\
 239          130    60
 240          131    61
 241          132    62
 242          133    63
 243          134    64
 244          135    65
 245          136    66
 246          137    67
 247          138    68
 248          139    69
 249          140    70 ];
 250
 251
 252 #               Certified Regression Statistics
 253 #
 254 #                                 Standard Deviation
 255 #               Estimate             of Estimate
 256 target = [ \
 257           0                    0
 258           2.07438016528926     0.165289256198347E-01 ];
 259
 260
 261 if 1
 262   disp("Eberhardt, K., NIST");
 263   do_test(1, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)),'origin');
 264 endif
 265
 266
 267 #Reference:     Wampler, R. H. (1970).
 268 #               A Report of the Accuracy of Some Widely-Used Least
 269 #               Squares Computer Programs.
 270 #               Journal of the American Statistical Association, 65, 549-565.
 271 #
 272 #Model:         Polynomial Class
 273 #               6 Parameters (B0,B1,...,B5)
 274 #
 275 #               y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
 276 #
 277 #               Certified Regression Statistics
 278 #
 279 #                                          Standard Deviation
 280 #     Parameter        Estimate               of Estimate
 281 target = [\
 282                 1.00000000000000        0.000000000000000
 283                 1.00000000000000        0.000000000000000
 284                 1.00000000000000        0.000000000000000
 285                 1.00000000000000        0.000000000000000
 286                 1.00000000000000        0.000000000000000
 287                 1.00000000000000        0.000000000000000 ];
 288
 289 #Data:            y     x
 290 data = [\
 291                  1     0
 292                  6     1
 293                 63     2
 294                364     3
 295               1365     4
 296               3906     5
 297               9331     6
 298              19608     7
 299              37449     8
 300              66430     9
 301             111111    10
 302             177156    11
 303             271453    12
 304             402234    13
 305             579195    14
 306             813616    15
 307            1118481    16
 308            1508598    17
 309            2000719    18
 310            2613660    19
 311            3368421    20 ];
 312
 313 if 1
 314   disp("Wampler1");
 315   do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 316 endif
 317
 318 ##Reference:     Wampler, R. H. (1970).
 319 ##               A Report of the Accuracy of Some Widely-Used Least
 320 ##               Squares Computer Programs.
 321 ##               Journal of the American Statistical Association, 65, 549-565.
 322 ##Model:         Polynomial Class
 323 ##               6 Parameters (B0,B1,...,B5)
 324 ##
 325 ##               y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
 326 ##
 327 ##               Certified Regression Statistics
 328 ##                                       Standard Deviation
 329 ## Parameter         Estimate               of Estimate
 330 target = [ \
 331                 1.00000000000000         0.000000000000000
 332                 0.100000000000000        0.000000000000000
 333                 0.100000000000000E-01    0.000000000000000
 334                 0.100000000000000E-02    0.000000000000000
 335                 0.100000000000000E-03    0.000000000000000
 336                 0.100000000000000E-04    0.000000000000000 ];
 337
 338
 339 #Data:          y       x
 340 data = [ \
 341             1.00000    0
 342             1.11111    1
 343             1.24992    2
 344             1.42753    3
 345             1.65984    4
 346             1.96875    5
 347             2.38336    6
 348             2.94117    7
 349             3.68928    8
 350             4.68559    9
 351             6.00000   10
 352             7.71561   11
 353             9.92992   12
 354            12.75603   13
 355            16.32384   14
 356            20.78125   15
 357            26.29536   16
 358            33.05367   17
 359            41.26528   18
 360            51.16209   19
 361            63.00000   20 ];
 362
 363 if 1
 364   disp("Wampler2");
 365   do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 366 endif
 367
 368
 369
 370
 371 ##Reference:   Wampler, R. H. (1970).
 372 ##             A Report of the Accuracy of Some Widely-Used Least
 373 ##             Squares Computer Programs.
 374 ##             Journal of the American Statistical Association, 65, 549-565.
 375 ##
 376 ##Model:       Polynomial Class
 377 ##             6 Parameters (B0,B1,...,B5)
 378 ##
 379 ##             y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
 380 ##
 381 ##             Certified Regression Statistics
 382 ##
 383 ##                                        Standard Deviation
 384 ##   Parameter          Estimate             of Estimate
 385 target = [\
 386                   1.00000000000000         2152.32624678170
 387                   1.00000000000000         2363.55173469681
 388                   1.00000000000000         779.343524331583
 389                   1.00000000000000         101.475507550350
 390                   1.00000000000000         5.64566512170752
 391                   1.00000000000000         0.112324854679312    ];
 392
 393 #Data:           y      x
 394 data = [ \
 395               760.     0
 396             -2042.     1
 397              2111.     2
 398             -1684.     3
 399              3888.     4
 400              1858.     5
 401             11379.     6
 402             17560.     7
 403             39287.     8
 404             64382.     9
 405            113159.    10
 406            175108.    11
 407            273291.    12
 408            400186.    13
 409            581243.    14
 410            811568.    15
 411           1121004.    16
 412           1506550.    17
 413           2002767.    18
 414           2611612.    19
 415           3369180.    20 ];
 416 if 1
 417   disp("Wampler3");
 418   do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 419 endif
 420
 421 ##Model:         Polynomial Class
 422 ##               6 Parameters (B0,B1,...,B5)
 423 ##
 424 ##               y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
 425 ##
 426 ##              Certified Regression Statistics
 427 ##
 428 ##                                          Standard Deviation
 429 ##     Parameter          Estimate             of Estimate
 430 target = [\
 431                   1.00000000000000         215232.624678170
 432                   1.00000000000000         236355.173469681
 433                   1.00000000000000         77934.3524331583
 434                   1.00000000000000         10147.5507550350
 435                   1.00000000000000         564.566512170752
 436                   1.00000000000000         11.2324854679312 ];
 437
 438 #Data:            y     x
 439 data = [\
 440               75901    0
 441             -204794    1
 442              204863    2
 443             -204436    3
 444              253665    4
 445             -200894    5
 446              214131    6
 447             -185192    7
 448              221249    8
 449             -138370    9
 450              315911   10
 451              -27644   11
 452              455253   12
 453              197434   13
 454              783995   14
 455              608816   15
 456             1370781   16
 457             1303798   17
 458             2205519   18
 459             2408860   19
 460             3444321   20 ];
 461
 462 if 1
 463   disp("Wampler4");
 464   do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 465 endif
 466
 467
 468
 469 ##Model:         Polynomial Class
 470 ##               6 Parameters (B0,B1,...,B5)
 471 ##
 472 ##               y = B0 + B1*x + B2*(x**2) + B3*(x**3)+ B4*(x**4) + B5*(x**5)
 473 ##
 474 ##               Certified Regression Statistics
 475 ##
 476 ##                                          Standard Deviation
 477 ##     Parameter          Estimate             of Estimate
 478 target = [\
 479                   1.00000000000000         21523262.4678170
 480                   1.00000000000000         23635517.3469681
 481                   1.00000000000000         7793435.24331583
 482                   1.00000000000000         1014755.07550350
 483                   1.00000000000000         56456.6512170752
 484                   1.00000000000000         1123.24854679312 ];
 485
 486 ##Data:            y     x
 487 data = [ \
 488              7590001     0
 489            -20479994     1
 490             20480063     2
 491            -20479636     3
 492             25231365     4
 493            -20476094     5
 494             20489331     6
 495            -20460392     7
 496             18417449     8
 497            -20413570     9
 498             20591111    10
 499            -20302844    11
 500             18651453    12
 501            -20077766    13
 502             21059195    14
 503            -19666384    15
 504             26348481    16
 505            -18971402    17
 506             22480719    18
 507            -17866340    19
 508             10958421    20 ];
 509 if 1
 510   disp("Wampler5");
 511   do_test(5, data(:,2),data(:,1),flipud(target(:,1)),flipud(target(:,2)));
 512 endif