$ontext Nonlinear Least Squares Regression example Erwin Kalvelagen, nov 2007 Reference: http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml -------------------------------------------------------------------------- Procedure: Nonlinear Least Squares Regression Description: These data are the result of a NIST study involving scanning electron microscope line with standards. Reference: Kirby, R., NIST (197?). Scanning electron microscope line width standards. Data: 1 Response (y) 1 Predictor (x) 151 Observations Average Level of Difficulty Observed Data Model: Rational Class (quadratic/quadratic) 5 Parameters (b1 to b5) y = (b1 + b2*x + b3*x**2) / (1 + b4*x + b5*x**2) + e Starting values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = 2 1.5 1.6745063063E+00 8.7989634338E-02 b2 = -0.1 -0.15 -1.3927397867E-01 4.1182041386E-03 b3 = 0.003 0.0025 2.5961181191E-03 4.1856520458E-05 b4 = -0.001 -0.0015 -1.7241811870E-03 5.8931897355E-05 b5 = 0.00001 0.00002 2.1664802578E-05 2.0129761919E-07 Residual Sum of Squares: 3.9050739624E+00 Residual Standard Deviation: 1.6354535131E-01 Degrees of Freedom: 146 Number of Observations: 151 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i151/; table data(i,*) y x i1 0.0082E0 9.65E0 i2 0.0112E0 10.74E0 i3 0.0149E0 11.81E0 i4 0.0198E0 12.88E0 i5 0.0248E0 14.06E0 i6 0.0324E0 15.28E0 i7 0.0420E0 16.63E0 i8 0.0549E0 18.19E0 i9 0.0719E0 19.88E0 i10 0.0963E0 21.84E0 i11 0.1291E0 24.00E0 i12 0.1710E0 26.25E0 i13 0.2314E0 28.86E0 i14 0.3227E0 31.85E0 i15 0.4809E0 35.79E0 i16 0.7084E0 40.18E0 i17 1.0220E0 44.74E0 i18 1.4580E0 49.53E0 i19 1.9520E0 53.94E0 i20 2.5410E0 58.29E0 i21 3.2230E0 62.63E0 i22 3.9990E0 67.03E0 i23 4.8520E0 71.25E0 i24 5.7320E0 75.22E0 i25 6.7270E0 79.33E0 i26 7.8350E0 83.56E0 i27 9.0250E0 87.75E0 i28 10.2670E0 91.93E0 i29 11.5780E0 96.10E0 i30 12.9440E0 100.28E0 i31 14.3770E0 104.46E0 i32 15.8560E0 108.66E0 i33 17.3310E0 112.71E0 i34 18.8850E0 116.88E0 i35 20.5750E0 121.33E0 i36 22.3200E0 125.79E0 i37 22.3030E0 125.79E0 i38 23.4600E0 128.74E0 i39 24.0600E0 130.27E0 i40 25.2720E0 133.33E0 i41 25.8530E0 134.79E0 i42 27.1100E0 137.93E0 i43 27.6580E0 139.33E0 i44 28.9240E0 142.46E0 i45 29.5110E0 143.90E0 i46 30.7100E0 146.91E0 i47 31.3500E0 148.51E0 i48 32.5200E0 151.41E0 i49 33.2300E0 153.17E0 i50 34.3300E0 155.97E0 i51 35.0600E0 157.76E0 i52 36.1700E0 160.56E0 i53 36.8400E0 162.30E0 i54 38.0100E0 165.21E0 i55 38.6700E0 166.90E0 i56 39.8700E0 169.92E0 i57 40.0300E0 170.32E0 i58 40.5000E0 171.54E0 i59 41.3700E0 173.79E0 i60 41.6700E0 174.57E0 i61 42.3100E0 176.25E0 i62 42.7300E0 177.34E0 i63 43.4600E0 179.19E0 i64 44.1400E0 181.02E0 i65 44.5500E0 182.08E0 i66 45.2200E0 183.88E0 i67 45.9200E0 185.75E0 i68 46.3000E0 186.80E0 i69 47.0000E0 188.63E0 i70 47.6800E0 190.45E0 i71 48.0600E0 191.48E0 i72 48.7400E0 193.35E0 i73 49.4100E0 195.22E0 i74 49.7600E0 196.23E0 i75 50.4300E0 198.05E0 i76 51.1100E0 199.97E0 i77 51.5000E0 201.06E0 i78 52.1200E0 202.83E0 i79 52.7600E0 204.69E0 i80 53.1800E0 205.86E0 i81 53.7800E0 207.58E0 i82 54.4600E0 209.50E0 i83 54.8300E0 210.65E0 i84 55.4000E0 212.33E0 i85 56.4300E0 215.43E0 i86 57.0300E0 217.16E0 i87 58.0000E0 220.21E0 i88 58.6100E0 221.98E0 i89 59.5800E0 225.06E0 i90 60.1100E0 226.79E0 i91 61.1000E0 229.92E0 i92 61.6500E0 231.69E0 i93 62.5900E0 234.77E0 i94 63.1200E0 236.60E0 i95 64.0300E0 239.63E0 i96 64.6200E0 241.50E0 i97 65.4900E0 244.48E0 i98 66.0300E0 246.40E0 i99 66.8900E0 249.35E0 i100 67.4200E0 251.32E0 i101 68.2300E0 254.22E0 i102 68.7700E0 256.24E0 i103 69.5900E0 259.11E0 i104 70.1100E0 261.18E0 i105 70.8600E0 264.02E0 i106 71.4300E0 266.13E0 i107 72.1600E0 268.94E0 i108 72.7000E0 271.09E0 i109 73.4000E0 273.87E0 i110 73.9300E0 276.08E0 i111 74.6000E0 278.83E0 i112 75.1600E0 281.08E0 i113 75.8200E0 283.81E0 i114 76.3400E0 286.11E0 i115 76.9800E0 288.81E0 i116 77.4800E0 291.08E0 i117 78.0800E0 293.75E0 i118 78.6000E0 295.99E0 i119 79.1700E0 298.64E0 i120 79.6200E0 300.84E0 i121 79.8800E0 302.02E0 i122 80.1900E0 303.48E0 i123 80.6600E0 305.65E0 i124 81.2200E0 308.27E0 i125 81.6600E0 310.41E0 i126 82.1600E0 313.01E0 i127 82.5900E0 315.12E0 i128 83.1400E0 317.71E0 i129 83.5000E0 319.79E0 i130 84.0000E0 322.36E0 i131 84.4000E0 324.42E0 i132 84.8900E0 326.98E0 i133 85.2600E0 329.01E0 i134 85.7400E0 331.56E0 i135 86.0700E0 333.56E0 i136 86.5400E0 336.10E0 i137 86.8900E0 338.08E0 i138 87.3200E0 340.60E0 i139 87.6500E0 342.57E0 i140 88.1000E0 345.08E0 i141 88.4300E0 347.02E0 i142 88.8300E0 349.52E0 i143 89.1200E0 351.44E0 i144 89.5400E0 353.93E0 i145 89.8500E0 355.83E0 i146 90.2500E0 358.32E0 i147 90.5500E0 360.20E0 i148 90.9300E0 362.67E0 i149 91.2000E0 364.53E0 i150 91.5500E0 367.00E0 i151 92.2000E0 371.30E0 ; * * extract data * parameter x(i),y(i); x(i) = data(i,'x'); y(i) = data(i,'y'); * * certified values * scalars cb1 'certified value for b1' / 1.6745063063E+00 / cb2 'certified value for b2' / -1.3927397867E-01 / cb3 'certified value for b3' / 2.5961181191E-03 / cb4 'certified value for b4' / -1.7241811870E-03 / cb5 'certified value for b5' / 2.1664802578E-05 / ce1 'certified std err for b1 ' / 8.7989634338E-02 / ce2 'certified std err for b2 ' / 4.1182041386E-03 / ce3 'certified std err for b3 ' / 4.1856520458E-05 / ce4 'certified std err for b4 ' / 5.8931897355E-05 / ce5 'certified std err for b5 ' / 2.0129761919E-07 / ; *----------------------------------------------------------------------------- * statistical model *----------------------------------------------------------------------------- variables sse 'sum of squared errors' b1 'coefficient to estimate' b2 'coefficient to estimate' b3 'coefficient to estimate' b4 'coefficient to estimate' b5 'coefficient to estimate' ; equations fit(i) 'the non-linear model' obj 'objective' ; obj.. sse =n= 0; fit(i).. y(i) =e= (b1 + b2*x(i) + b3*x(i)**2) / (1 + b4*x(i) + b5*x(i)**2); *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l = 2; b2.l = -0.1; b3.l = 0.003; b4.l = -0.001; b5.l = 0.00001; option nlp=nls; model nlfit /obj,fit/; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.l,b4.l,b5.l; abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3) +abs(b4.l-cb4)+abs(b5.l-cb5))>0.0001) "Accuracy problem"; abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3) +abs(b4.m-ce4)+abs(b5.m-ce5))>0.0001) "Accuracy problem"; *----------------------------------------------------------------------------- * second set of initial values *----------------------------------------------------------------------------- b1.l = 1.5; b2.l = -0.15; b3.l = 0.0025; b4.l = -0.0015; b5.l = 0.00002; solve nlfit minimizing sse using nlp; display sse.l,b1.l,b2.l,b3.l,b4.l,b5.l; abort$((abs(b1.l-cb1)+abs(b2.l-cb2)+abs(b3.l-cb3) +abs(b4.l-cb4)+abs(b5.l-cb5))>0.0001) "Accuracy problem"; abort$((abs(b1.m-ce1)+abs(b2.m-ce2)+abs(b3.m-ce3) +abs(b4.m-ce4)+abs(b5.m-ce5))>0.0001) "Accuracy problem";