$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: This problem was found to be difficult for some very good algorithms. See More, J. J., Garbow, B. S., and Hillstrom, K. E. (1981). Testing unconstrained optimization software. ACM Transactions on Mathematical Software. 7(1): pp. 17-41. Reference: Osborne, M. R. (1972). Some aspects of nonlinear least squares calculations. In Numerical Methods for Nonlinear Optimization, Lootsma (Ed). New York, NY: Academic Press, pp. 171-189. Data: 1 Response (y) 1 Predictor (x) 33 Observations Average Level of Difficulty Generated Data Model: Exponential Class 5 Parameters (b1 to b5) y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5] + e Starting values Certified Values Start 1 Start 2 Parameter Standard Deviation b1 = 50 0.5 3.7541005211E-01 2.0723153551E-03 b2 = 150 1.5 1.9358469127E+00 2.2031669222E-01 b3 = -100 -1 -1.4646871366E+00 2.2175707739E-01 b4 = 1 0.01 1.2867534640E-02 4.4861358114E-04 b5 = 2 0.02 2.2122699662E-02 8.9471996575E-04 Residual Sum of Squares: 5.4648946975E-05 Residual Standard Deviation: 1.3970497866E-03 Degrees of Freedom: 28 Number of Observations: 33 $offtext *----------------------------------------------------------------------------- * data *----------------------------------------------------------------------------- set i /i1*i33/; table data(i,*) y x i1 8.440000E-01 0.000000E+00 i2 9.080000E-01 1.000000E+01 i3 9.320000E-01 2.000000E+01 i4 9.360000E-01 3.000000E+01 i5 9.250000E-01 4.000000E+01 i6 9.080000E-01 5.000000E+01 i7 8.810000E-01 6.000000E+01 i8 8.500000E-01 7.000000E+01 i9 8.180000E-01 8.000000E+01 i10 7.840000E-01 9.000000E+01 i11 7.510000E-01 1.000000E+02 i12 7.180000E-01 1.100000E+02 i13 6.850000E-01 1.200000E+02 i14 6.580000E-01 1.300000E+02 i15 6.280000E-01 1.400000E+02 i16 6.030000E-01 1.500000E+02 i17 5.800000E-01 1.600000E+02 i18 5.580000E-01 1.700000E+02 i19 5.380000E-01 1.800000E+02 i20 5.220000E-01 1.900000E+02 i21 5.060000E-01 2.000000E+02 i22 4.900000E-01 2.100000E+02 i23 4.780000E-01 2.200000E+02 i24 4.670000E-01 2.300000E+02 i25 4.570000E-01 2.400000E+02 i26 4.480000E-01 2.500000E+02 i27 4.380000E-01 2.600000E+02 i28 4.310000E-01 2.700000E+02 i29 4.240000E-01 2.800000E+02 i30 4.200000E-01 2.900000E+02 i31 4.140000E-01 3.000000E+02 i32 4.110000E-01 3.100000E+02 i33 4.060000E-01 3.200000E+02 ; * * 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' / 3.7541005211E-01 / cb2 'certified value for b2' / 1.9358469127E+00 / cb3 'certified value for b3' / -1.4646871366E+00 / cb4 'certified value for b4' / 1.2867534640E-02 / cb5 'certified value for b5' / 2.2122699662E-02 / ce1 'certified std err for b1 ' / 2.0723153551E-03 / ce2 'certified std err for b2 ' / 2.2031669222E-01 / ce3 'certified std err for b3 ' / 2.2175707739E-01 / ce4 'certified std err for b4 ' / 4.4861358114E-04 / ce5 'certified std err for b5 ' / 8.9471996575E-04 / ; *----------------------------------------------------------------------------- * 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' r(i) 'residuals' ; equations fit1(i) 'the non-linear model, LS format' obj1 'objective, LS format' fit2(i) 'the non-linear model, NLP format' obj2 'objective, NLP format' ; obj1.. sse =n= 0; fit1(i).. y(i) =e= b1 + b2*exp[-x(i)*b4] + b3*exp[-x(i)*b5]; obj2.. sse =e= sum(i, sqr(r(i))); fit2(i).. y(i) =e= b1 + b2*exp[-x(i)*b4] + b3*exp[-x(i)*b5] + r(i); model nlfit1 /obj1,fit1/; model nlfit2 /obj2,fit2/; *----------------------------------------------------------------------------- * first set of initial values *----------------------------------------------------------------------------- b1.l = 50; b2.l = 150; b3.l = -100; b4.l = 1; b5.l = 2; * The starting point is bad. * NL2SOL will fail with "False Convergence". * We solve this by first using CONOPT to find an * optimal solution which we then pass on to NL2SOL. option nlp=conopt; solve nlfit2 minimizing sse using nlp; option nlp=nls; solve nlfit2 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 = 0.5; b2.l = 1.5; b3.l = -1; b4.l = 0.01; b5.l = 0.02; solve nlfit1 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";