options double crt; name NLS 'NIST StRD Nonlinear Least Squares benchmarks'; load; mform(nrow=27,ncol=5) tab; ? NB1 NB2 NB3 NSE P/F mform(nrow=27,ncol=3) ctab; ? conv1 conv2 conv3 mform(nrow=27,ncol=3) nstab; ? non-singular solution table const it,0 npass,0; supres smpl nob coef; title 'Misra1a'; frml eq yM = b1*(1-exp[-b2*xM]); const nc,2 nyx,0 nob,14; read(nrow=nc,ncol=4,type=gen) vals; 500 250 2.3894212918E+02 2.7070075241E+00 0.0001 0.0005 5.5015643181E-04 7.2668688436E-06 smpl 1,nob; read yM xM; 10.07 77.6 14.73 114.9 17.94 141.1 23.93 190.8 29.61 239.9 35.18 289.0 40.02 332.8 44.82 378.4 50.76 434.8 55.05 477.3 61.01 536.8 66.40 593.1 75.47 689.1 81.78 760.0 cert; title 'Chwirut2'; frml eq y = exp(-b1*x)/(b2+b3*x); const nc,3 nyx,0 nob,54; read(nrow=nc,ncol=4,type=gen) vals; 0.1 0.15 1.6657666537E-01 3.8303286810E-02 0.01 0.008 5.1653291286E-03 6.6621605126E-04 0.02 0.010 1.2150007096E-02 1.5304234767E-03 smpl 1,nob; read(file='nls.dat') y x; ? read first dataset from file outside of proc cert; title 'Chwirut1'; frml eq y = exp(-b1*x)/(b2+b3*x); const nc,3 nyx,2 nob,214; read(nrow=nc,ncol=4,type=gen) vals; 0.1 0.15 1.9027818370E-01 2.1938557035E-02 0.01 0.008 6.1314004477E-03 3.4500025051E-04 0.02 0.010 1.0530908399E-02 7.9281847748E-04 cert; title 'Lanczos3'; frml eq yL = b1*exp(-b2*xL) + b3*exp(-b4*xL) + b5*exp(-b6*xL); const nc,6 nyx,0 nob,24; read(nrow=nc,ncol=4,type=gen) vals; 1.2 0.5 8.6816414977E-02 1.7197908859E-02 0.3 0.7 9.5498101505E-01 9.7041624475E-02 5.6 3.6 8.4400777463E-01 4.1488663282E-02 5.5 4.2 2.9515951832E+00 1.0766312506E-01 6.5 4 1.5825685901E+00 5.8371576281E-02 7.6 6.3 4.9863565084E+00 3.4436403035E-02 smpl 1,nob; trend xL 0 .05; ? 0 .05 .10 .15 ... 1.15 read yL; ? 4-5 digits of accuracy 2.5134 2.0443 1.6684 1.3664 1.1232 .9269 .7679 .6389 .5338 .4479 .3776 .3197 .2720 .2325 .1997 .1723 .1493 .1301 .1138 .1000 .0883 .0783 .0698 .0624 cert; title 'Gauss1'; frml eq y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 ) + b6*exp( -(x-b7)**2 / b8**2 ); const nc,8 nyx,1 nob,250; read(nrow=nc,ncol=4,type=gen) vals; 97.0 94.0 9.8778210871E+01 5.7527312730E-01 0.009 0.0105 1.0497276517E-02 1.1406289017E-04 100.0 99.0 1.0048990633E+02 5.8831775752E-01 65.0 63.0 6.7481111276E+01 1.0460593412E-01 20.0 25.0 2.3129773360E+01 1.7439951146E-01 70.0 71.0 7.1994503004E+01 6.2622793913E-01 178.0 180.0 1.7899805021E+02 1.2436988217E-01 16.5 20.0 1.8389389025E+01 2.0134312832E-01 smpl 1,nob; trend x; cert; title 'Gauss2'; frml eq y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 ) + b6*exp( -(x-b7)**2 / b8**2 ); const nc,8 nyx,1 nob,250; read(nrow=nc,ncol=4,type=gen) vals; 96.0 98.0 9.9018328406E+01 5.3748766879E-01 0.009 0.0105 1.0994945399E-02 1.3335306766E-04 103.0 103.0 1.0188022528E+02 5.9217315772E-01 106.0 105.0 1.0703095519E+02 1.5006798316E-01 18.0 20.0 2.3578584029E+01 2.2695595067E-01 72.0 73.0 7.2045589471E+01 6.1721965884E-01 151.0 150.0 1.5327010194E+02 1.9466674341E-01 18.0 20.0 1.9525972636E+01 2.6416549393E-01 smpl 1,nob; trend x; cert; title 'DanielWood'; frml eq y = b1*x**b2; const nc,2 nyx,0 nob,6; read(nrow=nc,ncol=4,type=gen) vals; 1 0.7 7.6886226176E-01 1.8281973860E-02 5 4 3.8604055871E+00 5.1726610913E-02 smpl 1,nob; read y x; 2.138 1.309 3.421 1.471 3.597 1.490 4.340 1.565 4.882 1.611 5.660 1.680 cert; title 'Misra1b'; frml eq yM = b1 * (1-(1+b2*xM/2)**(-2)); const nc,2 nyx,0 nob,14; ? use data from Misra1a read(nrow=nc,ncol=4,type=gen) vals; 500 300 3.3799746163E+02 3.1643950207E+00 0.0001 0.0002 3.9039091287E-04 4.2547321834E-06 cert; title 'Kirby2'; frml eq y = (b1 + b2*x + b3*x**2) / (1 + b4*x + b5*x**2); const nc,5 nyx,2 nob,151; read(nrow=nc,ncol=4,type=gen) vals; 2 1.5 1.6745063063E+00 8.7989634338E-02 -0.1 -0.15 -1.3927397867E-01 4.1182041386E-03 0.003 0.0025 2.5961181191E-03 4.1856520458E-05 -0.001 -0.0015 -1.7241811870E-03 5.8931897355E-05 0.00001 0.00002 2.1664802578E-05 2.0129761919E-07 cert; title 'Hahn1'; frml eq y = (b1+b2*x+b3*x**2+b4*x**3) / (1+b5*x+b6*x**2+b7*x**3); const nc,7 nyx,2 nob,236; read(nrow=nc,ncol=4,type=gen) vals; 10 1 1.0776351733E+00 1.7070154742E-01 -1 -0.1 -1.2269296921E-01 1.2000289189E-02 0.05 0.005 4.0863750610E-03 2.2508314937E-04 -0.00001 -0.000001 -1.4262662514E-06 2.7578037666E-07 -0.05 -0.005 -5.7609940901E-03 2.4712888219E-04 0.001 0.0001 2.4053735503E-04 1.0449373768E-05 -0.000001 -0.0000001 -1.2314450199E-07 1.3027335327E-08 cert; title 'Nelson'; frml eq log[y] - (b1 - b2*x1 * exp[-b3*x2]); const nc,3 nyx,3 nob,128; read(nrow=nc,ncol=4,type=gen) vals; 2 2.5 2.5906836021E+00 1.9149996413E-02 0.0001 0.000000005 5.6177717026E-09 6.1124096540E-09 -0.01 -0.05 -5.7701013174E-02 3.9572366543E-03 cert; title 'MGH17'; ?frml eq y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5]; ? Note: order of b4 and b3 swapped, so that the order of coefs ? in the output @COEF and @SES is b1 b2 b3 b4 b5 frml eq y = b1 + b2*exp[-x*b3] + b4*exp[-x*b5]; const nc,5 nyx,0 nob,33; read(nrow=nc,ncol=4,type=gen) vals; ? Note: order of b4 and b3 also swapped below 50 0.5 3.7541005211E-01 2.0723153551E-03 150 1.5 1.9358469127E+00 2.2031669222E-01 1 0.01 1.2867534640E-02 4.4861358114E-04 -100 -1 -1.4646871366E+00 2.2175707739E-01 2 0.02 2.2122699662E-02 8.9471996575E-04 smpl 1,nob; trend x 0 10; ? 0 10 20 ... 320 read y; .844 .908 .932 .936 .925 .908 .881 .850 .818 .784 .751 .718 .685 .658 .628 .603 .580 .558 .538 .522 .506 .490 .478 .467 .457 .448 .438 .431 .424 .420 .414 .411 .406 nosupres coef; ? turn on output to diagnose bad results cert; supres coef; title 'Lanczos1'; frml eq yL = b1*exp(-b2*xL) + b3*exp(-b4*xL) + b5*exp(-b6*xL); const nc,6 nyx,0 nob,24; read(nrow=nc,ncol=4,type=gen) vals; 1.2 0.5 9.5100000027E-02 5.3347304234E-11 0.3 0.7 1.0000000001E+00 2.7473038179E-10 5.6 3.6 8.6070000013E-01 1.3576062225E-10 5.5 4.2 3.0000000002E+00 3.3308253069E-10 6.5 4 1.5575999998E+00 1.8815731448E-10 7.6 6.3 5.0000000001E+00 1.1057500538E-10 smpl 1,nob; read yL; ? 14 digits of accuracy 2.5134 2.044333373291 1.668404436564 1.366418021208 1.123232487372 .9268897180037 .7679338563728 .6388775523106 .5337835317402 .4479363617347 .3775847884350 .3197393199326 .2720130773746 .2324965529032 .1996589546065 .1722704126914 .1493405660168 .1300700206922 .1138119324644 .1000415587559 .08833209084540 .07833544019350 .06976693743449 .06239312536719 cert; title 'Lanczos2'; frml eq yL = b1*exp(-b2*xL) + b3*exp(-b4*xL) + b5*exp(-b6*xL); const nc,6 nyx,0 nob,24; read(nrow=nc,ncol=4,type=gen) vals; 1.2 0.5 9.6251029939E-02 6.6770575477E-04 0.3 0.7 1.0057332849E+00 3.3989646176E-03 5.6 3.6 8.6424689056E-01 1.7185846685E-03 5.5 4.2 3.0078283915E+00 4.1707005856E-03 6.5 4 1.5529016879E+00 2.3744381417E-03 7.6 6.3 5.0028798100E+00 1.3958787284E-03 smpl 1,nob; read yL; ? 6 digits of accuracy 2.51340 2.04433 1.66840 1.36642 1.12323 .926890 .767934 .638878 .533784 .447936 .377585 .319739 .272013 .232497 .199659 .172270 .149341 .130070 .113812 .100042 .0883321 .0783354 .0697669 .0623931 cert; title 'Gauss3'; frml eq y = b1*exp( -b2*x ) + b3*exp( -(x-b4)**2 / b5**2 ) + b6*exp( -(x-b7)**2 / b8**2 ) ; const nc,8 nyx,1 nob,250; read(nrow=nc,ncol=4,type=gen) vals; 94.9 96.0 9.8940368970E+01 5.3005192833E-01 0.009 0.0096 1.0945879335E-02 1.2554058911E-04 90.1 80.0 1.0069553078E+02 8.1256587317E-01 113.0 110.0 1.1163619459E+02 3.5317859757E-01 20.0 25.0 2.3300500029E+01 3.6584783023E-01 73.8 74.0 7.3705031418E+01 1.2091239082E+00 140.0 139.0 1.4776164251E+02 4.0488183351E-01 20.0 25.0 1.9668221230E+01 3.7806634336E-01 smpl 1,nob; trend x; cert; title 'Misra1c'; frml eq yM = b1 * (1-(1+2*b2*xM)**(-.5)); const nc,2 nyx,0 nob,14; read(nrow=nc,ncol=4,type=gen) vals; 500 600 6.3642725809E+02 4.6638326572E+00 0.0001 0.0002 2.0813627256E-04 1.7728423155E-06 cert; title 'Misra1d'; frml eq yM = b1*b2*xM*((1+b2*xM)**(-1)); const nc,2 nyx,0 nob,14; read(nrow=nc,ncol=4,type=gen) vals; 500 450 4.3736970754E+02 3.6489174345E+00 0.0001 0.0003 3.0227324449E-04 2.9334354479E-06 cert; title 'Roszman1'; set pi = 3.141592653589793238462643383279; frml eq y = b1 - b2*x - atan[b3/(x-b4)]/pi; const nc,4 nyx,2 nob,25; read(nrow=nc,ncol=4,type=gen) vals; 0.1 0.2 2.0196866396E-01 1.9172666023E-02 -0.00001 -0.000005 -6.1953516256E-06 3.2058931691E-06 1000 1200 1.2044556708E+03 7.4050983057E+01 -100 -150 -1.8134269537E+02 4.9573513849E+01 cert; title 'ENSO'; ?frml eq y = b1 + b2*cos( 2*pi*x/12 ) + b3*sin( 2*pi*x/12 ) ? + b5*cos( 2*pi*x/b4 ) + b6*sin( 2*pi*x/b4 ) ? + b8*cos( 2*pi*x/b7 ) + b9*sin( 2*pi*x/b7 ) ; ? Note: as in MGH17, we have to reorder the parameter names ? slightly, so that they will appear in @COEF and @SES in the ? expected ordering. Swap b4/b5 and b7/b8. frml eq y = b1 + b2*cos( 2*pi*x/12 ) + b3*sin( 2*pi*x/12 ) + b4*cos( 2*pi*x/b5 ) + b6*sin( 2*pi*x/b5 ) + b7*cos( 2*pi*x/b8 ) + b9*sin( 2*pi*x/b8 ) ; const nc,9 nyx,1 nob,168; read(nrow=nc,ncol=4,type=gen) vals; ? swap b4/b5 and b7/b8 11.0 10.0 1.0510749193E+01 1.7488832467E-01 3.0 3.0 3.0762128085E+00 2.4310052139E-01 0.5 0.5 5.3280138227E-01 2.4354686618E-01 -0.7 -1.5 -1.6231428586E+00 2.8078369611E-01 40.0 44.0 4.4311088700E+01 9.4408025976E-01 -1.3 0.5 5.2554493756E-01 4.8073701119E-01 -0.3 -0.1 2.1232288488E-01 5.1460022911E-01 25.0 26.0 2.6887614440E+01 4.1612939130E-01 1.4 1.5 1.4966870418E+00 2.5434468893E-01 smpl 1,nob; trend x; cert; title 'MGH09'; frml eq y = b1*(x**2+x*b2) / (x**2+x*b3+b4); const nc,4 nyx,0 nob,11; read(nrow=nc,ncol=4,type=gen) vals; 25 0.25 1.9280693458E-01 1.1435312227E-02 39 0.39 1.9128232873E-01 1.9633220911E-01 41.5 0.415 1.2305650693E-01 8.0842031232E-02 39 0.39 1.3606233068E-01 9.0025542308E-02 smpl 1,nob; read y x; .1957 4 .1947 2 .1735 1 .1600 .500 .0844 .250 .0627 .167 .0456 .125 .0342 .100 .0323 .0833 .0235 .0714 .0246 .0625 cert; title 'Thurber'; frml eq y = (b1 + b2*x + b3*x**2 + b4*x**3) / (1 + b5*x + b6*x**2 + b7*x**3); const nc,7 nyx,2 nob,37; read(nrow=nc,ncol=4,type=gen) vals; 1000 1300 1.2881396800E+03 4.6647963344E+00 1000 1500 1.4910792535E+03 3.9571156086E+01 400 500 5.8323836877E+02 2.8698696102E+01 40 75 7.5416644291E+01 5.5675370270E+00 0.7 1 9.6629502864E-01 3.1333340687E-02 0.3 0.4 3.9797285797E-01 1.4984928198E-02 0.03 0.05 4.9727297349E-02 6.5842344623E-03 cert; title 'BoxBOD'; frml eq y = b1*(1-exp[-b2*x]); const nc,2 nyx,0 nob,6; read(nrow=nc,ncol=4,type=gen) vals; 1 100 2.1380940889E+02 1.2354515176E+01 1 0.75 5.4723748542E-01 1.0455993237E-01 smpl 1,nob; read y x; 109 1 149 2 149 3 191 5 213 7 224 10 cert; title 'Ratkowsky2'; frml eq y = b1 / (1+exp[b2-b3*x]); const nc,3 nyx,0 nob,9; read(nrow=nc,ncol=4,type=gen) vals; 100 75 7.2462237576E+01 1.7340283401E+00 1 2.5 2.6180768402E+00 8.8295217536E-02 0.1 0.07 6.7359200066E-02 3.4465663377E-03 smpl 1,nob; read y x; 8.93 9 10.80 14 18.59 21 22.33 28 39.35 42 56.11 57 61.73 63 64.62 70 67.08 79 cert; title 'MGH10'; frml eq y = b1 * exp[b2/(x+b3)]; const nc,3 nyx,0 nob,16; read(nrow=nc,ncol=4,type=gen) vals; 2 0.02 5.6096364710E-03 1.5687892471E-04 400000 4000 6.1813463463E+03 2.3309021107E+01 25000 250 3.4522363462E+02 7.8486103508E-01 smpl 1,nob; trend x 50 5; ? 50 55 ... 125 read y; 34780 28610 23650 19630 16370 13720 11540 9744 8261 7030 6005 5147 4427 3820 3307 2872 nosupres coef; ? turn on output to diagnose bad results cert; supres coef; title 'Eckerle4'; frml eq y = (b1/b2) * exp[-0.5*((x-b3)/b2)**2]; const nc,3 nyx,2 nob,35; read(nrow=nc,ncol=4,type=gen) vals; 1 1.5 1.5543827178E+00 1.5408051163E-02 10 5 4.0888321754E+00 4.6803020753E-02 500 450 4.5154121844E+02 4.6800518816E-02 cert; title 'Ratkowsky3'; frml eq y = b1 / ((1+exp[b2-b3*x])**(1/b4)); const nc,4 nyx,0 nob,15; read(nrow=nc,ncol=4,type=gen) vals; 100 700 6.9964151270E+02 1.6302297817E+01 10 5 5.2771253025E+00 2.0828735829E+00 1 0.75 7.5962938329E-01 1.9566123451E-01 1 1.3 1.2792483859E+00 6.8761936385E-01 smpl 1,nob; trend x; read y; 16.08 33.83 65.80 97.20 191.55 326.20 386.87 520.53 590.03 651.92 724.93 699.56 689.96 637.56 717.41 cert; title 'Bennett5'; frml eq y = b1 * (b2+x)**(-1/b3); const nc,3 nyx,2 nob,154; read(nrow=nc,ncol=4,type=gen) vals; -2000 -1500 -2.5235058043E+03 2.9715175411E+02 50 45 4.6736564644E+01 1.2448871856E+00 0.8 0.85 9.3218483193E-01 2.0272299378E-02 cert; ?print tab,ctab; write(format="(/' NLS results'/ 9X,' ----NDigits Beta---- NDigits'/ ' Model ',' Start1 Start2 Start3 SE(B) Pass=1/Fail=0')"); smpl 1,27; unmake tab c1-c5; mmake tab cnames c1-c5; unmake ctab c1-c3; unmake nstab c4-c6; mmake ctab cnames c1-c6; write(format="(1X,A8,4F7.1,F5.0)") tab; write(format="(' NLS summary: ',F3.0,'/27 passed.')") npass; write(format="(/' Convergence flags , Non-singular flags')"); write(format="(1X,A8,3F3.0,',',3F3.0)") ctab; ?====================================================== Proc Cert; ? Procedure for running NIST StRD NLS benchmarks -- ? Run with 3 sets of starting values, record number of ? correct digits in table. ? Inputs: (global variables) ? nc = number of coefficients ? vals = nc x 4 matrix of starting and certified values ? nob = number of observations ? nyx = number of y, x series to read from nls.dat (can be zero) ? eq = nonlinear equation to estimate, containing parameters ? b1 b2 ... b(nc), as given in the NIST files ? tab = table of output results ? tabc = table of convergence flags ? it = pointer to row of tab and tabc for this estimation ? (maintained here) ? npass = sum of 5th column of tab (number of models that pass) ? by Clint Cummins 8/97 list(prefix=b,last=nc) bs; param bs; ? unmake matrix of start and certified values into columns ? Certified coef values are in column 3, which are also used ? as the 3rd set of starting values (to check iteration behavior ? when started at the optimum). smpl 1,nc; unmake vals s1s s2s s3s scs; mmake s1 s1s; mmake s2 s2s; mmake s3 s3s; mmake sc scs; ? read data if not already defined smpl 1,nob; if (nyx=1); then; read(file='nls.dat') y; if (nyx=2); then; read(file='nls.dat') y x; if (nyx=3); then; read(file='nls.dat') y x1 x2; if (nyx>3); then; title 'ERROR: nyx > 3'; ? prepare for estimation set it = it+1; ? pointer to row of tables print it; dot(index=is) 1 2 3; unmake s. bs; lsq(tol=1e-11,maxit=100,maxsqz=25,silent) eq; set ctab(it,is) = @ifconv; ? ? Alternative definition of convergence -- if ? norm of gradient is < 1.e-10 ? if .not.@ifconv; then; do; mat crit = @grad'@vcov*@grad; set ctab(it,is) = crit < 1.e-10; enddo; set nstab(it,is) = @ncoef=@ncid; ? equal for non-singular solution ? min number of digits for coefs mat dc = abs[(@coef-s3)/(s3 + (s3=0))]; ? denom is replaced by 1 if s3=0 mat dc = (dc >= 1e-15)%dc + (dc < 1e-15)*1e-15; ? avoid log of zero mat nb. = min( -log10(dc) ); set tab(it,is) = nb.; ? min number of digits for SEs mat dc = abs[(@ses-sc)/sc]; mat dc = (dc >= 1e-15)%dc + (dc < 1e-15)*1e-15; mat ns. = min( -log10(dc) ); copy @ssr ssr.; enddot; if (ssr2 < ssr1); then; do; copy nb2 nbopt; copy ns2 nsopt; enddo; else; do; copy nb1 nbopt; copy ns1 nsopt; enddo; set tab(it,4) = nsopt; ? N digits SE for min SSR solution set ipass = (nbopt >= 4) & (nsopt >= 4); ? Pass/Fail criteria set tab(it,5) = ipass; set npass = npass+ipass; endproc; ?====================================================== end; noprint; smpl 1,27; read(format='(A8)') names; Misra1a Chwirut2 Chwirut2 Lanczos3 Gauss1 Gauss2 DanielWo Misra1b Kirby2 Hahn1 Nelson MGH17 Lanczos1 Lanczos2 Gauss3 Misra1c Misra1d Roszman1 ENSO MGH09 Thurber BoxBOD Ratkows2 MGH10 Eckerle4 Ratkows3 Bennett5 mmake cnames names;